THE STABILITY OF TIDALLY DEFORMED NEUTRON STARS TO THREE- AND FOUR-MODE COUPLING

Consider a two-dimensional harmonic oscillator with characteristic frequencies ω₁ and ω₂, such that ω₁ ≫ ω₂. The potential energy of this system for a general displacement η from equilibrium, written in coordinates with a unit-mass normalized kinetic energy term ${\mathcal {T}}=({1}/{2}) \dot{\eta }^\top \dot{\eta }$, is given by

$$\begin{eqnarray} \mathcal {V}(\eta) = \frac{1}{2} \eta ^{\top} \mathcal {M} \eta, \quad {\rm where }\; \mathcal {M} ={\left(\begin{array}{@{} cc @{}}\omega _1^2 & 0 \\ 0 & \omega _2^2 \end{array}\right)}. \end{eqnarray} \tag{ 1 }$$

Consider some effect (like interaction with a third degree of freedom, for example) which rotates the coordinate basis, so that displacement vectors become η' = U(θ)η, where U(θ) is in the SO(2) representation of the rotation. In the new basis, the potential energy is given by

$$\begin{eqnarray} &&\mathcal {V}(\eta ^\prime) = \frac{1}{2} (\eta ^\prime)^\top \mathcal {M}^\prime \eta ^\prime, \quad{\rm where }\; \mathcal {M}^\prime = U \mathcal {M} U^\top. \end{eqnarray} \tag{ 2 }$$

Putting in the form of U(θ), we see that

$$\begin{eqnarray} \mathcal {M}^\prime (\theta) &= &{\left(\begin{array}{@{}cc@{}}\omega _1^2 \cos ^2{\theta } + \omega _2^2 \sin ^2{\theta } & \left(\omega _1^2 - \omega _2^2 \right) \cos {\theta } \sin {\theta } \\ \left(\omega _1^2 - \omega _2^2 \right) \cos {\theta } \sin {\theta } & \omega _1^2 \sin ^2{\theta } + \omega _2^2 \cos ^2{\theta } \end{array}\right)} \nonumber\\ &=& \mathcal {M} + \delta \mathcal {M}, \quad{\rm where} \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \delta \mathcal {M} &=& {\left(\begin{array}{@{}cc@{}}\left(\omega _2^2 - \omega _1^2 \right) \theta ^2 &\left(\omega _1^2 - \omega _2^2 \right) \theta \\ \left(\omega _1^2 - \omega _2^2 \right) \theta & \left(\omega _1^2 - \omega _2^2 \right) \theta ^2 \end{array}\right)} + {\cal O}(\theta ^3). \end{eqnarray} \tag{ 4 }$$

The change in the smaller eigenvalue of $\mathcal {M}$ due to this change can be formally calculated using second-order perturbation theory as

$$\begin{eqnarray} \omega _-^2 &=& \, \omega _2^2 + \delta \mathcal {M}_{22} + \frac{\delta \mathcal {M}_{12}^2}{\omega _2^2 - \omega _1^2} + \cdots = \omega _2^2 + \theta ^2 \left(\omega _1^2 - \omega _2^2 \right) \nonumber\\ && + \frac{\left[\left(\omega _1^2 - \omega _2^2 \right) \theta \right]^2}{\omega _2^2 - \omega _1^2} + \cdots = \omega _2^2. \end{eqnarray} \tag{ 5 }$$

We could have predicted this from the fact that the interactions act as a pure rotation, but this avenue of analysis highlights a fact which is useful when we do not have such global, non-perturbative information. When calculating perturbations to the eigenvalues, perturbations to the matrix elements along the diagonal ($\mathcal {M}_{22}$) enter into the analysis at the same order as perturbations to those off the diagonal ($\mathcal {M}_{12}$) squared. Hence, if the angle θ is a small angle such that $\omega _2^2/(\omega _1^2 - \omega _2^2) \simeq \omega _2^2/\omega _1^2 < \theta ^2 \ll 1$, the diagonal perturbation $\delta \mathcal {M}_{22}$ is much smaller than the off-diagonal perturbation $\delta \mathcal {M}_{12}$, but the latter cannot be ignored in the calculation of ω₋. Ignoring it would lead us to the conclusion that the deformed potential has an unstable direction, i.e., $\omega _-^2 = \omega _2^2 + \theta ^2 (\omega _2^2 - \omega _1^2) + \cdots < 0$. Figure 1 illustrates this point by plotting contours of constant potential $\mathcal {V}$ using both the full rotation, and using only the leading-order term in small angle θ. In the former case, it is clear that the origin remains a stable equilibrium, but in the latter case, neglecting the higher order terms leads the origin to become a saddle point.

Zoom In Zoom Out Reset image size

This example captures much of the physics describing the system we are interested in—a tidally deformed neutron star. In Section 2.4, we show that in the sub-matrix of a pair of high-order modes, the off-diagonal perturbation to the potential is given by the three-mode coupling between the modes and the tidal deformation, and the diagonal perturbations are given by appropriate four-mode couplings. Lessons learned from the toy model tell us that we need to evaluate both the three- and the four-mode couplings in order to determine the lowest order or ² (in dimensionless tidal strength ) perturbation to the eigenfrequencies.

The example also captures one other aspect of our analysis—what happens if the shallow direction of the potential changes with time. In the limit of ω₁ → ∞, this example corresponds to that of a bead in a restoring force ($-\omega _2^2x$) sliding on a wire at position angle θ. If θ varies, then the bead experiences a centrifugal or anti-restoring force ($\dot{\theta }^2x$). In Appendix E we consider the consequences of this effect when the tidal field is varying with time due to the motion of the binary.

The remainder of this paper is structured as follows. We review the theory of nonlinear perturbations of a star using variational techniques in Section 2. We first discuss the Lagrangian formulation of the dynamics in Section 2.1, review the expansion of the problem in the basis of the linear modes of the star in Section 2.2, and discuss the equilibrium tidal deformation in Section 2.3. We arrive at the full expression of the perturbations to the mode frequencies due to three- and four-mode coupling in Section 2.4. In Section 3 we compute the four-mode coupling terms using a novel technique. We describe this technique in Sections 3.1 and 3.2, apply it to recast the problem of perturbations of the mode frequencies in Sections 3.3 and 3.4, and show that largest potentially unstable terms cancel in Section 3.6. In Section 4 we estimate the size of the remaining perturbations, and show that they are small in the case of high-order daughter modes in the presence of the tidal deformation. Finally, we conclude with an overview and discussion of our results in Section 5. Some technical details that arise along the way are collected into the Appendices.

Consider a non-rotating, inviscid, fluid neutron star with total mass M, radius R_*, and a characteristic dynamical frequency $\omega _0^2 = G M/R_*^3$, which is perturbed by a weak tidal field. The tidal potential has the form U(x), where is the dimensionless tidal strength. Our particular expression for U is the leading-order tidal potential due to a distant companion star of mass m, held at a fixed separation a, so that $\epsilon = G m/(\omega _0^2 a^3)$. We consider the leading multipole of the corresponding tidal potential, which is

$$\begin{eqnarray} &&U = - \omega _0^2 r^2 P_2 (\cos \theta), \end{eqnarray} \tag{ 6 }$$

where P_l is a Legendre polynomial with unit normalization, and θ measures the angle from the line joining the star to its companion. The perturbing potential is axisymmetric in this coordinate system. This work can be generalized to include the higher, weaker multipoles of the potential in a straightforward manner.

Let ξ denote a general displacement field on the star. Each fluid element of the star is labeled by its initial coordinates, x (the Lagrangian coordinates of the element), so that elements initially at x are displaced to their actual (Eulerian) positions x' as x → x' = x + ξ. Whatever coordinate system we choose to use, whether it be x, x', or some other coordinate system X, note that the masses of the elements are invariant, so that d³xρ(x) = d³x'ρ'(x') = d³Xρ_X(X). From now on, quantities such as ρ are taken to be defined in the coordinate system being used at that stage.

The displacement field has a Lagrangian $\mathcal {L}$ given by

$$\begin{equation} \mathcal {L}({\boldsymbol \xi },\dot{\boldsymbol \xi })=\int d^3 x \, \rho ({\boldsymbol x}) \, \frac{1}{2} \dot{\boldsymbol \xi }({\boldsymbol x})\cdot \dot{\boldsymbol \xi }({\boldsymbol x}) - \mathcal {V} ({\boldsymbol \xi }), \end{equation} \tag{ 7 }$$

where the potential $\mathcal {V}(\xi)$ incorporates both the internal energy of the perturbed star and the gravitational energy. For a star with an external perturbing field, the potential energy with an arbitrary displacement field takes the form

$$\begin{eqnarray} \mathcal {V}({\boldsymbol \xi })& =& \int d^3 x^{\prime } \, \rho ({\boldsymbol x}^{\prime }) [ \mathcal {E}_{\rm int}({\boldsymbol x}^{\prime }) + \Phi _0({\boldsymbol x}^{\prime }) + \epsilon U({\boldsymbol x}^{\prime })] + \mathcal {C} \nonumber \\ &=& \int d^3 x \, \rho ({\boldsymbol x}) [ \mathcal {E}_{\rm int}({\boldsymbol x} + {\boldsymbol \xi }) + \Phi _0({\boldsymbol x} + {\boldsymbol \xi }) + \epsilon U({\boldsymbol x} + {\boldsymbol \xi })] + \mathcal {C} \nonumber \\ &=& \int d^3 x \, \rho ({\boldsymbol x}) \left[ \mathcal {E}_{\rm int}({\boldsymbol x}) + \Phi _0 ({\boldsymbol x}) + \frac{1}{2} {\boldsymbol \xi } \cdot {\boldsymbol C} \cdot {\boldsymbol \xi } \right. \nonumber\\ &&\left. + \,\frac{1}{3!}f_3({\boldsymbol \xi },{\boldsymbol \xi },{\boldsymbol \xi })+\frac{1}{4!}f_4({\boldsymbol \xi },{\boldsymbol \xi },{\boldsymbol \xi },{\boldsymbol \xi })+\cdots \right. \nonumber \\ &&\left. +\, \epsilon \, {\boldsymbol \xi } \cdot {\boldsymbol \nabla } U +\frac{1}{2} \epsilon \, {\boldsymbol \xi } \cdot ({\boldsymbol \xi } \cdot {\boldsymbol \nabla }) {\boldsymbol \nabla } U + \cdots \right] + \mathcal {C}. \end{eqnarray} \tag{ 8 }$$

In the case of linear perturbations, only the symmetric bilinear C and the gradient of U contribute to the dynamics of the displacement field. Physically, C · ξ is the linear restoring force opposing an infinitesimal displacement ξ. Lynden-Bell & Ostriker (1967) (see also, e.g., Schenk et al. 2002) derive a functional form for C. The functionals f_n, meanwhile, encode the nonlinear corrections to the restoring forces due to the internal and self-gravitational energy of the star. These functionals are symmetric and linear under addition and scalar multiplication of displacements, but not under multiplication of the displacements by scalar functions. Later, when we expand the displacement in terms of the linear modes, they give us the n-mode couplings.

The term $\mathcal {C}$ in this case contains the contributions to the energy from the gravitational field itself, which are fixed by the Cowling approximation and do not contribute to the dynamics of ξ. The issue of the appropriate division of gravitational potential energy between the field degrees of freedom and the interaction between the fluid elements and those fields is discussed in Appendix A. We also remark that for the tidal field given in Equation (6), ∫d³x ρ U = 0 due to the fact that the background ρ is spherically symmetric and U is dipolar, and so has not been written in Equation (8). In addition, since U is a quadratic function of x, all third and higher derivatives of U vanish, and as such have not been written.

The most convenient language in which to discuss nonlinear perturbations is in terms of an expansion in the orthonormal basis of linear modes of the star. We write the mode functions themselves are as ξ_a, and a general displacement can be expanded as

$$\begin{eqnarray} &&{\boldsymbol \xi } ({\boldsymbol x}) = \sum _a c_a {\boldsymbol \xi }_a ({\boldsymbol x}) \end{eqnarray} \tag{ 9 }$$

using this basis. In spherical coordinates, the mode functions have the form¹

$$\begin{eqnarray} &&{\boldsymbol \xi }_a = \xi _r Y_a \hat{\boldsymbol r} + \xi _h \left(\partial _\theta Y_a \hat{\boldsymbol \theta } + \frac{1}{\sin {\theta }} \partial _\phi Y_a \hat{\boldsymbol \phi } \right), \end{eqnarray} \tag{ 10 }$$

where ξ_r and ξ_h are functions of r and $Y_a = Y_{l_a,m_a} (\theta,\phi)$. The basis obeys the orthonormality condition

$$\begin{eqnarray} &&\int d^3 x \, \rho \, {\boldsymbol \xi }_a^* \cdot {\boldsymbol \xi }_b = \frac{E_0}{\omega _a^2} \delta _{ab}, \end{eqnarray} \tag{ 11 }$$

where E₀ = GM²/R_*. This normalization, which is the same as the one used by WAB but differs from that of Schenk et al. (2002), means that when a basis mode is excited to unit amplitude, it has energy E₀. By expanding in these basis functions, we can re-express our Lagrangian and potential as sums over mode amplitudes.

The square of the mode frequencies $\omega _a^2$ are the eigenvalues of the bilinear C in the mode basis,

$$\begin{eqnarray} &&{\boldsymbol C} \cdot {\boldsymbol \xi }_a = \omega _a^2 {\boldsymbol \xi }_a. \end{eqnarray} \tag{ 12 }$$

By defining the mode expansions of the remaining terms in the Lagrangian, we can write it entirely in terms of the mode amplitudes. The definitions we need are

$$\begin{eqnarray} U_a &=& {-} \frac{1}{E_0} \int d^3 x \, \rho \, {\boldsymbol \xi }_a^* \cdot {\boldsymbol \nabla } U, \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} U_{ab} &=& {-} \frac{1}{E_0} \int d^3 x \, \rho \, {\boldsymbol \xi }_a \cdot ({\boldsymbol \xi }_b \cdot {\boldsymbol \nabla }){\boldsymbol \nabla } U, \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} \kappa _{abc} &=& {-} \frac{1}{2E_0} \int d^3 x \, \rho \, f_3({\boldsymbol \xi }_a,{\boldsymbol \xi }_b,{\boldsymbol \xi }_c), \quad {\rm and} \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \kappa _{abcd} &=& {-} \frac{1}{3! E_0} \int d^3 x \, \rho \, f_4({\boldsymbol \xi }_a,{\boldsymbol \xi }_b,{\boldsymbol \xi }_c,{\boldsymbol \xi }_d). \end{eqnarray} \tag{ 16 }$$

Due to the symmetry of the functionals f_n, the coupling terms κ_abc and κ_abcd are symmetric in all their indices. When we expand the Lagrangian in this mode basis, we use the fact that the displacement ξ is real, and can be replaced by ξ* as and when needed. Using these definitions and the potential (8), we rewrite the Lagrangian (7) in the form

$$\begin{eqnarray} \mathcal {L} &= & E_0 \sum \left[ \frac{\dot{c}_a \dot{c}_a^{*}}{2 \omega _a^2} - \frac{1}{2} c_a c_a^{*} + \frac{1}{3} \kappa _{\bar{a} \bar{b} \bar{c}} c_a^{*} c_b^{*} c_c^{*} + \frac{1}{4} \kappa _{\bar{a} \bar{b} \bar{c} \bar{d}} c_a^{*} c_b^{*} c_c^{*} c_d^{*} \right. \nonumber\\ &&+ \left. \epsilon U_a c_a^{*} + \frac{1}{2} \epsilon U_{\bar{a} \bar{b}} c_a^{*} c_b^{*} + \cdots \right] - \mathcal {V}(0). \end{eqnarray} \tag{ 17 }$$

Here, our convention is that the sum runs over all the repeated indices in each term (including the index a in the first two terms). Bars indicate the use of complex conjugate wave-functions in the respective terms. Note that in the first and second terms, the amplitude of a mode with a given m comes up twice—in the terms corresponding to m and in those corresponding to −m, since the reality of ξ implies that they are related by $c_{a_{-m}}=(-1)^m c_{a_m}^\ast$. We recall that the potential term $\mathcal {V}$ with the displacement set to zero has the form

$$\begin{eqnarray} \mathcal {V}(0) = \int d^3 x \, \rho &\left[ \mathcal {E}_{\rm int}({\boldsymbol x}) + \Phi _0 ({\boldsymbol x}) \right] + \mathcal {C}. \end{eqnarray} \tag{ 18 }$$

We use the rule that sums run over repeated indices throughout this work, except where noted. Though we do not do so here, the equations of motion, as presented in, e.g., Schenk et al. (2002), can be derived by varying this Lagrangian with respect to the amplitude $c_a^*$.

With the formalism of the mode expansion established, we now investigate the perturbations excited by a static tidal field.

The star responds to the external tidal field U by deforming to a new equilibrium configuration. We denote this equilibrium tidal deformation χ and set ξ equal to χ in the preceding equations. The tidal displacement χ is such the internal restoring forces balance the external perturbing force,

$$\begin{eqnarray} &&\left. \frac{\delta \mathcal {V}}{\delta {\boldsymbol \xi }}\right|_{{\boldsymbol \xi } \rightarrow {\boldsymbol \chi }} = 0. \end{eqnarray} \tag{ 19 }$$

Physically, χ takes the spherically symmetric star to a static configuration where elements along contours of constant gravitational potential Φ₀ + U have the same density and pressure. The displacement χ has a part linear in the tidal strength , which is the so-called linear tide. There are also higher terms, which arise from the nonlinear restoring forces and the nonlinear tide. In order to consistently account for the changes in the eigenfrequencies ω_a of additional perturbations due to coupling with the tide, we need to keep terms up to O(²). As such, we write the expansion of χ in the mode basis as

$$\begin{eqnarray} &&{\boldsymbol \chi } =\sum \chi _a {\boldsymbol \xi }_a = \sum \left(\epsilon \chi ^{(1)}_a + \epsilon ^2 \chi ^{(2)}_a \right) {\boldsymbol \xi }_a. \end{eqnarray} \tag{ 20 }$$

Varying just the potential terms of the Lagrangian (17) with respect to a given mode amplitude $c_a^*$, and then substituting the amplitudes χ_a into the result gives an equation for χ_a,

$$\begin{eqnarray} &&\chi _a- \epsilon U_a -\! \sum [\kappa _{\bar{a} \bar{b} \bar{c}} \chi _b^{*}\chi _c^{*} + \kappa _{\bar{a} \bar{b} \bar{c} \bar{d}} \chi _b^{*}\chi _c^{*}\chi _d^{*} + \epsilon U_{\bar{a} \bar{b}} \chi _b^{*} +\cdots] = 0. \nonumber\\ \end{eqnarray} \tag{ 21 }$$

Solving order by order for χ_a, we find that

$$\begin{eqnarray} &&\chi ^{(1)}_a = U_a \quad {\rm and} \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} &&\chi ^{(2)}_a = \sum [ \kappa _{\bar{a} \bar{b} \bar{c}} U_b^{*} U_c^{*} + U_{\bar{a} \bar{b}} U_b^{*}]. \end{eqnarray} \tag{ 23 }$$

Given the tidal field U, the fluid elements respond by the nonlinear displacement χ. Note that since U is axisymmetric, χ is also axisymmetric and χ_a contains only terms where m = 0. This means that all of the terms in Equations (22) and (23) are actually real. The question we must answer is whether or not this deformed star is stable.

Now we allow the star to undergo further perturbations, so that

$$\begin{eqnarray} &&{\boldsymbol \xi } = {\boldsymbol \chi } + {\boldsymbol \eta }, \end{eqnarray} \tag{ 24 }$$

where η is an additional displacement field on the star away from the equilibrium configuration. We are interested in the interaction of pairs of daughter modes with the tidal perturbation, and so we expand our Lagrangian only to second order in the additional small perturbations, O(η²). We do not deal with the coupling of three or more non-tidal modes either to each other or to the tide, and so O(η³) and higher terms are dropped in what follows.

By expanding η as $\sum \eta _a {\boldsymbol \xi }_a$, and substituting the corresponding expansion of the displacement c_a = χ_a + η_a into the Lagrangian (17), we arrive at the following form of the Lagrangian, up to the order we are interested in,

$$\begin{eqnarray} \mathcal {L} &=& E_0 \sum \left[ \frac{\dot{\eta }_a \dot{\eta }_a^{*}}{2 \omega _a^2} - \frac{1}{2} (\chi _a \chi _a^{*} + \eta _a \eta _a^{*}) + \kappa _{\bar{a} \bar{b} \bar{c}} \eta _a^{*} \eta _b^{*} \chi _c^{*} \right. \nonumber\\ &&+ \left.\frac{3}{2} \kappa _{\bar{a} \bar{b} \bar{c} \bar{d}} \eta _a^{*} \eta _b^{*} \chi _c^{*} \chi _d^{*} + \epsilon U_a \chi _a^{*} + \frac{1}{2} \epsilon U_{\bar{a} \bar{b}} \eta _a^{*} \eta _b^{*} \right] - \mathcal {V}(0) \nonumber \\ &=& E_0 \sum \left[ \frac{\dot{\eta }_a \dot{\eta }_a^{*}}{2 \omega _a^2} - \frac{1}{2} \eta _a \eta _a^{*} + \left(\epsilon \left[\kappa _{\bar{a} \bar{b} \bar{c}} \chi _c^{(1)*} + \frac{1}{2} U_{\bar{a} \bar{b}} \right] \right.\right. \nonumber\\ &&+ \left.\left. \epsilon ^2 \left[\kappa _{\bar{a} \bar{b} \bar{c}} \chi _c^{(2)*} + \frac{3}{2} \kappa _{\bar{a} \bar{b} \bar{c} \bar{d}} \chi _c^{(1)*} \chi _d^{(1)*} \right] \right) \eta _a^{*} \eta _b^{*} \right] \nonumber\\ &&+ {\frac{E_0}{2 }\sum \epsilon ^2 |U_a|^2} - \mathcal {V}(0). \end{eqnarray} \tag{ 25 }$$

In the first line we have written only those terms which are of the correct order, and we have also neglected terms which have only a single factor of η_a, since these vanish upon substitution of our solution for χ_a. The terms linear in η vanish because the displacement χ takes the star to an equilibrium state. In the next line, we have substituted the decomposition (20), collected terms in orders of , and used Equation (22) to isolate the |U_a|² term. For now, we ignore the overall constant terms $\mathcal {V}(0)$ and E₀∑²|U_a|²/2, which follow along for the ride but do not affect the dynamics.

The O() term in Equation (25) is the three-mode interaction $\kappa _{\bar{a} \bar{b} \bar{c}} \eta _a^* \eta _b^* \chi _c^{(1)*}$. If the coupling $\kappa _{\bar{p} \bar{g} \bar{c}} \chi _c^{(1)*}$ happens to be large and positive for a particular pair of modes (p, g) then it can overcome the smallness of the tidal coupling strength . As first noted by WAB, $\kappa _{\bar{p} \bar{g} \bar{c}} \chi _c^{(1)*}$ is in fact large for certain high-order p-mode and g-mode pairs, due to a spatial resonance of the mode functions ξ_p and ξ_g. However, as we show below, this three-mode term perturbs the characteristic frequencies at O(²), as do all of the terms in Equation (25) multiplied by ². It is not immediately clear what role these terms play.

By varying the Lagrangian (25) with respect to the amplitudes η_a we can derive the equations of motion, and from there the characteristic frequencies. We instead adopt an equivalent, but hopefully more transparent strategy to obtain the characteristic frequencies. Defining rescaled amplitudes $\eta _a ^{\prime } = \eta _a/\omega _a$ leads us to the analogy to a system of coupled oscillators, for which we can rewrite the Lagrangian in matrix notation as

$$\begin{eqnarray} &&\mathcal {L} = \frac{E_0}{2} (\dot{\eta }^{\prime })^\dagger \dot{\eta }^{\prime } - \frac{E_0}{2} (\eta ^{\prime })^\dagger \mathcal {M} \eta ^{\prime }, \end{eqnarray} \tag{ 26 }$$

where η' is a vector of mode amplitudes and the matrix $\mathcal {M}$ contains the leading restoring terms, the effect of the tidal potential, and the mode–mode interactions. The rescaling introduces factors of ω_a at each instance of an η_a, and it is equivalent to normalizing the mode functions to all have the same moment of inertia $M R_*^2$ at unit amplitude rather than the same energy E₀, which is the choice of normalization used in Schenk et al. (2002).

The eigenvalues of $\mathcal {M}$ are the characteristic frequencies of this set of oscillators, and consideration of $\mathcal {M}$ is equivalent to writing out the equations of motion and solving the corresponding characteristic equation. Because the tidal potential is axisymmetric, there is an ordering of the basis modes where the matrix $\mathcal {M}$ is made up of 2 × 2 and 1 × 1 blocks along the diagonal, where each block consists of those modes which are allowed to interact given conservation of angular momentum. The blocks are all independent of each other, and analysis of one essentially applies to the rest.

We now focus on a particular pair of modes with m = 0, which we give indices p and g, and write the sub-block of $\mathcal {M}$ for this pair. We discuss the structure of $\mathcal {M}$ and the case of modes where m ≠ 0 further in Appendix B. The modes have unperturbed frequencies ω_p and ω_g, and we set ω_p > ω_g. The sub-block of $\mathcal {M}$ is composed of

$$\begin{eqnarray} \mathcal {M}_{pp} &= & \omega _p^2 - \epsilon \omega _p^2 \left(U_{pp} + \sum 2 \kappa _{app} \chi _a^{(1)} \right) \nonumber\\ &&- \epsilon ^2 \omega _p^2 \sum \left[ 2 \kappa _{app} \chi _a^{(2)} + 3 \kappa _{abpp} \chi _a^{(1)} \chi _b^{(1)} \right], \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} \mathcal {M}_{pg} &= & \mathcal {M}_{gp} = -\epsilon \, \omega _p \omega _g \left(U_{pg} + \sum 2 \kappa _{apg} \chi _a^{(1)} \right), \quad {\rm and} \nonumber\\ \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} \mathcal {M}_{gg} &= & \omega _g^2 - \epsilon \omega _g^2 \left(U_{gg} + \sum 2 \kappa _{agg} \chi _a^{(1)} \right) \nonumber\\ &&- \epsilon ^2 \omega _g^2 \sum \left[ 2 \kappa _{agg} \chi _a^{(2)} + 3 \kappa _{abgg} \chi _a^{(1)} \chi _b^{(1)} \right], \end{eqnarray} \tag{ 29 }$$

and note that since p and g denote particular modes, they are not summed over in the expressions for $\mathcal {M}$. Formally expanding in powers of , we can compute the perturbed eigenvalues of this matrix as

$$\begin{eqnarray} \omega _+^2 &= &\, \omega _p^2 - \epsilon \omega _p^2 \left(U_{pp} + \sum 2 \kappa _{app} \chi _a^{(1)} \right)\nonumber\\ && - \epsilon ^2 \omega _p^2 \sum \left[2 \kappa _{app} \chi _a^{(2)} + 3 \kappa _{abpp} \chi _a^{(1)} \chi _b^{(1)} \right] \nonumber\\ && + \epsilon ^2 \frac{ \omega _p^2 \omega _g^2 }{\omega _p^2 - \omega _g^2}\left(U_{pg} + \sum 2 \kappa _{apg} \chi _a^{(1)} \right)^2 \end{eqnarray} \tag{ 30 }$$

and

$$\begin{eqnarray} \omega _-^2 &= &\, \omega _g^2 - \epsilon \omega _g^2 \left(U_{gg} + \sum 2 \kappa _{agg} \chi _a^{(1)} \right) \nonumber\\ &&- \epsilon ^2 \omega _g^2 \sum \left[2 \kappa _{agg} \chi _a^{(2)} + 3 \kappa _{abgg} \chi _a^{(1)} \chi _b^{(1)} \right] \nonumber\\ && - \epsilon ^2 \frac{ \omega _p^2 \omega _g^2 }{\omega _p^2 - \omega _g^2} \left(U_{pg} + \sum 2 \kappa _{apg} \chi _a^{(1)} \right)^2. \end{eqnarray} \tag{ 31 }$$

Since ω_g is the smaller frequency to begin with, the negative-definite term involving κ_apg in Equation (31) is in danger of pushing ω_g to a negative value, especially if the mode is a high-order g-mode with $\omega _g^2 \ll \omega _0^2 = G M/R_*^3 \ll \omega _p^2$. This is the potential nonlinear instability described by WAB.

In order to make definite statements, we need to know how the other terms behave, in particular the four-mode interaction term $\kappa _{abgg} \chi _a^{(1)} \chi _b^{(1)}$. If they are large when the three-mode coupling is large, they can in principle prevent the instability. We now discuss a method for computing the interaction between two daughter modes and the tidal deformation at the level of the four-mode coupling. Along the way we also find a way to simply derive the three-mode terms. We find that the four-mode term serves to precisely cancel the κ_apg term.

To begin with, we apply the tidal perturbation U to a spherical star, so that the total gravitational potential becomes Φ₀ + U. In response the fluid elements of the star are displaced by the field χ as discussed in Section 2.3. Next, we consider a coordinate transform $\psi: {\boldsymbol x} = (r, \theta, \phi) \rightarrow {\boldsymbol X} = (R, \Theta, \phi)$. We require this transform to have the special properties that it is volume-preserving, and that it returns the star to spherical symmetry. These properties completely determine our mapping in terms of a power series expansion in . The resulting spherically symmetric star is not the original star, because χ itself is not volume-preserving; but as we will see the coordinate transform reverses χ at leading order in . Since our problem is axisymmetric, the ϕ coordinate is unchanged, and from here on we can safely ignore the coordinate ϕ.

We require an expansion of the coordinate transform in orders in , and so we represent ψ by a coordinate flow under an infinitesimal generator ζ(x). Explicitly, the finite transformation ψ generated by the infinitesimal transformation ζ is defined via a parameterized transform φ(s, x) by the ordinary differential equation $d\varphi (s,{\boldsymbol x})/ds = {\boldsymbol \zeta }(\varphi (s,{\boldsymbol x}))$, initial condition φ(0, x) = x, and assignment X = ψ(x) ≡ φ(1, x). The flow forward, ψ, gives us our new coordinates in the form of a Taylor expansion

$$\begin{eqnarray} &&{\boldsymbol X} = {\boldsymbol x}+ {\boldsymbol \zeta }({\boldsymbol x})+ \left. \frac{1}{2} ({\boldsymbol \zeta } \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta } \right|_{\boldsymbol x}, \end{eqnarray} \tag{ 32 }$$

up to order ζ². The volume-preserving requirement is equivalent to requiring the generator to be divergenceless, ∇ · ζ = 0.² It is actually more convenient to first get an explicit representation of ψ⁻¹: this is actually the transformation generated by the negative of the generator. (This can be seen from the integral curve definition of a generator if we reverse the sign of d/ds; for an explicit proof see Appendix C.) Thus the mapping from X to x via −ζ is

$$\begin{eqnarray} &&{\boldsymbol x} = {\boldsymbol X} + (-{\boldsymbol \zeta })({\boldsymbol X}) +\left. \frac{1}{2} [(-{\boldsymbol \zeta }) \cdot {\boldsymbol \nabla }](- {\boldsymbol \zeta }) \right|_{\boldsymbol X}. \end{eqnarray} \tag{ 33 }$$

Because of this, once we find an expression for the generator of the inverse flow −ζ in (R, Θ) coordinates, we have the desired generator ζ(x).³

The vector −ζ(X) has an expansion in powers of , which is

$$\begin{eqnarray} &&{\boldsymbol \zeta } = \epsilon {\boldsymbol \zeta }^{(1)} + \epsilon ^2 {\boldsymbol \zeta }^{(2)} + O (\epsilon ^3). \end{eqnarray} \tag{ 34 }$$

From Equation (33), the magnitude r of the coordinate vector x is

$$\begin{eqnarray} r &=& R- \epsilon \hat{\boldsymbol R} \cdot {\boldsymbol \zeta }^{(1)} +\epsilon ^2 \left[ \frac{ {\boldsymbol \zeta }^{(1)} \cdot {\boldsymbol \zeta }^{(1)}}{2 R} - \frac{ (\hat{\boldsymbol R} \cdot {\boldsymbol \zeta }^{(1)})^2}{2R} \right. \nonumber\\ &&+ \left.{\frac{\hat{\boldsymbol R} \cdot ({\boldsymbol \zeta }^{(1)} \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta }^{(1)}}{2}} - \hat{\boldsymbol R} \cdot {\boldsymbol \zeta }^{(2)}\right], \end{eqnarray} \tag{ 35 }$$

which can be further simplified using the vector identity $({\boldsymbol A} \cdot {\boldsymbol \nabla }) {\boldsymbol A} = {\boldsymbol \nabla } A^2 / 2- {\boldsymbol A} \times ({\boldsymbol \nabla } \times {\boldsymbol A})$, to yield

$$\begin{eqnarray} &&r = R - \epsilon \zeta ^{(1)}_R +\epsilon ^2 \left[ {\frac{ \zeta ^{(1)}_R \partial _R \zeta ^{(1)}_R}{2}} + {\frac{ \zeta ^{(1)}_\Theta \partial _\Theta \zeta ^{(1)}_R}{2R}} - \zeta ^{(2)}_R \right], \quad \end{eqnarray} \tag{ 36 }$$

with the definitions ${\boldsymbol \zeta }^{(i)} \cdot \hat{\boldsymbol R} = \zeta ^{(i)}_R$ and ${\boldsymbol \zeta }^{(i)} \cdot \hat{\boldsymbol \Theta } = \zeta ^{(i)}_\Theta$. Also useful is the expansion

$$\begin{eqnarray} P_2(\cos \theta) &=& (1 - \epsilon {\boldsymbol \zeta } \cdot {\boldsymbol \nabla }) P_2 (\cos \Theta) = P_2 (\cos \Theta) \nonumber\\ && - \epsilon \frac{\zeta ^{(1)}_\Theta \partial _\Theta P_2(\cos \Theta)}{R} \end{eqnarray} \tag{ 37 }$$

to the order needed. Now we are in place to solve for ζ, and at the same time to determine the form of the total gravitational potential Φ(R) out to order ². The tidally perturbed potential is

$$\begin{eqnarray} && \Phi _0(r) - \epsilon \omega _0^2 r^2 P_2 (\cos \theta) = \Phi _0 (R) - \epsilon \left[\zeta ^{(1)}_R \mathfrak {g} + \omega _0^2 R^2 P_2(\cos \Theta) \right] \nonumber\\ &&\quad +\, \epsilon ^2 \left[ \left({\frac{\zeta ^{(1)}_R \partial _R \zeta ^{(1)}_R}{2}} + {\frac{ \zeta ^{(1)}_\Theta \partial _\Theta \zeta ^{(1)}_R}{2R}} - \zeta ^{(2)}_R \right) \mathfrak {g} \right. \nonumber \\ &&\quad \left. + \frac{1}{2} \left(\zeta ^{(1)}_R \right)^2 \frac{d \mathfrak {g}}{dr} + 2 \omega _0^2 R \zeta ^{(1)}_R P_2(\cos \Theta) + \omega _0^2 R \zeta ^{(1)}_\Theta \partial _\Theta P_2(\cos \Theta) \right]\!. \nonumber\\ \end{eqnarray} \tag{ 38 }$$

For convenience here and later, we have denoted the local gravitational acceleration as $\mathfrak {g} = d \Phi _0/{dR}$.

By insisting that all Θ dependence in Φ vanishes, we find at first order

$$\begin{eqnarray} &&\zeta ^{(1)}_R = - \omega _0^2 \frac{R^2 }{\mathfrak {g}}P_2 (\cos \Theta), \end{eqnarray} \tag{ 39 }$$

and at second order

$$\begin{eqnarray} \zeta ^{(2)}_R & = &{\frac{\zeta ^{(1)}_R \partial _R \zeta ^{(1)}_R}{2}} + {\frac{ \zeta ^{(1)}_\Theta \partial _\Theta \zeta ^{(1)}_R}{2R}} + \frac{ \left(\zeta ^{(1)}_R \right)^2 }{2R} \frac{d \ln \mathfrak {g}}{d \ln R}\nonumber\\ && + \frac{\omega _0^2}{\mathfrak {g}}\left[ 2 R \zeta ^{(1)}_R P_2(\cos \Theta) + R\zeta ^{(1)}_\Theta \partial _\Theta P_2(\cos \Theta)\right]. \quad \end{eqnarray} \tag{ 40 }$$

In order to write this last equation entirely in terms of the gravitational potentials, we need an expression for $\zeta _\Theta ^{(1)}$. For this, we use the requirement that ψ be volume-preserving. The volume-preserving infinitesimal displacements are spanned by the vectors

$$\begin{eqnarray} {\boldsymbol f}^1_{lm} &=& l(l+1) \frac{u_{lm}}{R} \hat{\boldsymbol R} + \left(\frac{u_{lm}}{R} + \partial _R u_{lm} \right) R {\boldsymbol \nabla } Y_{lm} \quad{\rm and} \quad \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} {\boldsymbol f}^2_{lm} & =& w_{lm} {\boldsymbol R} \times {\boldsymbol \nabla } Y_{lm}. \end{eqnarray} \tag{ 42 }$$

This can be seen by considering the usual decomposition of vectors into vector spherical harmonics, and finding those combinations which satisfy ${\boldsymbol \nabla } \cdot {\boldsymbol f}_{lm} = 0$. These naturally decompose into the spheroidal displacements ${\boldsymbol f}_{lm}^1$ and the toroidal displacements ${\boldsymbol f}_{lm}^2$ (see also Blandford & Thorne 2011, chapter 12). We are interested in axially symmetric perturbations, so we only need to use the basis vectors ${\boldsymbol f}^1_{lm}$. We further absorb the normalization of the spherical harmonics into the definition of the u_lm and write (with m = 0 implicit)

$$\begin{eqnarray} {\boldsymbol \zeta } &=& \sum _l l (l+1) \frac{u_l}{R} P_l (\cos \Theta) \hat{\boldsymbol R} \nonumber\\ &&+ \left(\frac{u_l}{R} + \partial _R u_l \right) [\partial _\Theta P_l (\cos \Theta)] \hat{\boldsymbol \Theta }, \quad {\rm where} \end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray} u_{l} & =& \epsilon u^{(1)}_{l} + \epsilon ^2 u^{(2)}_{l} + O(\epsilon ^3). \end{eqnarray} \tag{ 44 }$$

We see that the ζ_Θ is simply related to ζ_R through the requirement that ζ be volume-preserving. Comparing Equation (39) to Equation (43), we immediately find that

$$\begin{eqnarray} &&u^{(1)}_l = - \omega _0^2 \frac{R^3 }{6 \mathfrak {g}} \delta _{l2}. \end{eqnarray} \tag{ 45 }$$

Substituting this into Equations (40) and (43) allows us to solve for the functions $u^{(2)}_l$. The angular terms are of the form (P₂)² and $(\partial _\Theta P_2)^2$, which when re-expanded in terms of Legendre polynomials couple only to l = 2, 4. We can pick off each of these terms by integrating Equation (40) with the appropriate Legendre polynomial and weight dμ = d[cos Θ]. Orthogonality then gives

$$\begin{eqnarray} &&u^{(2)}_2 = - \frac{\omega _0^4}{84}\frac{R^4 (8 - n)}{ \mathfrak {g}^2} \quad {\rm and}\; u^{(2)}_4 = \frac{3\omega _0^4}{350} \frac{ R^4(1- n)}{\mathfrak {g}^2}; \end{eqnarray} \tag{ 46 }$$

all the other $u^{(2)}_l$ vanish. We have defined $n = d \ln \mathfrak {g} / d\ln R$, which proves to be a convenient short hand in the estimates in Section 4. When $\mathfrak {g}$ has a simple power-law dependence on R, n is just its constant power law index. Note that the l = 0 contribution, which is Θ independent, represents the correction to the spherically symmetric potential Φ₀(R). Because of this we also have from the matching that

$$\begin{eqnarray} \Phi (R) & =& \Phi _0 (R) + \epsilon ^2 V (R) + O(\epsilon ^3),\quad {\rm where} \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} V (R)& =& {-}\frac{\omega _0^4}{10} \frac{R^3 (6-n) }{\mathfrak {g}}. \end{eqnarray} \tag{ 48 }$$

In order to have the generator of the transform ψ we simply need to replace the coordinates (R, Θ) with (r, θ), since we have already accounted for the sign change in reversing the flow. As such, we have completed the construction of ψ up to O(²).

We have an expression for the volume-preserving transform which takes the tidally perturbed star A into the spherically symmetric star B. We have seen that star B is more weakly perturbed than star A, and this allows us to match orders in perturbation in a way that gives useful relations. First though, it is useful to step back and consider the physical intuition underlying the transformation. Figure 3 illustrates the principles behind the transformation. Consider an element displaced from point $\mathcal {P}$ to point $\mathcal {Q}$ by the tidal displacement χ. As we noted before, this tidal displacement takes the star to a configuration where elements along contours of equal potential have the same density and pressure. Point ${\mathcal {Q}}$ lies on one such contour. The displacement field χ can be decomposed as the sum of certain virtual displacement fields, appropriate members of which are marked as 1, 2 and 3 in the figure. These virtual displacements are chosen to have some nice properties.

Zoom In Zoom Out Reset image size

First, each virtual displacement is chosen to be adiabatic. This is possible since their sum, χ, is adiabatic. Second, we choose 2 and 3 such that they both preserve the volume of fluid elements. The pressure and density of an element at $\mathcal {P}$ is not necessarily equal to the pressure and density of the element once it is displaced to $\mathcal {Q}$, as the tidal deformation χ is not volume-preserving beyond linear order. With this choice for 2 and 3, the entire volume change of an element, and the related pressure and density change, occurs during displacement 1. Third, the displacement 2 is chosen to be the linear part of χ (the so-called linear tide). This is possible because the linear tide is a volume-preserving deformation. Given this choice, the displacements 1 and 3 make up the nonlinear part of χ. Fourth, 3 is chosen to be that part of the nonlinear tide which preserves volumes, such that the remaining displacement, 1, takes the initial spherical potential contour that $\mathcal {P}$ lies on to another spherical contour. Hence, the displacement 1 is a radial displacement, which takes $\mathcal {P}$ to a point $\mathcal {R}$.

We are now faced with the question of what sort of perturbing fields we would need to apply to make the elements follow these virtual displacements. We have seen that 2 and 3 correspond to a pure coordinate transformation (they make up ψ⁻¹), so that the potential $\mathcal {V}$ of the element is frozen in during these displacements. The overall stellar structure naturally flows forward along with the transformation. Hence the potential at $\mathcal {R}$ is same as the potential at $\mathcal {Q}$. This differs from the background star's potential at $\mathcal {P}$, and the difference can be seen as a spherically symmetric perturbing potential causing the elements of the star to deform along 1.

We could judge this decomposition scheme solely on its power to illuminate the underlying components of the equilibrium tide χ, but its real value manifests in the study of small displacements on top of this, away from equilibrium. This is because these weak displacements, say η, which are originally on a star deformed to O(), when pulled back along 3 and 2 pick up additional corrections of O() and higher while the transformed star is more weakly deformed (perturbed at O(²)). We can explicitly compute these corrections to the displacements from the pullback. In order to get the ²η² part of the energy, we have to look at both of the functionals f₃ and f₄ in the original picture. Here we find that we only need f₃.

Now that we have an explicit form for the transformation between star A and star B, and a better understanding of the motivation behind this transform, we can consider the potential $\mathcal {V}$ of star B. Star B sits in a spherically symmetric gravitational potential, and its fluid elements have been adjusted by the combination of the displacement field χ and the coordinate transform. In total, the star has undergone a radial deformation σ (before called 1), sourced by the external potential ²V, so that ${\boldsymbol x} \rightarrow {\boldsymbol X} = {\boldsymbol x} + {\boldsymbol \sigma }$. As before, we have for the new configuration

$$\begin{eqnarray} \mathcal {V}({\boldsymbol \sigma }) &= & \int d^3 X \, \rho ({\boldsymbol X}) [ \mathcal {E}_{\rm int} ({\boldsymbol X}) + \Phi _0({\boldsymbol X}) + \epsilon ^2 V({\boldsymbol X})] \nonumber \\ & =& \int d^3 x \, \rho ({\boldsymbol x}) [ \mathcal {E}_{\rm int} ({\boldsymbol x}+{\boldsymbol \sigma }) + \Phi _0({\boldsymbol x}+{\boldsymbol \sigma }) + \epsilon ^2 V({\boldsymbol x}+{\boldsymbol \sigma })] \nonumber \\ &=& \mathcal {V}(0) + \int {d^3x} \, \rho ({\boldsymbol x}) \left[ \frac{1}{2} {\boldsymbol \sigma }\cdot {\boldsymbol C}\cdot {\boldsymbol \sigma } + \frac{1}{3!} f_3({\boldsymbol \sigma },{\boldsymbol \sigma },{\boldsymbol \sigma }) + \cdots \right.\nonumber \\ &&\left. +\, \epsilon ^2 V({\boldsymbol x}) + \epsilon ^2 {\boldsymbol \sigma }\cdot {\boldsymbol \nabla } V + \epsilon ^2 \frac{1}{2} {\boldsymbol \sigma } \cdot ({\boldsymbol \sigma } \cdot {\boldsymbol \nabla }) {\boldsymbol \nabla } V + \cdots \right]. \nonumber\\ \end{eqnarray} \tag{ 49 }$$

Expanding σ in the mode basis of the unperturbed star, ${\boldsymbol \sigma } = \sum \sigma _a {\boldsymbol \xi }_a$, we can write the potential as

$$\begin{eqnarray} \mathcal {V} &=& \mathcal {V}(0) + \int d^3 x \, \rho \, \epsilon ^2 V+ E_0 \sum \left[ \frac{1}{2} \sigma _a^2 - \frac{1}{3} \kappa _{abc} \sigma _a \sigma _b \sigma _c \right. \nonumber\\ &&\left.- \,\epsilon ^2 V_a \sigma _a - \epsilon ^2 \frac{1}{2} V_{ab} \sigma _a \sigma _b + \cdots \right], \end{eqnarray} \tag{ 50 }$$

with

$$\begin{eqnarray} V_a &=& {-}\frac{1}{E_0} \int d^3 x \, \rho \, {\boldsymbol \xi }^*_a \cdot {\boldsymbol \nabla } V \quad{\rm and } \nonumber\\ V_{ab} &=& {-}\frac{1}{E_0} \int d^3 x \, \rho \, {\boldsymbol \xi }_a \cdot ({\boldsymbol \xi }_b \cdot {\boldsymbol \nabla }) {\boldsymbol \nabla } V. \end{eqnarray} \tag{ 51 }$$

Defining an expansion in for σ_a which must begin at O(²) because the perturbing tidal field enters only at this order,

$$\begin{eqnarray} &&\sigma _a = \epsilon ^2 \sigma _a^{(2)} + O (\epsilon ^3), \end{eqnarray} \tag{ 52 }$$

we can again solve for σ_a for this static configuration by minimizing the variation in $\mathcal {V}$ with respect to the amplitude σ_a. At leading order, we simply have

$$\begin{eqnarray} &&\sigma _a^{(2)} = V_a. \end{eqnarray} \tag{ 53 }$$

It is worth noting here that since ${\boldsymbol \sigma }^{(1)} = 0$, we immediately have that

$$\begin{eqnarray} &&{\boldsymbol \chi }^{(1)} = - {\boldsymbol \zeta }^{(1)}. \end{eqnarray} \tag{ 54 }$$

To see this, we note that by definition the displacement of a point x in the unperturbed star by σ is the same as the composition of a displacement by χ and then the application of the coordinate transform ψ. Keeping only the leading-order terms in this gives

$$\begin{eqnarray} &&{\boldsymbol x} + {\boldsymbol \sigma }({\boldsymbol x}) = {\boldsymbol x} + O(\epsilon ^2) = {\boldsymbol x} + \epsilon {\boldsymbol \chi }^{(1)}({\boldsymbol x}) + \epsilon {\boldsymbol \zeta }^{(1)}({\boldsymbol x}) + O(\epsilon ^2), \nonumber\\ \end{eqnarray} \tag{ 55 }$$

from which Equation (54) follows.

As before, we now consider further perturbations to star B, which we denote ${\boldsymbol \eta }_S$. From the original unperturbed star we have a displacement

$$\begin{eqnarray} {\boldsymbol x} &\rightarrow & {\boldsymbol x} + {\boldsymbol \xi } = {\boldsymbol x} + {\boldsymbol \sigma } + {\boldsymbol \eta }_S \quad{\rm where} \end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray} {\boldsymbol \eta }_S &= & {\boldsymbol J}_\psi \cdot {\boldsymbol \eta }. \end{eqnarray} \tag{ 57 }$$

Here, η is the form of the corresponding oscillations of star A, which must be transformed into the coordinates X of the spherical star B if we are to consider further displacements beyond σ. This is the origin of the Jacobian transformation matrix ${\boldsymbol J}_\psi$ in Equation (57). In principle, Equation (57) should include terms of order ${\boldsymbol \eta }^2$; however, since star B is an equilibrium solution, the potential energy contains no terms of first order in ${\boldsymbol \eta }_S$ and the leading dependence is ${\boldsymbol \eta }_S^2$. Therefore to compute the potential energy to order ${\boldsymbol \eta }^2$, we need only obtain ${\boldsymbol \eta }_S$ to linear order in η.

Now we consider the potential $\mathcal {V}$ of star B when the additional perturbations ${\boldsymbol J}_\psi \cdot {\boldsymbol \eta }$ are present. We define the expansion of the Jacobian in orders of ,

$$\begin{eqnarray} &&{\boldsymbol J}_\psi = {\boldsymbol 1} + \epsilon {\boldsymbol J}^{(1)}_\psi + \epsilon ^2 {\boldsymbol J}^{(2)}_\psi + O (\epsilon ^3), \end{eqnarray} \tag{ 58 }$$

and its expansion in the mode basis ${\boldsymbol \xi }_a$ of the original star,

$$\begin{eqnarray} &&J_{ab} = \frac{\omega _a^2}{E_0} \int d^3 x \, \rho \, {{\boldsymbol \xi }^*_a} \cdot {\boldsymbol J}_\psi \cdot {\boldsymbol \xi }_b, \end{eqnarray} \tag{ 59 }$$

so that the basis coefficients are naturally transformed as

$$\begin{eqnarray} {\boldsymbol \eta }_S &=& \sum \eta _{S,a} {\boldsymbol \xi }_a = \sum J_{ab} \eta _b {\boldsymbol \xi }_a \nonumber\\ &=& \sum \left(\delta _{ab} + \epsilon J_{ab}^{(1)} + \epsilon ^{2} J_{ab}^{(2)} \right) \eta _b {\boldsymbol \xi }_a. \end{eqnarray} \tag{ 60 }$$

Taking all of this into account, we have for star B the mode expansion of the potential

$$\begin{eqnarray} \mathcal {V} ({\boldsymbol \sigma } + {\boldsymbol \eta }_S) &= & \mathcal {V}(0) + \int d^3 x \rho \epsilon ^2 V \nonumber\\ && + \frac{E_0}{2} \sum \left[ \eta _a^2 + \epsilon \left(J^{(1)}_{ab} + J^{(1)}_{ba}\right) \eta _a \eta _b \right.\nonumber\\ && + \epsilon ^2 \left(J^{(1)}_{ca} J^{(1)}_{cb} + J^{(2)}_{ab} + J^{(2)}_{ba}\right) \eta _a \eta _b \nonumber \\ &&\left.- 2\kappa _{abc} \eta _a \eta _b \sigma _c - \epsilon ^2 V_{ab} \eta _a \eta _b + \cdots \right] \nonumber \\ &=& \mathcal {V}(0) + \int d^3 x \rho \epsilon ^2 V \nonumber\\ && + \frac{E_0}{2} \sum \left[ \eta _a^2 + \epsilon \left(J^{(1)}_{ab} + J^{(1)}_{ba}\right) \eta _a \eta _b \right.\nonumber \\ &&+ \epsilon ^2 \left(J^{(1)}_{ca} J^{(1)}_{cb} +J^{(2)}_{ab} + J^{(2)}_{ba} \right.\nonumber\\ &&\left.\left.- 2\kappa _{abc} V_c -V_{ab}\right) \eta _a \eta _b \right], \end{eqnarray} \tag{ 61 }$$

where again in the first equality we have already eliminated the terms which vanish when the solution for σ_a is inserted and terms higher order in or η than we are considering. In the second equality we have substituted the solution (53) for σ_a and collected terms according to their order in .

We have computed an expression for $\mathcal {V}$ in the mode basis using two different methods, and we can now equate these expressions. For star A, $\mathcal {V}$ is given by the Lagrangian of Equation (17), through $\mathcal {V} = \sum \dot{\eta }_a^2/\omega _a^2- \mathcal {L}$. Matching this to Equation (61) order by order in and η, we get the following conditions on the three- and four-mode coupling terms.

At order , matching terms gives

$$\begin{eqnarray} &&U_{ab}+ \sum 2 \kappa _{abc} \chi _c^{(1)} = -\left(J^{(1)}_{ab} + J^{(1)}_{ba}\right), \end{eqnarray} \tag{ 62 }$$

which results in an expression for $\kappa _{abc} \chi _c^{(1)}$ in terms of the Jacobian transform and the tidal potential. Meanwhile, matching the ² terms gives

$$\begin{eqnarray} && \sum \left(2 \kappa _{abc} \chi _c^{(2)} + 3 \kappa _{abcd} \chi _c^{(1)} \chi _d^{(1)} \right)\nonumber\\ && \quad = {-} \sum \left(J^{(1)}_{ca} J^{(1)}_{cb} + J^{(2)}_{ab} + J^{(2)}_{ba} - 2\kappa _{abc} V_c -V_{ab} \right). \quad \end{eqnarray} \tag{ 63 }$$

In both of Equations (62) and (63), the indices (a, b) refer to particular modes and are not summed over.

Before using the matching conditions in Equations (62) and (63), we derive the explicit equations for the Jacobian that transforms the perturbations η on star A to the perturbations ${\boldsymbol \eta }_S$. This is achieved by considering the difference between two series of active transformations as depicted in Figure 2, using a fixed background coordinate system. We begin by considering a fluid element at a point x in the unperturbed star. We then consider the application of the tidal perturbation, which takes ${\boldsymbol x} \rightarrow {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})$, followed by the coordinate transform ψ acting on this element, which acts on the element at its displaced position ${\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})$. This takes us across the top row of Figure 2, and acts to transform the initial point as

$$\begin{eqnarray} &&{\boldsymbol x} \rightarrow {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x}) + \left. {\boldsymbol \zeta } \right|_{ {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})} + \frac{1}{2} \left. ({\boldsymbol \zeta }\cdot {\boldsymbol \nabla }){\boldsymbol \zeta }\right|_{ {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})}. \end{eqnarray} \tag{ 64 }$$

Next, consider the series of transforms which first takes x to the corresponding point in star A, ${\boldsymbol x} \rightarrow {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x}) + {\boldsymbol \eta }({\boldsymbol x})$. Since we wish to express all perturbations in Lagrangian coordinates, the η is written as a function of the position of the original, unperturbed element. We follow this pair of displacements by the coordinate transform ψ, which acts on the actual position of the element. This gives us the position of the element in star B,

$$\begin{eqnarray} &&{\boldsymbol x} \rightarrow {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x}) + {\boldsymbol \eta }({\boldsymbol x}) + \left. {\boldsymbol \zeta } \right|_{ {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})+ {\boldsymbol \eta }({\boldsymbol x})} + \frac{1}{2} \left. ({\boldsymbol \zeta }\cdot {\boldsymbol \nabla }){\boldsymbol \zeta }\right|_{ {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})+ {\boldsymbol \eta }({\boldsymbol x})}. \nonumber\\ \end{eqnarray} \tag{ 65 }$$

The difference between Equations (65) and (64) is simply ${\boldsymbol \eta }_S({\boldsymbol x}) = {\boldsymbol J}_\psi {\boldsymbol \eta }({\boldsymbol x})$. Taking this difference, and expanding our expressions in terms of small η, we have

$$\begin{eqnarray} {\boldsymbol J}_\psi {\boldsymbol \eta } &=& {\boldsymbol \eta }({\boldsymbol x}) + \left. {\boldsymbol \eta }({\boldsymbol x}) \cdot ({\boldsymbol \nabla } {\boldsymbol \zeta })\right|_{ {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})} \nonumber\\ &&+ \frac{1}{2} \left. {\boldsymbol \eta }({\boldsymbol x}) \cdot [ {\boldsymbol \nabla } ({\boldsymbol \zeta }\cdot {\boldsymbol \nabla }){\boldsymbol \zeta }\right]_{ {\boldsymbol x} + {\boldsymbol \chi }({\boldsymbol x})} + O(\eta ^2). \end{eqnarray} \tag{ 66 }$$

Further expanding out χ and ζ in terms of , this expression becomes

$$\begin{eqnarray} {\boldsymbol J}_\psi {\boldsymbol \eta } &=& {\boldsymbol \eta }+ \epsilon ({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta }^{(1)} \nonumber\\ &&+ \epsilon ^2 \left(({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta }^{(2)}+ {\boldsymbol \eta } \cdot [({\boldsymbol \chi }^{(1)} \cdot {\boldsymbol \nabla }){\boldsymbol \nabla } {\boldsymbol \zeta }^{(1)})] \right.\nonumber\\ &&+\left. \frac{1}{2} ({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) ({\boldsymbol \zeta }^{(1)}\cdot {\boldsymbol \nabla }){\boldsymbol \zeta }^{(1)} \right) + O (\eta ^3). \end{eqnarray} \tag{ 67 }$$

Here, all terms are evaluated at the base point x. This allows us to simply read off the Jacobian of Equation (58). It is useful to recall that ${\boldsymbol \chi }^{(1)} = - {\boldsymbol \zeta }^{(1)}$, which allows us to simplify the Jacobian somewhat:

$$\begin{eqnarray} {\boldsymbol J}^{(1)}_\psi {\boldsymbol \eta } &=& ({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta }^{(1)} \quad{\rm and } \; {\boldsymbol J}^{(2)}_\psi {\boldsymbol \eta } = ({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta }^{(2)} + [({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) {\boldsymbol \zeta }^{(1)}] \nonumber\\ && \cdot {\boldsymbol \nabla } {\boldsymbol \zeta }^{(1)} - \frac{1}{2} ({\boldsymbol \eta } \cdot {\boldsymbol \nabla }) ({\boldsymbol \zeta }^{(1)} \cdot {\boldsymbol \nabla }){\boldsymbol \zeta }^{(1)}. \end{eqnarray} \tag{ 68 }$$

To simplify Equation (67) we have temporarily resorted to a Cartesian basis in order to commute the covariant derivatives.

Now that we have an explicit expression for the Jacobian, we can express it in the mode basis using the expansion for η. The result, for the first- and second-order terms, is

$$\begin{eqnarray} J_{ab}^{(1)}&=&\frac{\omega _a^2}{E_0}\int d^3x \rho \, {\boldsymbol \xi }_a\cdot ({\boldsymbol \xi }_b\cdot {\boldsymbol \nabla }){\boldsymbol \zeta }^{(1)} \quad{\rm and} \end{eqnarray} \tag{ 69 }$$

$$\begin{eqnarray} J_{ab}^{(2)}&=&\frac{\omega _a^2}{E_0}\int d^3x \rho \, {\boldsymbol \xi }_a\cdot ({\boldsymbol \xi }_b\cdot {\boldsymbol \nabla })\left[{\boldsymbol \zeta }^{(2)}-\frac{1}{2}({\boldsymbol \zeta }^{(1)}\cdot { \boldsymbol \nabla }){\boldsymbol \zeta }^{(1)}\right] \nonumber\\ &&+ {\boldsymbol \xi }_a\cdot ([ ({\boldsymbol \xi }_b\cdot {\boldsymbol \nabla }){\boldsymbol \zeta }^{(1)}] \cdot {\boldsymbol \nabla }{\boldsymbol \zeta }^{(1)}). \end{eqnarray} \tag{ 70 }$$

For our initial check of the mode stability in Section 2.4, we note here that for a particular high-order mode, $\xi _a \sim \omega _a^{-1}$, and so we have the useful fact that for a particular pair of modes (a, b),

$$\begin{eqnarray} &&J^{(i)}_{ab} \sim \frac{\omega _a}{\omega _b} \end{eqnarray} \tag{ 71 }$$

so long as the angular integrations satisfy selection rules and the contraction of indices in Equations (69) and (70) do not lead to a much smaller value. This fact is made more explicit in Section 4.

With our matching results from Section 3.4, we can return to the expressions for the perturbed eigenvalues from Section 2.4. We focus on the case of interest for the instability proposed by WAB, that of a high-order (p, g) mode pair with comparable wave numbers and widely spaced frequencies, ω_p ≫ ω_g. Working in the rescaled mode amplitudes $\eta _a^{\prime } = \eta _a/\omega _a$, consider the eigenvalues of the matrix $\mathcal {M}$ in the Lagrangian (26). Substituting the matching conditions (62) and (63) into $\mathcal {M}$ as given by Equations (27)–(29), we have now that

$$\begin{eqnarray} \mathcal {M}_{pp} &= & \omega _p^2 + 2 \epsilon \omega _p^2 J^{(1)}_{pp} + \epsilon ^2 \omega _p^2 \left[ \sum \left(J^{(1)}_{ap} J^{(1)}_{ap} - 2\kappa _{ppa} V_a \right)\right. \nonumber\\ &&\left. +\, 2 J^{(2)}_{pp} -V_{pp} \right], \end{eqnarray} \tag{ 72 }$$

$$\begin{eqnarray} \mathcal {M}_{pg} &= & \mathcal {M}_{gp} = \epsilon \, \omega _p \omega _g \left(J^{(1)}_{pg} + J^{(1)}_{gp} \right), \quad {\rm and} \end{eqnarray} \tag{ 73 }$$

$$\begin{eqnarray} \mathcal {M}_{gg} &= & \omega _g^2 + 2 \epsilon \omega _g^2 J^{(1)}_{gg} \nonumber\\ &&+ \epsilon ^2 \omega _g^2 \left[ \sum \left(J^{(1)}_{ag} J^{(1)}_{ag} - 2\kappa _{gga} V_a \right) + 2 J^{(2)}_{gg} -V_{gg} \right]. \nonumber\\ \end{eqnarray} \tag{ 74 }$$

The resulting perturbed eigenvalues are similarly re-expressed in terms of the Jacobian and the potential V. Let us focus on the smaller frequency, which is perturbed to a lower (and possibly negative value):

$$\begin{eqnarray} \frac{\omega _-^2}{\omega _g^2} &= &1 + 2 \epsilon J_{gg}^{(1)}\nonumber\\ && + \epsilon ^2 \left[\left(J^{(1)}_{pg}\right)^2 +\left(J^{(1)}_{gg}\right)^2 + 2 J^{(2)}_{gg} -V_{gg} -\sum 2\kappa _{gga} V_a \right] \nonumber\\ && - \epsilon ^2 \frac{\omega _p^2}{\omega _p^2- \omega _g^2} \left(J^{(1)}_{pg} + J^{(1)}_{gp} \right)^2. \end{eqnarray} \tag{ 75 }$$

Keeping in mind the origin of each of the terms, and our estimate that $J^{(i)}_{ab} \sim \omega _a/\omega _b$, we can see how Equation (75) encodes a potential instability, and the particular manner in which it is actually canceled. The final O(²) term in Equation (75) is just the three-mode term $\kappa _{pga} \chi _a^{(1)}$ in Equation (31), expressed in the language of Jacobians. When the frequency of the p-mode is much larger than that of the g-mode, the size of this term is ${\sim} (J^{(1)}_{pg})^2 \sim (\omega _p^2/\omega _g^2) \gg 1$. In principle, this can overcome the ² suppression when the tidal field is relatively strong (but with ≪ 1), and overwhelm the order-unity contribution from the restoring force. However, when we express the perturbation to the diagonal terms (entering from the four-mode and $\kappa _{abc}\chi ^{(2)}_c$ terms in Equation (31)) in terms of Jacobians, we observe that it contains an identical contribution with the opposite sign.

Before we explicitly write down the remaining terms, it is worthwhile to pause and consider the impact of our coordinate transformation. The three- and four-mode terms in Equation (31) both affect the eigenvalues equally. Independently calculating their individual contributions would have involved careful book-keeping in order to accurately track the cancellation of large terms. Instead, our approach has confirmed the intuition developed from the toy model of Section 1.1, and illuminated the fact that large coupling terms can arise from rotations of modes into each other.

We now write down the terms that remain after the large cancellations. Expanding the prefactor of the final term in Equation (75),

$$\begin{eqnarray} &&\frac{\omega _p^2}{\omega _p^2- \omega _g^2} = \left(1 + \frac{\omega _g^2}{\omega _p^2} \right) + {\cal O}\left(\frac{\omega _g^4}{\omega _p^4}\right), \end{eqnarray} \tag{ 76 }$$

we can simplify Equation (75) to

$$\begin{eqnarray} \frac{\omega _-^2}{\omega _g^2} &= &1 + 2 \epsilon J^{(1)}_{gg} + \epsilon ^2 \left[\frac{\omega _g^2}{\omega _p^2}\left(J^{(1)}_{pg}\right)^2+\left(J^{(1)}_{gg}\right)^2 + 2 J^{(2)}_{gg} -V_{gg} \right.\nonumber\\ &&\left.-2J^{(1)}_{pg}J^{(1)}_{gp} -\sum 2\kappa _{gga} V_a \right]. \end{eqnarray} \tag{ 77 }$$

Our estimate for the size of the Jacobian terms shows us that all the terms are ${\cal O}(1)$ in terms of frequency, except perhaps the V_gg and κ_ggaV_a terms. We need to check that these terms are not large for the case of the high-order (p, g)-mode coupling. We show this formally and estimate the size of these terms in Section 4.

In both this section and in Section 2.4, we have considered the coupling of a single pair of modes (p, g), but in general $\mathcal {M}$ contains entries for all modes and their mode–mode coupling terms. Equations (72)–(77) extend naturally to this case. For example, if we consider a single g-mode, there is an off-diagonal matrix entry $\mathcal {M}_{a g}$ for every other mode η_a, and this entry couples the g-mode to η_a. In addition, the sum $\sum J^{(1)}_{ag} J^{(1)}_{ag}$ in $\mathcal {M}_{gg}$ contains a contribution from each of these other modes. When computing the perturbation to the g-mode frequency, the off-diagonal terms all enter as a sum of squared terms analogous to the last term in Equation (75). In those cases where a mode strongly couples to the g-mode, the same leading-order cancellation of large terms occurs mode-by-mode. This means that even in the case where many p-modes strongly couple to a single g-mode, cancellation between the large parts of the three-mode couplings (squared) and four-mode couplings prevents the g-mode from developing an instability.

It is remarkable that up until now, our investigation of the generation of the daughter modes from the tidal perturbation has been formal and general in nature, with no reference made to a particular stellar model. In our estimates we have only needed to note that ${\boldsymbol \xi }_a \sim \omega _a^{-1}$ given our normalization of the mode functions.

We begin by considering the amplitude $U_a = \chi ^{(1)}_a = - \zeta ^{(1)}_a$. It is conveniently obtained by considering the leading-order part of ζ,

$$\begin{eqnarray} &&{\boldsymbol \zeta }^{(1)} = - {\boldsymbol \chi }^{(1)} = - \omega _0^2\frac{r^2}{\mathfrak {g}} \left[ P_2(\cos \theta) \hat{\boldsymbol r} + \frac{4-n}{6} \partial _\theta P_2(\cos \theta) \hat{\boldsymbol \theta } \right].\nonumber\\ \end{eqnarray} \tag{ 78 }$$

When we recall the definition that $\epsilon = Gm/(a^3\omega _0^2)$ this recovers the linear tidal response,

$$\begin{eqnarray} &&\epsilon {\boldsymbol \chi }^{(1)} = \frac{Gm}{a^3} \frac{r^2}{\mathfrak {g}} \left[ P_2(\cos \theta) \hat{\boldsymbol r} + \frac{4-n}{6} \partial _\theta P_2(\cos \theta) \hat{\boldsymbol \theta } \right], \end{eqnarray} \tag{ 79 }$$

which can be compared to Equations (A12) and (A13) of Weinberg et al. (2012). Recall that our eigenmodes ${\boldsymbol \xi }_a$ have the form in Equation (10), and that here we have specialized to the m = 0 case. We see that ${\boldsymbol \chi }^{(1)}$ has the correct form for an l = 2 displacement.

Before continuing, it is useful to note that we have frequent need of integrals of the form

$$\begin{eqnarray} I_{abc} &= & \int d^3 x \, \rho \, {\boldsymbol \xi }_a \cdot ({\boldsymbol \xi }_b \cdot {\boldsymbol \nabla }) {\boldsymbol \xi _c}. \end{eqnarray} \tag{ 80 }$$

With our convention that ${\boldsymbol \xi }_a$ take the form in Equation (10), the integral (80) resolves to

$$\begin{eqnarray} I_{abc} &=& \int dr r^2 \rho (r) \left[ T_{abc} a_r b_r \frac{d c_r}{d r} + F_{a,bc} \frac{a_r b_h}{r} (c_r - c_h)\right. \nonumber\\ &&\left.+\, F_{b,ac} a_h b_r \frac{d c_h}{d r} + \frac{a_h b_h}{r}(G_{c,ab} c_h + F_{c,ab} c_r) \right]. \end{eqnarray} \tag{ 81 }$$

The angular integrals T_abc, F_{a, bc}, and G_{a, bc} are those originally defined in Wu & Goldreich (2001) and match those listed in Equations (A20)–(A22) of Weinberg et al. (2012). We give them in Appendix D, and they are determined by the angular indices (l_a, l_b, l_c) of the mode vectors in the integral. These integrals vanish unless the angular momentum indices obey the triangle inequality |l_b − l_c| ⩽ l_a ⩽ l_b + l_c, and in addition the sum l_a + l_b + l_c must be even.

Now we are ready to consider the three-mode coupling term for the case of a pair of (p, g) modes excited by the tide. From Equation (62), we have

$$\begin{eqnarray} &&\sum \kappa _{pgc} \chi _c^{(1)} = - \frac{1}{2} U_{pg} - \frac{1}{2} \left(J^{(1)}_{pg} + J^{(1)}_{gp}\right). \end{eqnarray} \tag{ 82 }$$

The coupling κ_pgc can be large in the case ω_p ≫ ω_g, and when the spatial resonance condition k_p ≃ k_g is met (Weinberg et al. 2013). This means that the (p, g) modes must be high-order modes with large wave numbers, and as such we use the WKB approximation (Unno et al. 1989; Weinberg et al. 2012) to get analytic forms for the eigenfunctions ${\boldsymbol \xi }_p$ and ${\boldsymbol \xi }_g$. In this case, the radial functions (p_r, p_h) and (g_r, g_h) are

$$\begin{eqnarray} (p_r, p_h)&\simeq & \frac{A_p}{\omega _p}\left(\cos {k_p r},\frac{c_s \sin {k_p r}}{\omega _p r}\right)\quad{\rm and} \end{eqnarray} \tag{ 83 }$$

$$\begin{eqnarray} (g_r, g_h)&\simeq & \frac{A_g}{\omega _g}\left(\frac{\omega _g \sin {k_g r}}{N},\frac{\cos {k_g r}}{\Lambda _g}\right), \end{eqnarray} \tag{ 84 }$$

where c_s is the adiabatic sound speed, N is the Brünt–Väisälä frequency,

$$\begin{eqnarray} N^2&=&{-}\left(\frac{1}{\rho }\frac{d\rho }{dr}-\frac{1}{\Gamma _1 P}\frac{dP}{dr}\right)\mathfrak {g}, \end{eqnarray} \tag{ 85 }$$

Γ₁ is the adiabatic index, $\Lambda _a^2=l_a(l_a+1)$, and we set

$$\begin{eqnarray} (k_p,k_g) &\simeq& \left(\frac{\omega _p}{c_s},\frac{{\Lambda _g} N}{r\omega _g}\right), \quad A_{p,g}=\sqrt{\frac{E_0\alpha _{p,g}}{\rho r^2}}, \nonumber \\ \alpha _p &=& \left(c_s\int c_s^{-1} dr\right)^{-1}, \quad {\rm and} \nonumber\\ \alpha _g &=& \frac{N}{r}\left(\int N d\ln {r}\right)^{-1}. \end{eqnarray} \tag{ 86 }$$

In the WKB approximation, k_p ≃ k_g implies that ω_pω_g ≃ Λ_gNc_s/r. The key consideration for estimating the largest terms is that for the p-mode, p_r is the dominant component of the mode (the acoustic modes are mostly radial) and for the g-mode, g_h is the dominant component (the gravity wave modes are mostly horizontal).

When we have need for a specific stellar model, we use rough approximations to the neutron star model (Steiner & Watts 2009) used by WAB, generated with the Skyrme–Lyon equation of state (SLy4; Chabanat et al. 1998). In this model the density is approximately constant throughout the core of the neutron star, where the coupling κ_pgc is large. When this is true the gravitational acceleration grows linearly with the radius, $\mathfrak {g} \simeq 4 \pi G \rho r/3$. In addition, c_s is roughly constant throughout the core of the star and so c_s ≃ ω₀R_*. We also use the expression for N from Reisenegger & Goldreich (1992), derived to leading order in small electron fraction Y_e,

$$\begin{eqnarray} &&N \simeq \frac{\mathfrak {g}}{c_s}\sqrt{\frac{Y_e}{2}}. \end{eqnarray} \tag{ 87 }$$

This leads us to see that Nc_s/r is constant, and moreover that N ∼ ω₀r/R_*, since $\omega _0^2 \simeq G \rho$ for a nearly constant density. Taken together with the condition that the wave numbers for the (p, g) modes are nearly equal, we have that $\omega _p \omega _g \sim \omega _0^2$. In this model, Equation (78) shows that $\chi ^{(1)}_{r}\sim \chi ^{(1)}_{h} \sim r$. These approximations are expected to hold for a variety of neutron star models (see WAB; Section 3.3).

We begin our computation of the three-mode coupling constant by evaluating U_pg,

$$\begin{eqnarray} U_{pg} &=& {-} \frac{\omega _0^2}{E_0} \int d^3 x \, \rho \, {\boldsymbol \xi }_p \cdot ({\boldsymbol \xi }_g \cdot {\boldsymbol \nabla }){\boldsymbol \nabla } [r^2 P_2(\cos \theta) ] \nonumber \\ &=& {-} \frac{\omega _0^2}{E_0} \sqrt{\frac{4\pi }{5}} \int dr \, r^2 \rho [ 2 p_r g_r T_{pg2} + p_r g_h F_{p,g2} \nonumber\\ &&+ p_h g_r F_{g,p2} + p_h g_h (2 F_{2,pg}+ G_{2,pg})] \nonumber \\ &\simeq & {-} \frac{\omega _0^2}{E_0} \sqrt{\frac{4\pi }{5}} \int dr r^2 \rho \, p_r g_h F_{p,g2} \sim \frac{\omega _0^2}{\omega _p \omega _g} \sim {\cal O}(1), \nonumber\\ \end{eqnarray} \tag{ 88 }$$

where in the third line terms have been dropped since they are higher order in ω_g/ω₀ or ω₀/ω_p. We see that U_pg contains no large terms.

Similarly, we can consider the leading-order Jacobian terms, using Equations (69) and (81) to write them as,

$$\begin{eqnarray} J_{ab}^{(1)} &= & \frac{\omega _a^2}{E_0} I_{a b \zeta } = -\frac{\omega _a^2}{E_0} I_{a b \chi ^{(1)}}, \end{eqnarray} \tag{ 89 }$$

where the subscripts ζ and χ⁽¹⁾ indicate the substitution of the corresponding displacement for the mode function ${\boldsymbol \xi }_c$ in the integral I_abc. The radial functions χ_r and χ_h are defined as in Equation (10), with l = 2 and m = 0, so that the largest contributions to the Jacobian terms are

$$\begin{eqnarray} J_{pg}^{(1)} &= & {-} \frac{\omega _p^2}{E_0} \int dr \, r^2 \rho \, F_{p,g2} \frac{p_r g_h}{r} \left(\chi _r^{(1)} - \chi _h^{(1)}\right) \sim \frac{\omega _p}{\omega _g} \quad {\rm and} \nonumber\\ \end{eqnarray} \tag{ 90 }$$

$$\begin{eqnarray} J_{gp}^{(1)} &= & {-} \frac{\omega _g^2}{E_0} \int dr \, r^2 \rho \, F_{g,p2} \, g_h p_r \frac{d\chi _h^{(1)}}{dr} \sim \frac{\omega _g}{\omega _p}. \end{eqnarray} \tag{ 91 }$$

At leading order, then,

$$\begin{eqnarray} &&\sum \kappa _{pgc} \chi _c^{(1)} \sim \frac{\omega _p^2}{2E_0} \int dr \, r^2 \rho \, F_{p,g2} \frac{p_r g_h}{r} \left(\chi _r^{(1)} - \chi _h^{(1)} \right) \sim \frac{\omega _p}{\omega _g}. \nonumber\\ \end{eqnarray} \tag{ 92 }$$

This matches the leading-order terms in WAB, and arises solely from our consideration of the volume-preserving transform. As we have shown, a correction involving the four-mode coupling term cancels its influence on the eigenfrequencies. We now show that none of the remaining terms are large enough to produce a potential instability.

We have calculated the correction to the frequency of the almost-neutrally stable g-modes in the presence of the tidal deformation, in Equation (77). We can write the relative corrections as

$$\begin{eqnarray} && \frac{\omega _-^2 - \omega _g^2}{\omega _g^2} = 2 \epsilon J^{(1)}_{gg} \nonumber\\ &&\quad + \epsilon ^2 \left[ \left\lbrace \frac{\omega _g^2}{\omega _p^2}\left(J^{(1)}_{pg}\right)^2+\left(J^{(1)}_{gg}\right)^2 -2J^{(1)}_{pg}J^{(1)}_{gp} + 2 J^{(2)}_{gg}\right\rbrace \right.\nonumber\\ &&\left.\quad - \left\lbrace V_{gg} + \sum 2\kappa _{gga} V_a \right\rbrace \right]. \end{eqnarray} \tag{ 93 }$$

The corrections can be divided into two classes - a set of Jacobian terms, and a set of potential terms.

We have previously encountered Jacobian terms while estimating the three-mode coupling. For high-order modes a and b, the leading order dependence of the first-order Jacobian $J^{(1)}_{ab}$ on the frequencies is $J^{(1)}_{a b} \sim {\cal O}(\omega _a/\omega _b)$. From this observation, we see that all the terms involving a first-order Jacobian in Equation (93) are ${\cal O}(1)$ in the large frequency ratio ω_p/ω_g.

The second-order Jacobian term $J^{(2)}_{ab}$ is given by Equation (70). We observe that it has the same frequency dependence as the first-order Jacobian, with a prefactor of $\omega _a^2/E_0$ and the two eigenfunctions ${\boldsymbol \xi }_a$ and ${\boldsymbol \xi }_b$ present in the integrand. As we did for the first-order Jacobian, we recall that for high-order modes the WKB eigenfunctions have a size $\xi _a \sim \sqrt{E_0}/ \omega _a$, as seen in Equations (83) and (84). Hence, the second-order Jacobian term is ${\cal O}(1)$ in the frequency ratio ω_p/ω_g.

Having established that all the Jacobian terms in Equation (77) are of the form or ² times terms of ${\cal O}(1)$ in the frequency ratio, we turn to the potential terms in Equation (93). Using Equation (53), we can write

$$\begin{eqnarray} -\epsilon ^2 \left(V_{gg} + \sum 2\kappa _{gga}V_a \right) &=& {-}\epsilon ^2 \left(V_{gg} + \sum 2\kappa _{gga}\sigma ^{(2)}_a \right) \nonumber\\ &=& {-}\epsilon ^2 (V_{gg} + 2\kappa _{gg\sigma }), \end{eqnarray} \tag{ 94 }$$

where in the last equality we have used the notation introduced in Equation (89), and further dropped the superscript on σ with the understanding that we are using the leading ² term in the expansion of the displacement. The key point is that ²κ_ggσ can be computed by substituting the displacement σ in place of the eigenfunction ${\boldsymbol \xi }_c$ in the definition (15) of the three-mode coupling, as κ_abc is linear in its arguments.

The nonlinear tidal term V_gg is given by Equation (51)

$$\begin{eqnarray} V_{gg} &= & {-} \frac{1}{E_0} \int d^3 x \, \rho \, {\boldsymbol \xi }_g \cdot ({\boldsymbol \xi }_g \cdot {\boldsymbol \nabla }) {\boldsymbol \nabla } V \nonumber \\ &= & {-} \sqrt{4\pi } \frac{1}{E_0} \int dr\, r^2 \rho (r) \nonumber\\ &&\times \left(g_r g_r \frac{d^2 V}{d r^2} T_{\sigma gg} + \frac{g_h g_h}{r} \frac{dV}{dr} F_{\sigma, gg} \right) \nonumber \\ &= & {-} \frac{1}{E_0} \int dr\, r^2 \rho (r) \left(g_r g_r \frac{d^2 V}{d r^2} + \Lambda _g^2 \frac{g_h g_h}{r} \frac{dV}{dr} \right). \nonumber\\ \end{eqnarray} \tag{ 95 }$$

In this expression, we have used expressions for the angular integrals from Wu & Goldreich (2001), together with the fact that the radial displacement σ has the angular quantum numbers l = 0, m = 0. The integrals are elaborated upon in Appendix D. Note that in this expression for V_gg, we have a potentially large part ∼g_hg_h due to the large size of the horizontal displacement g_h of the g-mode. We see that this term cancels exactly with a part of the three-mode term.⁴

The three-mode term κ_ggσ has a covariant form (e.g., Schenk et al. 2002, Equation (4.20)). It can be evaluated in terms of the radial and angular parts of the mode eigenfunctions (Wu & Goldreich 2001; Weinberg et al. 2012). Before continuing, we should note a complication which arises in the case we are dealing with, namely that of a radial displacement coupled to nonradial displacements.

The process of expanding and simplifying the expression for the three-mode coupling involves using the equations of motion for the constituent displacements. The equation for the divergence of a displacement of the form of Equation (10) is

$$\begin{eqnarray} ({\boldsymbol \nabla } \cdot {\boldsymbol \xi }_a)_r &= & \frac{d a_r}{d r} + \frac{2}{r}a_r - \frac{\Lambda _a^2}{r} a_h, \end{eqnarray} \tag{ 96 }$$

where we have written the divergence as $({\boldsymbol \nabla } \cdot {\boldsymbol \xi }_a)_r$ because we have omitted its spherical harmonic angular dependence, which is absorbed in the angular integrals. The equations of motion for a mode with angular index l > 0 are

$$\begin{eqnarray} \Gamma _1 P({\boldsymbol \nabla } \cdot {\boldsymbol \xi }_a)_r &= & {-} (\delta P)_r = \rho \mathfrak {g} a_r - \omega _a^2 \rho r a_h + \rho \Phi ^\prime \quad{\rm and} \nonumber\\ \end{eqnarray} \tag{ 97 }$$

$$\begin{eqnarray} && \frac{d}{dr} [\Gamma _1 P({\boldsymbol \nabla } \cdot {\boldsymbol \xi }_a)_r] = \frac{\Lambda _a^2}{r}\rho \mathfrak {g} a_h \nonumber\\ &&\quad - \left(\omega _a^2 + \frac{2\mathfrak {g}}{r} - \frac{d \mathfrak {g}}{d r} \right) \rho a_r + \rho \frac{d \Phi ^\prime }{d r}. \end{eqnarray} \tag{ 98 }$$

The Eulerian perturbation of the potential is denoted by Φ'. Within the Cowling approximation, it has a contribution only from the external driving potential V. Spherical symmetry demands that it only has a contribution from the part of the driving potential with the same spherical harmonic dependence as the displacement. For the g-mode, the external driving potential V_g vanishes, because the mode's eigenfunction is nonradial, and the potential is spherically symmetric.

We cannot use Equations (96)–(98) for the radial displacement σ just by setting the angular displacement σ_h to zero. Deriving Equation (97) involves using the angular parts of the equations of motion, which do not exist for a radial displacement. We could still have used Equation (97) for the radial displacement had we been operating within the hydrostatic approximation, with the Eulerian pressure perturbation P' = ρgσ_r = −ρΦ' and the Lagrangian pressure perturbation δP = 0. We cannot use the approximation, as it is not consistent with our construction of the radial displacement as that part of the tidal displacement χ which changes the volume of the elements. The root of the difference is that radial and angular modes have different analytic structures. For instance, the Lagrangian and Eulerian pressure perturbations δP and P' do not necessarily vanish at the center of the star for a radial displacement. This fact is noted in Cox (1980, Section 17.6). Another point of view is that for just a single degree of freedom, Equations (96)–(98) are overdetermined, and would need to satisfy a consistency condition.

We can still use Equation (96) (the definition of divergence) and Equation (98) (the radial equation of motion) for the radial displacement σ, with the Eulerian perturbation to the potential Φ' given by the external potential V on star B. Equation (98) still holds even though σ does not represent a normal mode of the star because it represents a force balance: star B is in hydrostatic equilibrium in the modified potential and hence we may use Equation (98) with ω_σ = 0.

Keeping this in mind, the starting point for simplifying the three-mode term κ_ggσ is Equations (A27)–(A30) of Weinberg et al. (2012). Our paths diverge from the point where the equations of motion are used. The end goal of our simplification is to get an expression which consists of terms which have dominant contributions of the form g_rg_r and g_hg_h, since given the WKB forms of the mode function, these are the forms which pick up a growing contribution as we integrate through the star. In order to do so, we repeatedly use the equations of motion for the nonradial mode ${\boldsymbol g}$ and the radial displacement σ, and integrate by parts wherever needed.

The final expression for the combination of coupling terms needed is

$$\begin{eqnarray} && \epsilon ^2 (V_{g g} + 2 \kappa _{\sigma g g}) = \frac{1}{E_0} \int dr \left[ r^2 P \left\lbrace \Gamma _1 (\Gamma _1 + 1) + \left(\frac{\partial \Gamma _1}{\partial \ln \rho } \right)_s \right\rbrace \right.\nonumber \\ &&\quad \times ({\boldsymbol \nabla } \cdot {\boldsymbol \sigma }) ({\boldsymbol \nabla } \cdot {\boldsymbol g})_r ({\boldsymbol \nabla } \cdot {\boldsymbol g})_r - 4 r \sigma _r \Gamma _1 P ({\boldsymbol \nabla } \cdot {\boldsymbol g})_r ({\boldsymbol \nabla } \cdot {\boldsymbol g})_r \nonumber\\ &&\quad - 2 r^2 \left(\Lambda _g^2 \omega _g^2 \rho r g_h g_h + 2 g_r \Gamma _1 P ({\boldsymbol \nabla } \cdot {\boldsymbol g})_r \right) \frac{d}{d r} \left(\frac{\sigma _r}{r} \right) \nonumber\\ &&\quad - \rho \mathfrak {g} r^3 g_r g_r \frac{d^2}{d r^2} \left(\frac{\sigma _r}{r} \right) \nonumber \\ &&\quad + \left(- \rho \mathfrak {g}^2 r \frac{d}{d r} \left(\frac{\sigma _r}{ \mathfrak {g} } \right) + \rho r \epsilon ^2 \frac{d V}{d r} \right) \nonumber\\ &&\left.\quad \times \left(2 r g_r ({\boldsymbol \nabla } \cdot {\boldsymbol g})_r + g_r g_r \frac{d \ln {\rho }}{d \ln {r}} \right) \right]. \end{eqnarray} \tag{ 99 }$$

The process of simplifying the terms from their canonical forms, and the cancellation of the large term in the inhomogeneous driving V_gg with the three-mode coupling are demonstrated in more detail in Appendix D.

In estimating the size of the remaining terms in the above expression, we find it useful to approximate the divergence of the g-mode by using Equation (97) in the following way. First we note that $\omega _g^2 g_h$ is small, and that

$$\begin{eqnarray} &&({\boldsymbol \nabla } \cdot {\boldsymbol g})_r \simeq \frac{\rho }{\Gamma _1 P} \mathfrak {g} \, g_r = \frac{\mathfrak {g}}{c_s^2} g_r. \end{eqnarray} \tag{ 100 }$$

The radial eigenfunction for the g-mode is given by Equation (84) within the WKB approximation, which involves the Brünt–Väisälä frequency N. Recall that for typical equations of state N and the acceleration due to gravity $\mathfrak {g}$ grow nearly linearly with radius till well outside the core, N ∼ ω₀r/R_* and $\mathfrak {g} \sim r \omega _0^2$, and the sound speed c_s is nearly constant, with c_s ∼ ω₀R_*. The radial displacement σ is regular near the center; in fact σ_r/r is an analytic function everywhere including around the center r → 0. From these observations, we can check that none of the potential terms given by Equation (99) pick up large contributions as we integrate through the star, and that their contribution is of the order ${\sim} {\cal O}(1)$ in the large frequency ratio ω_p/ω_g.

To sum up, we have shown that all the corrections to the frequency of the almost-neutrally stable g-mode due to inhomogeneous driving and the lowest nonlinear interactions are small. Specifically, they are the form or ² times terms of ${\cal O}(1)$ in the frequency ratio ω_p/ω_g. Since is a small parameter [$\epsilon = \Omega ^2/\omega _0^2 = (m/M)(R_*/a)^3$, where Ω is the orbital frequency of the binary] during the early part of the in-spiral phase, the interaction of internal modes with the equilibrium tide does not cause them to go unstable in this region of parameter space.

In this case, we further suppose that $\omega ^{\prime }_e$ and $\omega ^{\prime }_d$ are the alternative roots associated with the same diagonal entries $A^{(0)}_{ee}$ and $A^{(0)}_{dd}$. To do this, we take ω − ω_e, ω − ω_d, and δA as perturbations and expand to second order. The determinant is given by the usual formula,

$$\begin{eqnarray} &&\det {\bf A}(\omega) = \sum _\pi (-1)^\pi A_{1\pi (1)}(\omega) A_{2\pi (2)}(\omega) \ldots A_{N\pi (N)}(\omega), \end{eqnarray} \tag{ E8 }$$

where the sum is over the N! permutations π of {1...N} and (− 1)^π denotes whether the permutation is odd or even. Inspection shows that if at zeroth order (in ω − ω_e, ω − ω_d, and δA) A(ω) is diagonal and has zero entries in the ee and dd slots, then $\det {\bf A}(\omega)$ is second order and the only terms at second order have π(h) = h for h∉{e, d}. Then $\det {\bf A}(\omega)$ is the product of the diagonal entries with h∉{e, d}, times the determinant of the 2 × 2 sub-block,

$$\begin{eqnarray} {\bf A}(\omega) \ni \left(\begin{array}{@{} cc @{}}(\omega -\omega _e)(\omega ^{\prime }_e-\omega) + \delta A_{ee} & \delta A_{ed} \\ \delta A_{de} & (\omega -\omega _d)(\omega ^{\prime }_d-\omega)+\delta A_{dd} \end{array}\right). \end{eqnarray} \tag{ E9 }$$

Moreover, at second order $\omega ^{\prime }_e-\omega$ in the above determinant may be replaced with $\omega ^{\prime }_e-\omega _e$. The latter gives a quadratic equation for ω:

$$\begin{eqnarray} &&0 = (\omega ^{\prime }_e-\omega _e)(\omega ^{\prime }_d-\omega _d)(\omega -\omega _e)(\omega -\omega _d) + (\omega ^{\prime }_d-\omega _d) \delta A_{ee} (\omega -\omega _d) + (\omega ^{\prime }_e-\omega _e) \delta A_{dd} (\omega -\omega _e) + \delta A_{ee}\delta A_{dd} - \delta A_{ed}^2, \end{eqnarray} \tag{ E10 }$$

and from the discriminant (treating the above equation as a quadratic in ω − ω_e) we get the criterion for an instability:

$$\begin{eqnarray} &&(4\varpi _e\varpi _d\Delta - 2\varpi _d \delta A_{ee} - 2 \varpi _e \delta A_{dd})^2 -16\varpi _e\varpi _d \left(\delta A_{ee}\delta A_{dd} - 2\varpi _d \delta A_{ee} \Delta - \delta A_{ed}^2 \right) <0, \end{eqnarray} \tag{ E11 }$$

where we have defined Δ = ω_e − ω_d and $2\varpi _e=\omega _e-\omega ^{\prime }_e$ (so that ϖ_e = ±ω_p for p-modes and ±ω_g for g-modes; in general ϖ denotes an inertial-frame unperturbed frequency). Algebraic simplification leads to

$$\begin{eqnarray} &&\left(\Delta + \frac{ \delta A_{ee} }{2\varpi _e} - \frac{ \delta A_{dd}}{2\varpi _d} \right)^2 + \frac{\delta A_{ed}^2}{\varpi _e\varpi _d} <0. \end{eqnarray} \tag{ E12 }$$

This is in fact the familiar criterion for parametric resonance. It requires first that ϖ_e and ϖ_d have opposite signs; if we take ϖ_e to be positive, then ω_e = ϖ_e + m_eΩ, ω_d = ϖ_d + m_dΩ, and hence the resonance criterion is |ϖ_e| + |ϖ_d| ≈ (m_d − m_e)Ω, thus it occurs only when the unperturbed frequencies sum to a harmonic of the orbital frequency. The coupling strength must be at least |δA_ed| > |ϖ_eϖ_d|^1/2|Δ|, with a correction to the detuning if the applied perturbation leads to first-order corrections to the mode frequencies (δA_ee and δA_dd).

This time we are interested in the case where two frequencies associated with the same mode are "near" each other and may merge—say ω_e and $\omega ^{\prime }_e$. Here the frequency difference is $2\varpi _e=\omega _e-\omega ^{\prime }_e$ and is treated as small; we have ω_e = m_eΩ + ϖ_e and $\omega ^{\prime }_e = m_e\Omega -\varpi _e$. We suppose that ω is near these frequencies, i.e., ω ≈ m_eΩ, and we work to second order in ϖ_e, ω − m_eΩ, and $\delta {\cal M}$.

Evaluation of $\det {\bf A}(\omega)$ is more subtle than the case of e ≠ d because many more terms in the determinant are important. In Equation (E8), we have two types of terms—those with π(e) = e and those with π(e) ≠ e. Since A_ee(ω) is at least first order, as are all off-diagonal entries, the only surviving term with π(e) = e is where π is the identity permutation (i.e., corresponding to the product of all diagonal entries in A). If π(e) = h ≠ e, then A_eh is first order, and so such a term can only survive if there is at most one other off-diagonal entry in the product, i.e., if π(h) = e and π(i) = i for all i∉{e, h}. These terms lead to the approximation

$$\begin{eqnarray} &&\det {\bf A}(\omega) \approx A_{ee}(\omega) \prod _{i\ne e} A_{ii}(\omega) - \sum _{h\ne e} \delta A_{eh}^2 \prod _{i\notin \lbrace e,h\rbrace } A_{ii}(\omega). \end{eqnarray} \tag{ E13 }$$

We now set this to zero and divide by ∏_{i ≠ e}A_ii(ω) to get

$$\begin{eqnarray} &&0 \approx A_{ee}(\omega) - \sum _{h\ne e} \frac{\delta A_{eh}^2}{A_{hh}(\omega)}. \end{eqnarray} \tag{ E14 }$$

Finally we see that A_hh(ω) can be approximated by its zeroth-order value (since it already appears multiplying a second-order perturbation), which is $(\omega -\omega _h)(\omega ^{\prime }_h-\omega) \approx (m_e\Omega -\omega _h)(\omega ^{\prime }_h- m_e\Omega)$. Also we substitute for A_ee(ω):

$$\begin{eqnarray} &&0 \approx \varpi _e^2 -(\omega - m_e\Omega)^2 + \delta A_{ee} - \sum _{h\ne e} \frac{\delta A_{eh}^2}{(m_e\Omega -\omega _h)(\omega ^{\prime }_h- m_e\Omega)}. \end{eqnarray} \tag{ E15 }$$

Finally, replacing ω_h and $\omega ^{\prime }_h$ by ±ϖ_h + m_hΩ, we see that the roots of this equation become complex—i.e., unstable—when

$$\begin{eqnarray} &&\varpi _e^2 + \delta A_{ee} - \sum _{h\ne e} \frac{\delta A_{eh}^2}{\varpi _h^2 - (m_e-m_h)^2\Omega ^2} = (\omega -m_e\Omega)^2 < 0. \end{eqnarray} \tag{ E16 }$$

This resembles the quasi-static instability criterion including both A_ee (which through second order includes both three-mode and four-mode couplings of mode e to the tidal field) and the square of A_eh (three-mode coupling of e and h to the tidal field), modified by the rotating reference frame term, (m_e − m_h)²Ω². In the case of interest, where e is a g-mode and there is a perturbation coupling to p-modes h, the instability criterion becomes

$$\begin{eqnarray} &&\omega _g^2 + \delta A_{{gg}} - \frac{1}{\omega _p^2} \sum _{h\; {\rm is}\; p{\rm -mode}} \delta A_{{g}h}^2 \left[ 1 - \frac{(m_{{g}}-m_h)^2\Omega ^2}{\omega _p^2} \right]^{-1} = (\omega -m_{{g}} \Omega)^2 < 0. \end{eqnarray} \tag{ E17 }$$

The last term can be expanded to leading order in $\Omega ^2/\omega _p^2$ to give

$$\begin{eqnarray} &&\omega _g^2 + \delta A_{{gg}} - \frac{1}{\omega _p^2} \sum _{h\; {\rm is}\; p{\rm -mode}} \delta A_{{g}h}^2 - \frac{\Omega ^2}{\omega _p^4} \sum _{h\; {\rm is}\; p{\rm -mode}} (m_{{g}}-m_h)^2\delta A_{{g}h}^2 = (\omega -m_{{g}}\Omega)^2 < 0. \end{eqnarray} \tag{ E18 }$$

We see that the instability criterion—that the left-hand side be negative—differs from the quasi-static case only in the introduction of a "centrifugal correction" of order

$$\begin{eqnarray} &&\frac{\Omega ^2 \delta A_{gp}^2}{\omega _p^4} \sim \left(\frac{\Omega {\cal M}_{gp}}{\omega _p^2} \right)^2 \sim \left(\frac{\Omega \epsilon \omega _p \omega _g J^{(1)}_{pg}}{\omega _p^2} \right)^2 \sim \left[ \frac{\Omega \epsilon \omega _p \omega _g (\omega _p/\omega _g)}{\omega _p^2} \right]^2\sim \epsilon ^2 \Omega ^2, \end{eqnarray} \tag{ E19 }$$

where we have used that $\delta A_{gp}\sim {\cal M}_{gp}$, used Equation (73) for ${\cal M}_{pg}$, and the value of the Jacobian. This correction goes in the direction of de-stabilizing the star. It has a simple interpretation in the language of the two-dimensional oscillator of Section 1.1: if the eigenvectors of ${\cal M}$ rotate at some rate $\dot{\theta }$, then a particle moving in the shallow direction in the potential well experiences a "centrifugal force" correction $-\dot{\theta }^2$ to $\omega _-^2$. In our case where the external tidal field is rotating, the shallow direction varies by an angle ${\cal O}(\epsilon)$ over a timescale Ω⁻¹, hence its rotation rate squared is ∼²Ω². The interpretation as such a term is clear since in Equation (E18), the "rotation angle" of mode e into mode h is $\delta A_{eh}/\omega _p^2$, and it oscillates at a rate (m_e − m_h)Ω. The second summation then represents the sum of squared angular rates.

For the highest-order g-modes with frequencies ω_g ≲ Ω, it is possible that the centrifugal term may dominate and lead to an instability whose growth timescale would be t_cen ∼ (Ω)⁻¹.