format: text

text:
# Nonlinear mode coupling in rotating stars and the 𝑟-mode instability in neutron stars

## I Introduction

Various authors have proposed that $`r`$-modes might be linearly unstable inside fairly rapidly rotating neutron stars, depending on a competition between gravitational radiation reaction, which drives the instability, and viscous effects, which inhibit it; several other investigators have explored implications for the spin evolution of rapidly rotating newborn and accreting neutron stars . An open question is how the instability saturates, and at what perturbation amplitude. This question is not only important from a theoretical standpoint, but also must be answered in order to assess whether or not the gravitational radiation emitted during the development of the instability is detectable. The saturation amplitude also determines the final spin rate of the neutron star when it exits the instability region of parameter space. While turbulence generated in the strong shear layer at the interface between the fluid core and solid crust has been shown to saturate the $`r`$-mode , a crust will not be present initially for newly born neutron stars, and other hydrodynamical saturation mechanisms are still being investigated and may in fact be more important.

There have been, to date, two fully nonlinear numerical calculations of the development of the instability, but neither completely settles the question of saturation, principally because there are practical limitations on what it is feasible to compute. The first calculation, by Stergioulas and Font , considered the evolution of the fluid inside a neutron star with no buoyancy force and fixed spacetime geometry. Their calculations began with a large amplitude excitation in the mode expected to be most unstable, the $`l=m=2`$ $`r`$-mode, and followed the subsequent hydrodynamics for about 20 stellar rotation periods, approximately 25 ms for their calculations. During this time no substantial change in the amplitude of the modal excitation was seen, and no evidence for significant nonlinear excitation of other stellar normal modes was detected. Lindblom, Tohline and Vallisneri  subsequently performed simulations based on the Newtonian equations of hydrodynamics and gravitation, including an artificially enhanced gravitational radiation force consisting of the correct force  multiplied by a factor of about 4500 . The linear instability growth timescale in their simulations was about 13 rotation periods, rather than the correct value of $`10^4`$ rotation periods. The initial $`r`$mode amplitude chosen by Lindblom, Tohline and Vallisneri was small (about 0.1), and, as expected, grew via the radiation-driven instability to large amplitude, where shocks ultimately formed, accompanied by a decrease in $`r`$mode amplitude. Both calculations show that nonlinear mode-mode coupling is comparatively weak for $`r`$modes, producing little effect even at large modal amplitude on timescales $`1020`$ stellar rotation periods.

A second, more analytical approach to this problem is possible. One assumes, at the outset, that the modes develop only moderately nonlinear amplitudes, and nonlinear effects are modeled perturbatively, via $`n`$mode couplings. When amplitudes are small, relatively low order couplings, such as three-mode interactions, suffice, and (when dissipation and radiation reaction are ignored) the equations of motion for the fluid can be modeled as an infinite-dimensional Hamiltonian system with a polynomial interaction potential (which is cubic when only three-mode couplings are kept). This approach has the advantage of allowing explicit calculation of the strength of the lowest order coupling of a given $`r`$mode to other modes of the star, but has the disadvantage of failing when mode amplitudes grow very large. Although it is possible that the $`r`$mode instability saturates only in the fully nonlinear regime, neither of the numerical calculations establish this, because it is possible that modal couplings act on timescales longer than covered by the simulations, but shorter than the growth time ($`10^4`$ rotation periods) of the gravitational-radiation instability. The analytic approach we are taking should apply if the instability saturates before attaining very large amplitudes. Moreover, by computing the three-mode coupling coefficients for $`r`$modes, we should be able to shed some light on why nonlinear coupling is so inefficient for them on dynamical timescales in the simulations.

In this paper, we develop the tools needed to calculate the nonlinear evolution of the unstable modes in the weakly nonlinear regime. We concentrate on three-mode interactions for the most part, although some of the results we derive could also apply to interactions involving $`n>3`$ modes of the star.

We begin with a review of the linear perturbation theory in Sec. II. A principal goal of that section of the paper is to obtain equations of motion for modes of a rotating star acted on by an external force that are uncoupled from one another at linear order. The method we develop is based on a series of papers by Schutz and collaborators , but has not been applied consistently to $`r`$modes before. The key results, for a reader wishing to skip the details, are the mode decomposition formulae (74) and (75) and the equation of motion for each mode (76). In Sec. III, we review the slow rotation approximation to mode functions and frequencies, with particular attention to the rotational modes, whose frequencies are zero in the limit of zero rotation frequency, $`\mathrm{\Omega }`$. The key new result in this section is the fact that the rotational mode functions are orthogonal to leading order in the stellar angular velocity.

In Sec. IV, we discuss the second order Lagrangian perturbation theory at the heart of our calculation of the evolution of modes in the weakly nonlinear regime. The key results of this section are (i) the coupled nonlinear equations of motion (128) for the mode coefficients, and (ii) the formula (134) for the three-mode coupling coefficients that can be evaluated once the modes have been found for a given background star model. In Sec. V, we discuss properties of the coupling coefficients. We show how general parity arguments and selection rules can be used to explain why some modes do not couple to one another, and also show that there are other couplings that are suppressed, in that they may vanish to lowest order in $`\mathrm{\Omega }`$, or to lowest order in compressibility, but are nonzero more generally. In Sec. VI, we develop a new and efficient way of computing coupling coefficients that employs spin-weighted spherical harmonics . The method greatly simplifies the calculation of the coupling coefficients for nonrotating stars as well as rotating stars, and should be useful for other astrophysical applications (see, e.g., Refs. ).

This paper concentrates on formal issues primarily, postponing numerical integration of the time-dependent evolution of a network of weakly interacting modes to a subsequent publication. From a formal viewpoint, the most substantial result of this paper is a formulation of the equations of motion that, at linear order, separates the responses of individual modes to an external force, and, at the first nonlinear order, involves coupling coefficients that we can calculate relatively economically. From an astrophysical viewpoint, the most interesting result is that many interactions of $`r`$modes with themselves can be shown to be either forbidden entirely by selection rules, or suppressed to lowest order in some parameter, such as rotation frequency or compressibility. By itself, this is not quite enough to explain the smallness of the coupling seen in numerical simulations, but we may already be seeing a hint of why the transfer of energy from $`r`$modes to other modes is so inefficient on dynamical timescales.

## II Linear Perturbation Theory for Rotating Stars: The Schutz Formalism

### A Overview

In this section we develop the formalism we need in the linear Lagrangian perturbation theory of rotating stars. There are two aspects to the perturbation theory that will be essential for our purposes:

(i) Solving for the frequencies $`\omega _\alpha `$ and mode functions $`𝝃_\alpha (𝐱)`$. This aspect is well understood and highly developed; see for example Ch. VI of Unno et al , as well as Refs. .

(ii) Decomposing general motions in a star into sums over modes, and obtaining equations of motion for the mode coefficients that are uncoupled at linear order when an external force acts on the star.

The second aspect of the theory — deriving equations of motion for the mode coefficients — has not been well-understood in the past. The necessary theory is implicit in a series of papers by Schutz and collaborators , but the formalism has not been applied widely in astrophysics. For example, a classic computation for which modal equations of motion would be useful is the excitation of the modes of a rotating star by the tidal field of a binary companion. While there have been many papers on this topic, we are not aware of any that have used the correct uncoupled equations of motion. For example, Refs.  and  use a type of mode expansion \[Eq. (77) below\] for which the equations of motion of the various mode amplitudes are coupled, and Refs. ,  and  use the same mode expansion but simply drop the coupling terms between the different modes. Some authors have resorted to numerically solving the linearized hydrodynamic equations in the frequency domain, after factoring out a factor of $`e^{im\phi i\omega t}`$, rather than using a mode expansion .

This section of the paper is devoted to presenting and explaining the equations of motion that result from Schutz’s formalism when one makes a number of simplifying assumptions that usually are valid in applications. We have attempted to make the presentation transparent so as to be accessible to a wide audience. All the derivations are relegated to an appendix. We start in Sec. II B by summarizing the governing equations. In Sec. II C we describe the form of the mode expansion and equations of motion for the mode coefficients.

### B Governing equations

We assume that the unperturbed background star is uniformly rotating with angular velocity $`𝛀`$. The fluid equations of motion in the corotating frame are

$$\frac{\rho }{t}+(\rho 𝐮)=0$$
(1)
and

$$\frac{𝐮}{t}+(𝐮)𝐮+2𝛀\times 𝐮+𝛀\times \left(𝛀\times 𝐱\right)=\frac{p}{\rho }\varphi +𝐚_{\mathrm{ext}},$$
(2)
where $`\rho `$ is the density, $`p`$ the pressure, $`𝐮`$ the velocity and $`𝐚_{\mathrm{ext}}`$ is any acceleration due to external forces. The gravitational potential $`\varphi `$ is given by $`^2\varphi =4\pi G\rho `$. In the background solution, $`𝐮`$ vanishes and $`\rho `$ is time independent.

Now consider linearized perturbations characterized by Eulerian perturbations $`\delta 𝐮`$, $`\delta \rho `$ and $`\delta p`$. The first two of these are related to the linearized Lagrangian displacement $`𝝃(𝐱,t)`$ by

$$\delta \rho =(\rho 𝝃)$$
(3)
and $`\delta 𝐮=\dot{𝝃}`$ (since $`𝐮=0`$). We assume that the perturbations are characterized by an adiabatic index $`\mathrm{\Gamma }_1`$, so that

$$\frac{\mathrm{\Delta }p}{p}=\mathrm{\Gamma }_1\frac{\mathrm{\Delta }\rho }{\rho },$$
(4)
where $`\mathrm{\Delta }\rho =\delta \rho +(𝝃\mathbf{})\rho `$ and $`\mathrm{\Delta }p=\delta p+(𝝃\mathbf{})p`$ are the Lagrangian perturbations of density and pressure***If the background star and the perturbations both obey the same one-parameter equation of state $`p=p(\rho )`$, then $`\mathrm{\Gamma }_1=(\rho /p)dp/d\rho `$ depends only on the background density $`\rho `$. More generally, we can regard $`\mathrm{\Gamma }_1=\mathrm{\Gamma }_1(𝐱)`$ as a function of position determined by the structure of the background star. A very general definition of the adiabatic index $`\mathrm{\Gamma }_1=\mathrm{\Gamma }_1(𝐱)`$ that is consistent with the phenomenological relation (4) is given in Eq. (J27) below.. By combining these relations with linearized versions of the continuity and Euler equations (1) and (2),

$$\delta \dot{\rho }+(\rho \delta 𝐮)=0,$$
(5)
and

$$\delta \dot{𝐮}+2𝛀\times \delta 𝐮=\frac{\delta p}{\rho }+\frac{p}{\rho ^2}\delta \rho \delta \varphi +𝐚_{\mathrm{ext}},$$
(6)
where $`^2\delta \varphi =4\pi G\delta \rho `$, we obtain the equations of motion for linearized perturbations in the form

$`\ddot{𝝃}+𝐁\dot{𝝃}+𝐂𝝃=𝐚_{\mathrm{ext}}(𝐱,t).`$ (7)
Here the operator $`𝐁`$ is given by

$`𝐁𝝃2𝛀\times 𝝃,`$ (8)
and the operator $`𝐂`$ is given by

$`\rho \left(𝐂𝝃\right)_i`$ $`=`$ $`_i(\mathrm{\Gamma }_1p_j\xi ^j)+_ip_j\xi ^j_jp_i\xi ^j`$ (10)
$`+\rho \xi ^j_j_i\varphi +\rho \xi ^j_j_i\varphi _{\mathrm{rot}}+\rho _i\delta \varphi ,`$
where

$$\varphi _{\mathrm{rot}}(𝐱)=\frac{1}{2}(𝛀\times 𝐱)^2.$$
(11)
The $`𝐂`$ operator can be broken up into two pieces

$`𝐂=𝐂_a+𝐂_b`$ (12)
where $`𝐂_a`$ is proportional to the vectorial Schwarzschild discriminant

$$𝐀=\frac{\mathbf{}\rho }{\rho }\frac{1}{\mathrm{\Gamma }_1}\frac{\mathbf{}p}{p},$$
(13)
and $`𝐂_b`$ is independent of $`𝐀`$. We find

$$𝐂_a𝝃=\mathbf{}\left(\frac{\mathrm{\Gamma }_1p}{\rho }\frac{\delta \rho }{\rho }+\delta \varphi \right)$$
(14)
and

$`𝐂_b𝝃`$ $`=`$ $`\mathbf{}\left({\displaystyle \frac{\mathrm{\Gamma }_1p}{\rho }}𝝃𝐀\right)+{\displaystyle \frac{\delta p}{\rho }}𝐀+{\displaystyle \frac{\mathbf{}p}{\rho }}𝝃𝐀`$ (15)
$`=`$ $`{\displaystyle \frac{p}{\rho ^2}}\delta \rho \mathrm{\Gamma }_1𝐀+{\displaystyle \frac{1}{\rho }}\left(\mathrm{\Gamma }_1p𝝃𝐀\right).`$ (16)
The buoyancy force described by the operator $`𝐂_b`$ gives rise to $`g`$ modes whose frequency scale is set by the Brunt-Väisälä frequency $`N`$ (Brunt for short), defined by

$$N^2=\frac{p}{\rho }𝐀.$$
(17)
For much of this paper we will specialize to perturbations of stars without buoyancy, i.e., situations where $`𝐀=0`$For such stars $`dpd\rho =0`$ and so the background star satisfies $`p=p(\rho )`$ for some equation of state. Furthermore the perturbations obey the same equation of state, so the star is barotropic. \[A special case of this is a star of uniform specific entropy, i.e., an isentropic star.\] For more general stars $`dp`$ and $`d\rho `$ need not be proportional.. For example, if the background star and the perturbations obey one parameter equations of state with adiabatic indices $`\mathrm{\Gamma }`$ and $`\mathrm{\Gamma }_1`$, respectively, then $`𝐀\mathrm{\Gamma }\mathrm{\Gamma }_1`$, and so the zero-buoyancy case is $`\mathrm{\Gamma }=\mathrm{\Gamma }_1`$. For ordinary stars, the Brunt frequency in the convection zone is approximately zero when the superadiabatic gradient is small. For neutron stars, $`𝐀=0`$ means a cold, zero-entropy gas in which the composition adjusts instantly to changes in density. In real neutron stars the adjustment speed is limited by beta and inverse-beta decays, which gives rise to $`g`$ modes  (see Appendix J for a related discussion). Our analysis of $`r`$ modes will make the zero-buoyancy approximation $`𝐀=0`$, which is expected to be good in the regime $`N\mathrm{\Omega }`$. However, this inequality will only be marginally satisfied in newly born neutron stars, so it will be important for future analyses to analyze nonlinear couplings of $`r`$ modes in nonzero-buoyancy stars. We choose to focus on the zero-buoyancy case for simplicity.

The dynamical equation (7) forms the starting point for our discussions. We can find a large class of solutions for the case $`𝐚_{\mathrm{ext}}=0`$ of no forcing if we make the ansatz

$`𝝃(𝐱,t)=e^{i\omega t}𝝃(𝐱),`$ (18)
where $`\omega `$ is the rotating-frame frequency. This yields

$`\left[\omega ^2i\omega 𝐁+𝐂\right]𝝃=0.`$ (19)
The quadratic eigenvalue equation (19) is the standard equation that one solves to obtain the eigenfunctions $`𝝃(𝐱)`$ and eigenfrequencies $`\omega `$ for rotating stars. A mode of the star will consist of a pair $`(𝝃,\omega )`$ that satisfy Eq. (19).

### C Mode decomposition formalism

#### 1 Nonrotating stars

We start by recalling the standard mode decomposition formalism for non-rotating stars for which $`𝐁=0`$, in order to contrast it with the rotating case below. For non-rotating stars Eq. (19) reduces to

$$𝐂𝝃=\omega ^2𝝃,$$
(20)
which represents a standard eigenvalue problem for $`\omega ^2`$. We define the inner product on the space $``$ of complex vector functions $`𝝃=𝝃(𝐱)`$ by

$`𝝃,𝝃^{}={\displaystyle d^3x\rho (𝐱)𝝃(𝐱)^{}𝝃^{}(𝐱)}.`$ (21)
The operator $`𝐂`$ is Hermitian with respect to this inner productThroughout this paper we assume that the space $``$ is finite dimensional, as it will be in any numerical computations, to avoid discussing subtleties related to infinite dimensional spaces.: for any elements $`𝝃`$ and $`𝝃^{}`$ of $``$,

$`𝝃,𝐂𝝃^{}=𝐂𝝃,𝝃^{},`$ (22)
while $`𝐁`$ is anti-Hermitian

$`𝝃,𝐁𝝃^{}=𝐁𝝃,𝝃^{}.`$ (23)
It follows that one can find a set $`\{𝝃_\alpha \}`$ of eigenvectors of $`𝐂`$ with eigenvalues $`\omega _\alpha ^2`$ which are orthonormal

$$𝝃_\alpha ,𝝃_\beta =\delta _{\alpha \beta }$$
(24)
and which also form a basis for the vector space $``$. Hence we can decompose any Lagrangian displacement $`𝝃(𝐱,t)`$ as §§§ The sum over $`\alpha `$ in Eq. (25) is over all distinct solutions $`𝝃_\alpha `$ to Eq. (20), not over the larger set of distinct pairs $`(𝝃,\omega )`$, as there is only one term in the sum (25) corresponding to the two solutions $`(𝝃_\alpha ,\omega _\alpha )`$ and $`(𝝃_\alpha ,\omega _\alpha )`$. This will be important below.

$$𝝃(𝐱,t)=\underset{\alpha }{}q_\alpha (t)𝝃_\alpha (𝐱),$$
(25)
where

$$q_\alpha (t)=𝝃_\alpha ,𝝃(t).$$
(26)
Finally, by inserting the expansion (25) into the dynamical equation (7) with $`𝐁=0`$ and by contracting with $`𝝃_\alpha `$, one obtains the equation of motion

$$\ddot{q}_\alpha (t)+\omega _\alpha ^2q_\alpha (t)=𝝃_\alpha ,𝐚_{\mathrm{ext}}(t),$$
(27)
which is just the forced harmonic oscillator equation.

#### 2 Rotating stars

Turn, now, to the corresponding formalism for rotating stars. We make a number of simplifying assumptions:

* The background star is stable, so all frequencies are real.
* All the modes are non-degenerate: for each frequency $`\omega `$ the dimension of the space of eigenmodes is one. This assumption will be relaxed in Sec. II C 3 below.

A more general treatment, without these assumptions, is given in Appendix A.

The eigenvalue equation (19) can be reformulated as a standard eigenvalue equation \[Eq. (A15) below\] for $`\omega `$ involving a non-Hermitian operator. Non-Hermitian operators in general do not posess sufficiently many right eigenvectors to form a basis. However, for such operators there is a standard procedure described in linear algebra textbooks to obtain a basis by adding to the eigenvectors certain additional vectors, sometimes called generalized eigenvectors or Jordan chain vectors. Each actual eigenvector may have one or more generalized eigenvectors associated with it, forming a so-called Jordan chain.

There are two classes of modes $`(𝝃,\omega )`$ in a rotating star: modes associated with non-trivial Jordan chains , which we shall call Jordan-chain modes, and modes not associated with such Jordan chains (i.e., the length of the chain is zero; see Appendix A for more details). A mode that is just at the point of becoming unstable is always a Jordan-chain mode . However, even for stable rotating stars, Jordan-chain modes of zero frequency are always present. One set of such modes corresponds to displacements of the star to nearby equilibria with slightly different angular velocities

$$\mathrm{\Omega }\mathrm{\Omega }+\delta \mathrm{\Omega }(r_{}),$$
(28)
where $`r_{}`$ is the distance from the rotation axis (see Appendix D for a proof of this in the zero-buoyancy case)There are also Jordan chain modes corresponding to tilting the axis of rotation of the star.. Therefore, one cannot assume that no Jordan chains occur. In this section we shall nevertheless, for simplicity, describe the formalism that would apply if there were no Jordan-chain modes. This formalism should be useful, for example, in describing situations in which the Jordan-chain modes are unimportant dynamically. In the remainder of this paper we shall assume that the Jordan chain modes are unimportant for the process of saturation of $`r`$-modesHowever, note that several studies have found that the gravitational radiation reaction and/or hydrodynamic nonlinearities induce significant differential rotation in the star . A set of coupled equations for the growth of the $`r`$-mode as well as for the differential rotation can be formulated using the results of Appendix A. Since the differential rotation corresponds to a set of Jordan chains of length one (Appendix D), Eqs. (A29) and (A30) can be used to evolve the differential rotation in the linear regime; however the linear approximation breaks down after a time $`1/\delta \mathrm{\Omega }`$.. A more complete treatment of the formalism, allowing Jordan-chain modes, is detailed in Appendix A.

Given these assumptions, there are two obstacles to obtaining a mode decomposition formalism similar to that for non-rotating stars. First, distinct modes $`(𝝃,\omega )`$ and $`(𝝃^{},\omega ^{})`$ satisfying Eq. (19) with $`\omega \omega ^{}`$ need not be orthogonal with respect to the inner product (21) . Although the operator

$$𝐋(\omega )=\omega ^2i\omega 𝐁+𝐂$$
(29)
is Hermitian for real $`\omega `$ , the vectors $`𝝃`$ and $`𝝃^{}`$ are eigenvectors of the two different operators $`𝐋(\omega )`$ and $`𝐋(\omega ^{})`$ and so need not be orthogonal. Nevertheless, it is possible to find a set $`\{𝝃_\alpha \}`$ of solutions to Eq. (19) that form a basis of $``$, so that every Lagrangian displacement can be uniquely decomposed as a superposition of the form

$$𝝃(𝐱,t)=\underset{\alpha }{}q_\alpha (t)𝝃_\alpha (𝐱).$$
(30)
The second obstacle is the following: There is in general no choice of basis $`\{𝝃_\alpha \}`$ of eigenvectors for which the equations of motion of the coefficients $`q_\alpha (t)`$ defined by Eq. (30) are uncoupled from one another.

To circumvent this obstacle, it is necessary to use a phase space mode expansion instead of a configuration space mode expansion, as pointed out by Dyson and Schutz . Let’s label the distinct non-Jordan-chain solutions\**\**\**Here we define two solutions $`(𝝃,\omega )`$ and $`(𝝃^{},\omega ^{})`$ to be distinct if $`\omega \omega ^{}`$ or if $`𝝃`$ and $`𝝃^{}`$ are linearly independent. $`(𝝃,\omega )`$ to Eq. (19) as $`(𝝃_A,\omega _A)`$, so that

$$\left[\omega _A^2i\omega _A𝐁+𝐂\right]𝝃_A=0.$$
(31)
The set of vectors

$`\left[\begin{array}{c}𝝃_A\\ i\omega _A𝝃_A\end{array}\right]`$ (34)
forms a basis for the space $``$ of pairs of vectors $`[𝝃,𝝃^{}]`$. We deliberately use a different type of index here — capital Roman indices rather than lower case Greek indices — as a reminder that the number of distinct $`A`$’s is twice the dimension of $``$, whereas the number of distinct $`\alpha `$’s in the sums (25) and (30) is just the dimension of $``$.

The mode functions $`𝝃_A`$ are not orthogonal in general, that is, $`𝝃_A,𝝃_B`$ need not vanish when $`AB`$. We give an explicit example demonstrating the non-orthogonality in Sec. III below. However, as shown by Ref. , the modes $`(𝝃_A,\omega _A)`$ do obey a modified type of orthogonality relation, which is ††††††The left hand side of Eq. (35) is proportional to the symplectic product $`W(𝝃_1,𝝃_2)`$ defined by Ref. , where $`𝝃_1(𝐱,t)=𝝃_A(𝐱)\mathrm{exp}[i\omega _At]`$ and $`𝝃_2(𝐱,t)=𝝃_B(𝐱)\mathrm{exp}[i\omega _Bt]`$. An alternative proof of the orthogonality relation (35) from a phase space construction was given by Ref. , which we review in Appendix A.

$$𝝃_A,i𝐁𝝃_B+(\omega _A+\omega _B)𝝃_A,𝝃_B=0$$
(35)
for $`AB`$.

Now, at any time $`t`$, we can form the vector

$`\left[\begin{array}{c}𝝃(t)\\ \dot{𝝃}(t)\end{array}\right]`$ (38)
from the Lagrangian displacement $`𝝃(𝐱,t)`$ and its time derivative $`\dot{𝝃}(𝐱,t)`$. We can expand this vector in the basis (34) as‡‡‡‡‡‡Note that it follows from Eq. (43) that
$$\underset{A}{}(\dot{c}_A+i\omega _Ac_A)𝝃_A=0.$$
However, it does not follow that $`\dot{c}_A+i\omega _Ac_A=0`$ for each $`A`$ since the set $`\{𝝃_A\}`$ is not a basis for $``$; it has twice too many basis elements.

$`\left[\begin{array}{c}𝝃(t)\\ \dot{𝝃}(t)\end{array}\right]={\displaystyle \underset{A}{}}c_A(t)\left[\begin{array}{c}𝝃_A\\ i\omega _A𝝃_A\end{array}\right].`$ (43)
Using the orthogonality relation (35), the expansion (43) can be inverted (see Appendix A) to obtain

$$c_A(t)=\frac{1}{b_A}𝝃_A,\omega _A𝝃(t)+i\dot{𝝃}(t)+i𝐁𝝃(t),$$
(44)
where the constant $`b_A`$ is given by \**\**\**Note that adjusting the normalization of $`𝝃_A`$ will change the magnitude but not the sign of $`b_A`$. Hence we cannot find a basis that achieves $`b_A=1`$ for all $`A`$, since some of the $`b_A`$’s are negative and some are positive.

$$b_A=𝝃_A,i𝐁𝝃_A+2\omega _A𝝃_A,𝝃_A.$$
(45)
Since the operator $`i𝐁`$ is Hermitian, $`b_A`$ is real although it need not be positive. In Appendix A we show that for $`\omega _A0`$, the constant $`b_A`$ is related to the rotating-frame mode energy $`\epsilon _A`$ at unit amplitude by $`b_A=\epsilon _A/\omega _A`$. Finally, the equation of motion for the mode coefficient $`c_A(t)`$ is the first-order-in-time equation

$$\dot{c}_A(t)+i\omega _Ac_A(t)=\frac{i}{b_A}𝝃_A,𝐚_{\mathrm{ext}}(t).$$
(46)
Using the mode decomposition (43), its inverse (44), and the equation of motion for each mode (46), it is straightforward to compute the response of the star to any externally applied acceleration $`𝐚_{\mathrm{ext}}(𝐱,t)`$. See Appendix A for justification of all the claims and formulae of this subsection.

#### 3 Degeneracy

Consider next the situation where there is degeneracy, that is, where distinct modes $`(𝝃_A,\omega _A)`$ and $`(𝝃_B,\omega _B)`$ have identical frequencies $`\omega _A=\omega _B`$. We introduce the following index notation: we write the distinct eigenfrequencies as $`\omega _a`$, and we let the associated eigenvectors be $`𝝃_{a,k}`$ for $`1kn_a`$, where $`n_a`$ is the degeneracy associated with $`\omega _a`$. Thus, we identify the index $`A`$ with the pair of indices $`(a,k)`$, and we make the identifications

$`\omega _A`$ $`=`$ $`\omega _{a,k}=\omega _a,`$ (47)
$`𝝃_A`$ $`=`$ $`𝝃_{a,k}.`$ (48)
We can also rewrite sums over $`A`$ as

$$\underset{A}{}=\underset{a}{}\underset{k=1}{\overset{n_a}{}}.$$
(49)
Consider now the matrix

$$_{AB}𝝃_A,i𝐁𝝃_B+(\omega _A+\omega _B)𝝃_A,𝝃_B.$$
(50)
Using the notation $`A=(a,k)`$ and $`B=(b,l)`$, we can write the matrix as $`_{AB}=_{ak,bl}`$. In Appendix A we show that the matrix $``$ is always block diagonal in the sense that $`_{ak,bl}=0`$ for $`ab`$. In other words, for $`\omega _A\omega _B`$ we have $`_{AB}=0`$, which is a generalization of the orthogonality relation (35) above.

Within a given degenerate subspace, however, we have

$$_{ak,al}=𝝃_{a,k},i𝐁𝝃_{a,l}+2\omega _a𝝃_{a,k},𝝃_{a,l},$$
(51)
which is an $`n_a\times n_a`$ matrix that depends on the choice of basis $`𝝃_{a,k}`$. Since we are assuming in this section that all the eigenfrequencies are real, the matrix (51) is Hermitian and can be diagonalized by an appropriate chance of basis $`𝝃_{a,k}T_k^l𝝃_{a,l}`$ within the degenerate subspace. A choice of basis $`𝝃_{a,k}`$ that diagonalizes the matrix (51) is the analog for rotating stars of orthonormal bases for non-rotating stars. For such bases, the formulae of Sec. II C 2 above are valid in the degenerate as well as the non-degenerate cases, namely the mode decomposition (43), its inverse (44), the definition (45) of the constant $`b_A`$, and the equation of motion (46).

#### 4 Reality conditions

So far, the formalism is valid for complex Lagrangian displacements $`𝝃(𝐱,t)`$. We now describe simplifications that occur in the physically relevant case of real $`𝝃(𝐱,t)`$.

If $`(𝝃_A,\omega _A)`$ is a solution of the quadratic eigenvalue equation (19), then another solution $`(𝝃_A^{},\omega _A)`$ can be obtained by complex-conjugating the mode function and reversing the sign of the frequency . Under this transformation, the normalization constant $`b_A`$ flips sign, from Eq. (45) \*†\*†\*†The inertial-frame frequency $`\omega +m\mathrm{\Omega }`$, where $`m`$ is the azimuthal quantum number, also flips sign under this transformation.. In Appendix A we show that all non Jordan-chain modes have $`b_A0`$. Hence, all the distinct non Jordan-chain modes occur in pairs $`(𝝃,\omega )`$ and $`(𝝃^{},\omega )`$. \*‡\*‡\*‡If a mode satisfies (i) $`𝝃_A`$ is purely real, and (ii) $`\omega _A=0`$, then from Eqs. (45) and (23) that mode also has $`b_A=0`$ and hence is a Jordan-chain mode.

Let us now focus attention on the set of modes with $`b_A>0`$ For modes with $`\omega _A0`$, we have
$$b_A=\omega _A𝝃_A,𝝃_A+\frac{1}{\omega _A}𝝃_A,𝐂𝝃_A,$$
from Eqs. (19) and (45). Therefore positive frequency $`\omega _A>0`$ would correspond to $`b_A>0`$ if the operator $`𝐂`$ were positive. However, $`𝐂`$ need not always be a positive operator. Nevertheless we suspect that for modes of nonzero frequency, positive frequency always corresponds to $`b_A>0`$, although we have been unable to prove this.. We introduce the following index notation: we write the distinct modes with $`b_A>0`$ as $`(𝝃_\alpha ,\omega _\alpha )`$. We identify the index $`A`$ with the pair of indices $`(\alpha ,ϵ)`$, where $`ϵ`$ takes on the values $`ϵ=+`$ and $`ϵ=`$, and we make the identifications $`𝝃_A=𝝃_{\alpha ,ϵ}`$, $`\omega _A=\omega _{\alpha ,ϵ}`$, and $`b_A=b_{\alpha ,ϵ}`$ with

$`\omega _{\alpha ,+}`$ $`=`$ $`\omega _\alpha `$ (52)
$`\omega _{\alpha ,}`$ $`=`$ $`\omega _\alpha `$ (53)
$`𝝃_{\alpha ,+}`$ $`=`$ $`𝝃_\alpha `$ (54)
$`𝝃_{\alpha ,}`$ $`=`$ $`𝝃_\alpha ^{}`$ (55)
$`b_{\alpha ,+}`$ $`=`$ $`b_\alpha `$ (56)
$`b_{\alpha ,}`$ $`=`$ $`b_\alpha .`$ (57)
We can also rewrite sums over $`A`$ as

$$\underset{A}{}=\underset{\alpha }{}\underset{ϵ=\pm }{}.$$
(58)
Using these notations, the mode expansion (43) can be written as

$`\left[\begin{array}{c}𝝃(t)\\ \dot{𝝃}(t)\end{array}\right]={\displaystyle \underset{\alpha }{}}c_{\alpha ,+}(t)\left[\begin{array}{c}𝝃_\alpha \\ i\omega _\alpha 𝝃_\alpha \end{array}\right]+c_{\alpha ,}(t)\left[\begin{array}{c}𝝃_\alpha ^{}\\ i\omega _\alpha 𝝃_\alpha ^{}\end{array}\right],`$ (65)
and its inverse (44) can be written as

$$c_{\alpha ,+}(t)=\frac{1}{b_\alpha }𝝃_\alpha ,\omega _\alpha 𝝃(t)+i\dot{𝝃}(t)+i𝐁𝝃(t),$$
(66)
and

$$c_{\alpha ,}(t)=\frac{1}{b_\alpha }𝝃_\alpha ^{},\omega _\alpha 𝝃(t)i\dot{𝝃}(t)i𝐁𝝃(t).$$
(67)
It follows that $`𝝃(t)`$ will be real if and only if $`c_{\alpha ,}(t)=c_{\alpha ,+}(t)^{}`$ for all $`\alpha `$. In this case, writing $`c_\alpha c_{\alpha ,+}`$, the mode expansion (65) becomes

$`\left[\begin{array}{c}𝝃(t)\\ \dot{𝝃}(t)\end{array}\right]={\displaystyle \underset{\alpha }{}}c_\alpha (t)\left[\begin{array}{c}𝝃_\alpha \\ i\omega _\alpha 𝝃_\alpha \end{array}\right]+c_\alpha (t)^{}\left[\begin{array}{c}𝝃_\alpha ^{}\\ i\omega _\alpha 𝝃_\alpha ^{}\end{array}\right],`$ (74)
the inverse mode expansion (66) becomes

$$c_\alpha (t)=\frac{1}{b_\alpha }𝝃_\alpha ,\omega _\alpha 𝝃(t)+i\dot{𝝃}(t)+i𝐁𝝃(t),$$
(75)
and the equation of motion (46) becomes

$$\dot{c}_\alpha (t)+i\omega _\alpha c_\alpha (t)=\frac{i}{b_\alpha }𝝃_\alpha ,𝐚_{\mathrm{ext}}(t).$$
(76)
Note that the number of distinct $`𝝃_\alpha `$’s is the dimension of $``$. In Appendix B we show that the $`𝝃_\alpha `$’s are linearly independent and thus form a basis of $``$ More precisely, they would form a basis for $``$ if there were no Jordan-chain modes.. Hence, one could perform a configuration-space mode expansion of the form

$$𝝃(𝐱,t)=\underset{\alpha }{}q_\alpha (t)𝝃_\alpha (𝐱),$$
(77)
instead of the phase space mode expansion (74). However, the coefficients $`q_\alpha (t)`$ defined by Eq. (77) are not related in any simple way to the coefficients $`c_\alpha (t)`$ \*∥\*∥\*∥For the special case of a non-rotating star and real $`𝝃(𝐱,t)`$, there is a simple relation between the two sets of coefficients, namely $`2c_{\alpha ,+}=q_\alpha +i\dot{q}_\alpha /\omega _\alpha `$ and $`2c_{\alpha ,}=q_\alpha ^{}i\dot{q}_\alpha ^{}/\omega _\alpha `$., and furthermore the equations of motion for the coefficients $`q_\alpha `$ will not be uncoupled from one another.

In Appendix C we give an alternative form of the equations of motion, valid for arbitrary complex $`𝝃(𝐱,t)`$, which combines the two first order equations for $`c_{\alpha ,+}(t)`$ and $`c_{\alpha ,}(t)`$ in a single second order equation. That alternative form, when specialized to a non-rotating star and to a real basis of eigenmodes $`𝝃_\alpha (𝐱)`$, coincides with the standard formalism (27) for a non-rotating star.

Finally, in Appendix K we give expressions for the energy and angular momentum of perturbations in terms of the mode amplitudes.

## III The Slow Rotation Expansion

### A Overview

In this section we discuss the approximation method of solving for the mode functions $`𝝃_A(𝐱)`$ and frequencies $`\omega _A`$ by using a perturbative expansion in the star’s angular velocity $`\mathrm{\Omega }`$. The main result of this section is that the leading order (in $`\mathrm{\Omega }`$) mode functions for the inertial or hybrid modes in a slowly rotating zero-buoyancy star are orthogonal with respect to the inner product (21) of the non-rotating star.

There are two general classes of modes in rotating stars :

* Modes for which the frequency goes to zero in the non-rotating limit. These modes can be called rotational since they have non-zero frequency only in the presence of rotation. For zero-buoyancy stars, these modes are the hybrid modes of Ref. , also called inertial modes. For nonzero-buoyancy stars, they are purely axial.
* Modes with a finite frequency in the non-rotating limit, such as $`f`$ and $`p`$ modes, and $`g`$ modes in nonzero-buoyancy stars.

There is a substantial literature on computing mode functions and frequencies perturbatively in powers of $`\mathrm{\Omega }`$. For nonzero-buoyancy stars, the rotational or Rossby modes have been computed to leading order in $`\mathrm{\Omega }`$ by Refs. ,  and , and to all orders in $`\mathrm{\Omega }`$ by Ref. . Non-rotational modes in the nonzero-buoyancy case have been computed by Ref. , and by Ref.  beyond the leading order in $`\mathrm{\Omega }`$.

For zero-buoyancy stars, there are two classes of rotational modes. The so-called pure or classical $`r`$-modes  exist only for $`l=|m|`$ and are purely axial. The remaining rotational modes have both polar and axial pieces and have been called generalized $`r`$-modes, hybrid modes or inertial modes. They have been obtained to leading order in $`\mathrm{\Omega }`$ by Ref.  and by Ref.  for Maclaurin spheroids, to higher order in $`\mathrm{\Omega }`$ by Refs.  and , to all orders in $`\mathrm{\Omega }`$ by Ref. , and in relativistic stars by Ref. . The hybrid modes consist of mixtures of the $`r`$-modes of nonzero-buoyancy stars together with the zero-buoyancy limit of $`g`$-modes ; see also Ref.  for hints in this direction.

In this section we will re-derive the equations obtained by Ref.  that describe the leading order mode functions of the hybrid modes \[Eqs. (108)–(109) below\], by using a formal perturbation theory expansion of the quadratic eigenvalue equation (19). The main benefit of the analysis is that it shows that all of the rotational modes are orthogonal with respect to the inner product of the non-rotating star. In addition, it clarifies why, in the analysis of Ref. , it is sufficient to consider the curl of the perturbed Euler equation and to neglect its longitudinal part.

We also note that Ref.  shows that if a background star has Jordan chains of length $`p`$, and is subject to a perturbation of order $`\epsilon `$, then generically the leading order change in mode frequencies scales as $`\epsilon ^{1/(p+1)}`$. Now, non-rotating stars always have Jordan-chain modes of length 1 (see Appendix D), and the leading order effects of rotation on the mode dynamics is linear in the star’s angular velocity $`\mathrm{\Omega }`$. Hence, one might suspect the existence of modes in rotating stars whose frequencies scale as $`\sqrt{\mathrm{\Omega }}`$ as $`\mathrm{\Omega }0`$ \***\***\***Modes whose frequencies $`\omega `$ are proportional to $`\sqrt{\mathrm{\Omega }}`$ at large $`\mathrm{\Omega }`$ have been found by Ref. , but for these modes $`\omega \omega _00`$ as $`\mathrm{\Omega }0`$.. We show explicitly in Appendix F that no such modes exist, and explain why the argument of Ref.  does not apply in Appendix D.

### B Normalization convention

The formalism of Sec. II is invariant under changes of normalization of the basis functions $`𝝃_A`$; changes of the normalization of the modes are compensated for by changes of the normalization of the mode coefficients $`c_A(t)`$. In general one can choose any convenient normalization convention. However, in the context of the slow rotation expansion, statements about the behavior of the mode functions $`𝝃_A`$ as $`\mathrm{\Omega }0`$ are dependent on the normalization convention\*††\*††\*††For example, if one chooses to normalize an $`r`$ mode in such a way that an amplitude $`c_A`$ of order unity corresponds to a mode energy of order the stellar binding energy, then $`𝝃_A1/\mathrm{\Omega }`$ as $`\mathrm{\Omega }0`$.. Throughout the rest of this paper we shall adopt the convention

$$𝝃_A,𝝃_A=MR^2,$$
(78)
where $`M`$ is the stellar mass and $`R`$ the stellar radius. Then the mode coefficients are dimensionless, and order unity amplitudes $`c_A1`$ correspond to Lagrangian displacements of order size of the star. With this convention, all of the mode functions $`𝝃_A`$ will have finite limits as $`\mathrm{\Omega }0`$.

### C Perturbation expansion

We denote by $`_0`$ the subspace of the $``$ of Lagrangian displacements $`𝝃(𝐱)`$ spanned by modes of the non-rotating star of zero frequency. In the zero-buoyancy case, the space $`_0`$ consists of those perturbations $`𝝃(𝐱)`$ for which $`\delta \rho =(\rho 𝝃)`$ vanishes, and in the nonzero-buoyancy case it consists of purely axial vectors (under reasonable assumptions). These characterizations of the space of zero-frequency modes have been proved for fully relativistic stars by Ref. . In Appendix E we give Newtonian versions of the proofs of Ref. , and also derive some other properties of the space $`_0`$.

The operators $`𝐁`$ and $`𝐂`$ can be expanded as power series in the angular velocity $`\mathrm{\Omega }`$ as \*‡‡\*‡‡\*‡‡More precisely, the expansion parameter is the dimensionless quantity $`\mathrm{\Omega }\sqrt{R^3/(GM)}`$, where $`M`$ and $`R`$ are the stellar mass and radius.

$$𝐁=\mathrm{\Omega }𝐁^{(1)}$$
(79)
and

$$𝐂=𝐂^{(0)}+\mathrm{\Omega }^2𝐂^{(2)}+O(\mathrm{\Omega }^4).$$
(80)
Also the inner product (21) depends on $`\mathrm{\Omega }`$ through the background density $`\rho `$ and can be expanded as

$$𝝃,𝝃^{}=𝝃,𝝃^{}_0+\mathrm{\Omega }^2𝝃,𝝃^{}_2+O(\mathrm{\Omega }^4).$$
(81)
Since the operators $`𝐂`$ and $`i𝐁`$ are Hermitian with respect to $`,`$, it follows that $`𝐂^{(0)}`$ and $`i𝐁^{(1)}`$ (although not $`𝐂^{(2)}`$) are Hermitian with respect to $`,_0`$.

The operator $`𝐂^{(0)}`$ governs the modes of the non-rotating star. Let $`\omega _A^{(0)}`$, $`\widehat{𝝃}_A`$, for $`A=1,2,3\mathrm{}`$, be a complete set of solutions of the eigenvalue equation for the spherical star:

$`𝐂^{(0)}\widehat{𝝃}_A=\omega _A^{(0)\mathrm{\hspace{0.17em}2}}\widehat{𝝃}_A.`$ (82)
As in Sec. II C 4 above we can take $`A=(\alpha ,ϵ)`$ with $`ϵ=+`$ or $`ϵ=`$, and with $`\widehat{𝝃}_{\alpha ,+}=\widehat{𝝃}_{\alpha ,}\widehat{𝝃}_\alpha `$, and $`\omega _{\alpha ,+}^{(0)}=\omega _{\alpha ,}^{(0)}\omega _\alpha ^{(0)}0`$. Since $`𝐂^{(0)}`$ is a Hermitian operator, the basis $`\{\widehat{𝝃}_\alpha \}`$ is a complete, orthonormal basis of $``$. The basis $`\widehat{𝝃}_\alpha `$ can be chosen arbitrarily within each degenerate subspace of the operator $`𝐂^{(0)}`$. For example, the basis can be chosen arbitrarily within the subspace $`_0`$. Below we will introduce a different, preferred basis $`𝝃_\alpha ^{(0)}`$ of eigenvectors of $`𝐂^{(0)}`$ that are the $`\mathrm{\Omega }0`$ limit of modes of the rotating star.

Consider now trying to solve for the one parameter family or families of modes $`𝝃_\alpha (\mathrm{\Omega }),\omega _\alpha (\mathrm{\Omega })`$ of the rotating star for which $`\omega _\alpha (\mathrm{\Omega })\omega _\alpha ^{(0)}`$ as $`\mathrm{\Omega }0`$ †\*†\*†\*It suffices to consider the case $`ϵ=+`$ and to seek a one-parameter family of modes whose limiting frequency is $`\omega _{\alpha ,+}^{(0)}`$, since a one parameter family whose limiting frequency is $`\omega _{\alpha ,}^{(0)}`$ can then be obtained using the transformation $`(𝝃,\omega )(𝝃^{},\omega )`$; see Sec. II C 4 above.. As explained above, there are two types of such modes, rotational modes for which $`\omega _\alpha ^{(0)}=0`$, and non-rotational modes for which $`\omega _\alpha ^{(0)}0`$. For both types of modes, we make the ansatz that the frequency and mode function can be expanded as

$`\omega _\alpha (\mathrm{\Omega })=\omega _\alpha ^{(0)}+\mathrm{\Omega }\omega _\alpha ^{(1)}+\mathrm{\Omega }^2\omega _\alpha ^{(2)}+O(\mathrm{\Omega }^3),`$ (83)
and

$`𝝃_\alpha (\mathrm{\Omega })=𝝃_\alpha ^{(0)}+\mathrm{\Omega }𝝃_\alpha ^{(1)}+\mathrm{\Omega }^2𝝃_\alpha ^{(2)}+O(\mathrm{\Omega }^3).`$ (84)
If we substitute the expansions (83)–(84) and also the expansions (79)–(80) for the operators $`𝐁`$ and $`𝐂`$ into the quadratic eigenvalue equation (19) we get a series of equations. First, at order $`O(\mathrm{\Omega }^0)`$ we get

$`\left[\omega _\alpha ^{(0)\mathrm{\hspace{0.17em}2}}+𝐂^{(0)}\right]𝝃_\alpha ^{(0)}=0,`$ (85)
which is the usual eigenvalue equation for spherical stars. At order $`O(\mathrm{\Omega })`$ we get

$`\left[\omega _\alpha ^{(0)\mathrm{\hspace{0.17em}2}}+𝐂^{(0)}\right]𝝃_\alpha ^{(1)}\omega _\alpha ^{(0)}\left[2\omega _\alpha ^{(1)}+i𝐁^{(1)}\right]𝝃_\alpha ^{(0)}=0.`$ (86)
Finally, at order $`O(\mathrm{\Omega }^2)`$, we get the equation

$`\left[\omega _\alpha ^{(0)\mathrm{\hspace{0.17em}2}}+𝐂^{(0)}\right]𝝃_\alpha ^{(2)}\omega _\alpha ^{(0)}\left[2\omega _\alpha ^{(1)}+i𝐁^{(1)}\right]𝝃_\alpha ^{(1)}`$ (87)
$`+\left[2\omega _\alpha ^{(0)}\omega _\alpha ^{(2)}\omega _\alpha ^{(1)\mathrm{\hspace{0.17em}2}}i\omega _\alpha ^{(1)}𝐁^{(1)}+𝐂^{(2)}\right]𝝃_\alpha ^{(0)}=0.`$ (88)
We next express these perturbation equations in the basis of modes of the non-rotating star. We expand the mode functions $`𝝃_\alpha ^{(0)}`$, $`𝝃_\alpha ^{(1)}`$ and $`𝝃_\alpha ^{(2)}`$ in terms of the eigenbasis $`\{\widehat{𝝃}_\alpha \}`$ of the non-rotating star as

$`𝝃_\alpha ^{(0)}={\displaystyle \underset{\beta }{}}c_{\alpha \beta }^{(0)}\widehat{𝝃}_\beta ,`$ (89)
$`𝝃_\alpha ^{(1)}={\displaystyle \underset{\beta }{}}c_{\alpha \beta }^{(1)}\widehat{𝝃}_\beta ,`$ (90)
and

$`𝝃_\alpha ^{(2)}={\displaystyle \underset{\beta }{}}c_{\alpha \beta }^{(2)}\widehat{𝝃}_\beta .`$ (91)
Note that we cannot assume that $`c_{\alpha \beta }^{(0)}=\delta _{\alpha \beta }`$ because of the possibility of degeneracies. Just as in degenerate perturbation theory in quantum mechanics, the modes of the rotating star define a preferred basis $`𝝃_\alpha ^{(0)}`$ of each degenerate subspace. Since we do not know this basis before we solve for the modes of the rotating star, we must start our computation with some arbitrary choice of basis $`\widehat{𝝃}_\alpha `$ in each degenerate subspace.

Next, we substitute the expansion (89) of $`𝝃_\alpha ^{(0)}`$ into the zeroth-order eigenvalue equation (85). This gives, using Eq. (82),

$`{\displaystyle \underset{\beta }{}}c_{\alpha \beta }^{(0)}\left[\omega _\alpha ^{(0)\mathrm{\hspace{0.17em}2}}\omega _\beta ^{(0)\mathrm{\hspace{0.17em}2}}\right]\widehat{𝝃}_\beta =0.`$ (92)
It follows that $`c_{\alpha \beta }^{(0)}=0`$ when $`\omega _\alpha ^{(0)}\omega _\beta ^{(0)}`$. In other words, the matrix $`c_{\alpha \beta }^{(0)}=0`$ is block diagonal, with each block corresponding to a degenerate subspace. Hence the zeroth order mode function can be written as, from Eq. (89),

$`𝝃_\alpha ^{(0)}={\displaystyle \underset{\beta \mathrm{with}\omega _\beta ^{(0)}=\omega _\alpha ^{(0)}}{}}c_{\alpha \beta }^{(0)}\widehat{𝝃}_\beta ,`$ (93)
and thus lies within the degenerate subspace corresponding to the eigenvalue $`\omega _\alpha ^{(0)}`$. In particular, for the modes with $`\omega _\alpha ^{(0)}=0`$, the zeroth order mode function $`𝝃_\alpha ^{(0)}`$ lies inside the space $`_0`$.

The analysis now divides into two cases for the two different types of mode discussed above.

### D Modes with finite frequency in the non-rotating limit

If we substitute the expansions (89) and (90) for $`𝝃_\alpha ^{(0)}`$ and $`𝝃_\alpha ^{(1)}`$ into the first-order eigenvalue equation (86), and then take the zeroth order inner product $`,_0`$ with $`\widehat{𝝃}_\beta `$, we get

$`\left[(\omega _\alpha ^{(0)})^2(\omega _\beta ^{(0)})^2\right]c_{\alpha \beta }^{(1)}`$ $`+2\omega _\alpha ^{(0)}\omega _\alpha ^{(1)}c_{\alpha \beta }^{(0)}`$ (95)
$`+i\omega _\alpha ^{(0)}{\displaystyle \underset{\gamma }{}}B_{\beta \gamma }^{(1)}c_{\alpha \gamma }^{(0)}=0,`$
where

$`B_{\beta \gamma }^{(1)}\widehat{𝝃}_\beta ,𝐁^{(1)}\widehat{𝝃}_\gamma _0.`$ (96)
If we specialize Eq. (95) to values of $`\beta `$ for which $`\omega _\beta ^{(0)}=\omega _\alpha ^{(0)}`$, and divide across by $`\omega _\alpha ^{(0)}`$ we get

$`2\omega _\alpha ^{(1)}c_{\alpha \beta }^{(0)}+i{\displaystyle \underset{\gamma }{}}B_{\beta \gamma }^{(1)}c_{\alpha \gamma }^{(0)}=0,`$ (97)
which is a standard eigenvalue equation. Thus, the zeroth order modes are given by diagonalizing the operator $`B_{\beta \gamma }^{(1)}`$ within the degenerate subspace of modes $`\beta `$ with $`\omega _\beta ^{(0)}=\omega _\alpha ^{(0)}`$. The eigenvalues of $`B_{\beta \gamma }^{(1)}`$ give the first order changes $`\omega _\alpha ^{(1)}`$ to the eigenfrequency, just as in degenerate perturbation theory in quantum mechanics.

Suppose now we switch to the basis that diagonalizes $`B_{\alpha \beta }^{(1)}`$ within each degenerate subspace, so that $`c_{\alpha \beta }^{(0)}=\delta _{\alpha \beta }`$. Then, applying Eq. (95) for values of $`\beta `$ with $`\omega _\beta ^{(0)}\omega _\alpha ^{(0)}`$ gives for the first order change in the mode functions

$`c_{\alpha \beta }^{(1)}={\displaystyle \frac{i\omega _\beta ^{(0)}}{\omega _\alpha ^{(0)\mathrm{\hspace{0.17em}2}}\omega _\beta ^{(0)\mathrm{\hspace{0.17em}2}}}}B_{\beta \alpha }^{(1)}.`$ (98)
Using the expansions (84), (89), and (90) together with Eq. (98) shows that the inner product of two mode functions with $`\omega _\alpha ^{(0)}\omega _\beta ^{(0)}`$ is

$$𝝃_\alpha (\mathrm{\Omega }),𝝃_\beta (\mathrm{\Omega })=\frac{i\mathrm{\Omega }}{\omega _\alpha ^{(0)}+\omega _\beta ^{(0)}}B_{\alpha \beta }^{(1)}+O(\mathrm{\Omega }^2).$$
(99)
This explicitly demonstrates that modes in a rotating star are not orthogonal in general, since we can find pairs of modes $`\alpha `$, $`\beta `$ of the non-rotating star for which the matrix element $`B_{\alpha \beta }^{(1)}`$ is non vanishing.

### E Modes with vanishing frequency in the non-rotating limit

Things work somewhat differently for the second class of modes. If we substitute $`\omega _\alpha ^{(0)}=0`$ into Eq. (86) we get

$`𝐂^{(0)}𝝃_\alpha ^{(1)}=0.`$ (100)
Thus, the first order equation (86) does not determine $`𝝃_\alpha ^{(1)}`$, except to dictate that $`𝝃_\alpha ^{(1)}`$ lie in the space $`_0`$ of zero-frequency modes of the non-rotating star. In addition, the ambiguity in the choice of basis of $`𝝃_\alpha ^{(0)}`$ of $`_0`$ is not resolved at this linear order.

Hence we must use the second order equation (88). For $`\omega _\alpha ^{(0)}=0`$, this equation reduces to

$`𝐂^{(0)}𝝃_\alpha ^{(2)}+\left[\omega _\alpha ^{(1)\mathrm{\hspace{0.17em}2}}i\omega _\alpha ^{(1)}𝐁^{(1)}+𝐂^{(2)}\right]𝝃_\alpha ^{(0)}=0.`$ (101)
If we substitute the expansions (89) and (91) into Eq. (101) and take the zeroth order inner product $`,_0`$ with $`\widehat{𝝃}_\beta `$, we get

$`\omega _\alpha ^{(1)\mathrm{\hspace{0.17em}2}}c_{\alpha \beta }^{(0)}+\omega _\beta ^{(0)\mathrm{\hspace{0.17em}2}}c_{\alpha \beta }^{(2)}+{\displaystyle \underset{\gamma }{}}\left[i\omega _\alpha ^{(1)}B_{\beta \gamma }^{(1)}+C_{\beta \gamma }^{(2)}\right]c_{\alpha \gamma }^{(0)}=0,`$ (102)
where

$`C_{\beta \gamma }^{(2)}\widehat{𝝃}_\beta ,𝐂^{(2)}\widehat{𝝃}_\gamma _0.`$ (103)
We now specialize Eq. (102) to values of $`\beta `$ for which $`\omega _\beta ^{(0)}=0`$, i.e., project both sides of the equation into $`_0`$ ††††††The component of Eq. (102) orthogonal to $`_0`$ is not needed to determine the leading order quantities $`𝝃_\alpha ^{(0)}`$ and $`\omega _\alpha ^{(1)}`$. It determines the higher order correction $`𝝃_\alpha ^{(2)}`$ to the mode function.. This gives

$`\omega _\alpha ^{(1)\mathrm{\hspace{0.17em}2}}c_{\alpha \beta }^{(0)}+{\displaystyle \underset{\gamma }{}}\left[i\omega _\alpha ^{(1)}B_{\beta \gamma }^{(1)}+C_{\beta \gamma }^{(2)}\right]c_{\alpha \gamma }^{(0)}=0,`$ (104)
which is a quadratic eigenvalue equation that determines both $`\omega _\alpha ^{(1)}`$ and $`c_{\alpha \beta }^{(0)}`$. Equation (104) can be rewritten as

$`𝐏_0\left[\omega _\alpha ^{(1)\mathrm{\hspace{0.17em}2}}i\omega _\alpha ^{(1)}𝐁^{(1)}+𝐂^{(2)}\right]𝐏_0𝝃_\alpha ^{(0)}=0,`$ (105)
where $`𝐏_0`$ is the orthogonal projection operator (with respect to $`,_0`$) that projects into $`_0`$ (see appendix E).

It is at this point in the analysis that the difference between stars with zero and non-zero buoyancy enters. We treat the two cases separately.

#### 1 Zero-buoyancy stars

When $`𝐀=0`$, the term $`𝐏_0𝐂^{(2)}𝐏_0𝝃_\alpha ^{(0)}`$ in Eq. (105) vanishes. To see this, let

$`𝜿=𝐂𝝃_\alpha ^{(0)}=𝐂𝐏_0𝝃_\alpha ^{(0)}.`$ (106)
It follows from the formula (14) for the operator $`𝐂`$ that $`𝜿(𝐱)`$ must be a pure gradient, so that $`\times 𝜿=0`$. However, Eq. (E6) of Appendix E then implies that $`𝐏_0𝜿=0`$, and expanding this last equation to second order in $`\mathrm{\Omega }`$ gives the desired result.

Equation (105) now simplifies to

$`\left[𝐏_0i𝐁^{(1)}𝐏_0\right]𝝃_\alpha ^{(0)}=\omega _\alpha ^{(1)}𝝃_\alpha ^{(0)},`$ (107)
which is a standard eigenvalue equation in $`_0`$. Thus, the zeroth order mode functions $`𝝃_\alpha ^{(0)}`$ for the hybrid modes are simply the eigenvectors in $`_0`$ of the Hermitian operator $`i𝐏_0𝐁^{(1)}𝐏_0`$. It follows that these zeroth order mode functions are orthogonal to each other with respect the inner product $`,_0`$ of the non-rotating star. They are also orthogonal to the zeroth order mode functions of the modes with $`\omega _\alpha ^{(0)}0`$, since those mode functions $`𝝃_\alpha ^{(0)}`$ lie completely inside degenerate subspaces that are orthogonal to $`_0`$, by Eq. (93).

We now rewrite the equations defining the hybrid modes $`𝝃_\alpha ^{(0)}(𝐱)`$ in a more accessible notation. First, from Eq. (93) we know that $`𝝃_\alpha ^{(0)}`$ lies in $`_0`$, which from Appendix E is equivalent to

$`(\rho 𝝃_\alpha ^{(0)})=0.`$ (108)
Second, Eq. (107) can be rewritten using the definition (8) of $`𝐁`$ and Eq. (E6) as

$`\times \left[i\omega _\alpha ^{(1)}𝝃_\alpha ^{(0)}+2𝛀\times 𝝃_\alpha ^{(0)}\right]=0.`$ (109)
Equations (108) and (109) are the equations that are solved to obtain the hybrid modes in Ref. .

Finally we note that going from Eq. (105) to Eq. (107) entailed dividing by $`\omega _\alpha ^{(1)}`$. There is another solution to Eq. (105) with $`\omega _\alpha ^{(1)}=0`$. Examining the perturbation expansion to higher order for this solution shows that $`\omega _\alpha ^{(2)}=0`$ also. We suspect that this solution is a zero-frequency solution to all orders in $`\mathrm{\Omega }`$. These modes are not Jordan-chain modes, but they are pure gauge in the sense of Ref. , at least to linear order in the mode amplitude. We note that it is impractical to include these modes (and also the similar modes of Appendix D) in nonlinear mode evolution calculations since the perturbation expansion in mode amplitude for these modes typically breaks down on a timescale $`1/\mathrm{\Omega }`$.

#### 2 Stars with buoyancy

For $`𝐀0`$ (stars for which there is a nonzero buoyancy force), the space $`_0`$ is the space of purely axial vectors (see Appendix E). In this case the quadratic eigenvalue equation (105) for the mode functions can be written as

$$𝐏_0\left[\omega _\alpha ^{(1)\mathrm{\hspace{0.17em}2}}𝝃_\alpha ^{(0)}2i\omega _\alpha ^{(1)}𝛀\times 𝝃_\alpha ^{(0)}+𝐂_a^{(2)}𝝃_\alpha ^{(0)}+𝐂_b^{(2)}𝝃_\alpha ^{(0)}\right]=0,$$
(110)
where $`𝐏_0𝐯`$ is now the axial part of $`𝐯`$, and $`𝐂_a^{(2)}`$ and $`𝐂_b^{(2)}`$ are the second order in $`\mathrm{\Omega }`$ pieces of the operators (14) and (16). The term $`𝐂_a^{(2)}𝝃_\alpha ^{(0)}`$ drops out since it is a gradient and thus has no axial part, and the term $`𝐂_b^{(2)}𝝃_\alpha ^{(0)}`$ drops out since it depends only on $`𝝃_\alpha ^{(0)}`$ and $`𝝃_\alpha ^{(0)}𝐞_r`$ which vanishes as $`𝝃_\alpha ^{(0)}`$ is purely axial. Using Eq. (E11), the equation now reduces to

$`𝐫\times \left[i\omega _\alpha ^{(1)}𝝃_\alpha ^{(0)}+2𝛀\times 𝝃_\alpha ^{(0)}\right]=0.`$ (111)
In Ref.  it is shown that Eq. (111) determines the frequency $`\omega _\alpha ^{(1)}`$ but not the eigenfunction $`𝝃_\alpha ^{(0)}`$, and that one must go to higher order in the $`\mathrm{\Omega }`$ expansion in order to determine the eigenfunction; see also Refs.  and .

## IV Second Order Lagrangian Perturbation Theory in Rotating Stars

In this section we develop a second order Lagrangian perturbation theory for rotating stars, extending previous analyses of non-rotating stars (see, e.g., Refs. ). The main new result is that the expression for the three-mode coupling coefficient in terms of the mode functions \[Eq. (134) below\] in the same for rotating stars as for non-rotating stars. Although we restrict attention in this paper to second order perturbation theory, we lay the foundations in Appendix J for computing mode coupling coefficients at third and higher orders.

### A Second order equation of motion

To describe the fluid perturbation we use the rotating-frame Lagrangian displacement $`𝝃(𝐱,t)`$. The definition of this quantity is the same in a nonlinear context as in linear theory: a fluid element at rotating-frame location $`𝐱`$ in the unperturbed star is moved to rotating-frame location $`𝐱+𝝃(𝐱,t)`$ by the perturbation. We assume that the pressure $`p`$ in the background star is a function $`p(\rho ,\stackrel{}{\mu })`$ of the the background mass density $`\rho `$ and of a vector $`\stackrel{}{\mu }`$ of other fluid variables such as entropy or composition, which can vary with position. We also assume that the fluid variables $`\stackrel{}{\mu }`$ of each fluid element are preserved by the perturbation.

Under these assumptions, the equation of motion for $`𝝃`$ is derived in Appendix I to second order in $`𝝃`$ via direct perturbation of the Euler equation, and in Appendix J to any order in $`𝝃`$ from a variational principle. The result has the form

$`\ddot{𝝃}+𝐁\dot{𝝃}+𝐂𝝃=𝐚[𝝃],`$ (112)
where the left hand side contains all the linear terms \[cf. Eq. (7) above\] and the right hand side $`𝐚[𝝃]`$ is a nonlinear acceleration. That acceleration can be written as $`𝐚=𝐚_P+𝐚_G`$, where $`𝐚_P`$ and $`𝐚_G`$ are the contributions from pressure gradients and gravity; the corresponding force densities are $`𝐟=\rho 𝐚`$, $`𝐟_P=\rho 𝐚_P`$, and $`𝐟_G=\rho 𝐚_G`$.

We expand the acceleration as

$$𝐚[𝝃]=𝐚^{(2)}[𝝃,𝝃]+O(𝝃^3)=𝐚_P^{(2)}[𝝃,𝝃]+𝐚_G^{(2)}[𝝃,𝝃]+O(𝝃^3),$$
(113)
where $`𝐚^{(2)}[𝝃,𝝃^{}]`$, $`𝐚_P^{(2)}[𝝃,𝝃^{}]`$, and $`𝐚_G^{(2)}[𝝃,𝝃^{}]`$ are symmetric bilinear functions of their arguments. In order to give the explicit results for the second order pieces $`𝐚_P^{(2)}[𝝃,𝝃]`$ and $`𝐚_G^{(2)}[𝝃,𝝃]`$ of the pressure gradient and gravitational accelerations, we need to introduce some notations. The first and second order Eulerian perturbations $`\delta ^{(1)}\varphi `$ and $`\delta ^{(2)}\varphi `$ to the Newtonian potential are defined as functionals of $`𝝃`$ by the equations \[see Eq. (I37) below\]

$$\frac{1}{4\pi G}^2\delta ^{(1)}\varphi =(\rho 𝝃)$$
(114)
and

$$\frac{1}{4\pi G}^2\delta ^{(2)}\varphi =\frac{1}{2}_i_j(\rho \xi ^i\xi ^j).$$
(115)
We define the tensors $`\mathrm{\Theta }_j^i`$, $`\mathrm{\Xi }_j^i`$ and $`\chi _j^i`$ by

$`\mathrm{\Theta }_j^i`$ $``$ $`\xi _{;j}^i\xi _{;k}^k,`$ (116)
$`\mathrm{\Xi }_j^i`$ $``$ $`\xi _{;k}^i\xi _{;j}^k,`$ (117)
and

$$\chi _j^i\xi _{;l}^i\xi _{;k}^l\xi _{;j}^k.$$
(118)
Here semicolons denote covariant derivatives, for example $`\xi _{;j}^i_j\xi ^i`$. Finally, we denote by $`\mathrm{\Gamma }_1`$ the generalized adiabatic index governing the perturbations, defined by Eq. (I20) or Eq. (J27).

The second order gravitational acceleration is \[Eqs. (I59) and (J35) below\]

$$a_{Gj}^{(2)}[𝝃,𝝃]=_j\delta ^{(2)}\varphi \xi ^k_k_j\delta ^{(1)}\varphi \frac{1}{2}\xi ^k\xi ^l_k_l_j\varphi .$$
(119)
The second order pressure term is \[Eqs. (I59), (J26) and (J31) below\]

$$a_{Pi}^{(2)}[𝝃,𝝃]=\frac{1}{\rho }_j\left[p(\mathrm{\Gamma }_11)\mathrm{\Theta }_i^j+p\mathrm{\Xi }_i^j+\mathrm{\Psi }\delta _i^j\right],$$
(120)
where

$$\mathrm{\Psi }=\frac{1}{2}p\mathrm{\Theta }\left[(\mathrm{\Gamma }_11)^2+\frac{\mathrm{\Gamma }_1}{\mathrm{ln}\rho }\right]+\frac{1}{2}p(\mathrm{\Gamma }_11)\mathrm{\Xi },$$
(121)
and $`\mathrm{\Theta }`$ and $`\mathrm{\Xi }`$ are the traces of $`\mathrm{\Theta }_i^j`$ and $`\mathrm{\Xi }_i^j`$. More general expressions that are valid to all orders in $`𝝃`$ are given in Eqs. (J32) below.

A key feature of the equation of motion (112) is that there is no dependence on the angular velocity $`\mathrm{\Omega }`$ of the rotating frame in the nonlinear terms on the right hand side. The dependence on $`\mathrm{\Omega }`$ occurs only in the linear terms on the left hand side. This can be understood most easily from the variational principle that underlies the hydrodynamical equations. As explained in Appendix J, the only term in the action that is not explicitly invariant under a change of rotational frame is the kinetic energy term, and that kinetic energy term is quadratic in the Lagrangian displacement $`𝝃`$ and its time derivative $`\dot{𝝃}`$.

Finally, we note that in Lagrangian perturbation theory only pressure and gravity can contribute to the nonlinear force. In spite of the fact that the nonlinear advection term $`𝐯\mathbf{}𝐯`$ appears in the Euler equation, no nonlinear advection terms are present in the Lagrangian framework, since the kinetic energy is explicitly quadratic in $`𝝃`$.

### B Coupled equations of motion for the mode coefficients

We now use the equation of motion given by Eqs. (112), (113), (119) and (120) to derive the coupled, nonlinear ordinary differential equations satisfied by the mode expansion coefficients $`c_A(t)`$. The basic equation (7) of Sec. II can be applied with the externally applied acceleration $`𝐚_{\mathrm{ext}}(𝐱,t)`$ replaced by the nonlinear acceleration $`𝐚[𝝃]`$. Therefore we can directly carry over the equation of motion (46), which gives

$$\dot{c}_A(t)+i\omega _Ac_A(t)=\frac{i}{b_A}𝝃_A,𝐚^{(2)}[𝝃,𝝃]+O(𝝃^3).$$
(122)
Next, we take the complex conjugate of the phase space mode expansion (43) and use the fact that $`𝝃(𝐱,t)`$ is real to obtain

$$𝝃(𝐱,t)=\underset{A}{}c_A(t)^{}𝝃_A(𝐱)^{}.$$
(123)
Substituting this expansion into the right hand side of Eq. (122) gives

$$\dot{c}_A(t)+i\omega _Ac_A(t)=\frac{i}{b_A}\underset{B,C}{}\kappa _{ABC}^{}c_B(t)^{}c_C(t)^{}+O(c^3).$$
(124)
where the three-mode coupling coefficient is

$$\kappa _{ABC}=𝝃_A^{},𝐚^{(2)}[𝝃_B,𝝃_C]$$
(125)
which is completely symmetric in the indices $`A`$, $`B`$ and $`C`$. For $`\omega _A0`$, this equation can also be written as

$$\dot{c}_A(t)+i\omega _Ac_A(t)=i\omega _A\underset{B,C}{}\frac{\kappa _{ABC}^{}}{\epsilon _A}c_B(t)^{}c_C(t)^{}+O(c^3),$$
(126)
where we have used Eq. (K25). This form of the equation has a simple physical significance, since $`\kappa _{ABC}/\epsilon _A`$ is the ratio of the nonlinear interaction energy at unit amplitude to the energy of the mode at unit amplitude.

The equation of motion (124) is valid for arbitrary complex $`𝝃(𝐱,t)`$. However, for the physically relevant case of real $`𝝃(𝐱,t)`$, the coefficients $`c_A(t)`$ will not all be independent, as explained in Sec. II C 4 above. The modes will occur in complex conjugate pairs $`(𝝃,\omega )`$ and $`(𝝃^{},\omega )`$, and it is convenient to use the set of modes $`𝝃_\alpha `$ defined in Sec. II C 4 consisting of one mode from each complex-conjugate pair. For this set of modes, the mode expansion is given by Eq. (74):

$$𝝃(t)=\underset{\alpha }{}\left[c_\alpha (t)𝝃_\alpha +c_\alpha (t)^{}𝝃_\alpha ^{}\right],$$
(127)
and the inverse mode expansion by Eq. (75). For this expansion, the variables $`c_\alpha (t)`$ are all independent. If we now use the expansion (127) in the equation of motion (76) we find as a replacement for Eq. (124) the equation

$$\dot{c}_\alpha +i\omega _\alpha c_\alpha =\frac{i}{b_\alpha }\underset{\beta ,\gamma }{}\left[\kappa _{\overline{\alpha }\beta \gamma }c_\beta c_\gamma +\kappa _{\overline{\alpha }\overline{\beta }\gamma }c_\beta ^{}c_\gamma +\kappa _{\overline{\alpha }\beta \overline{\gamma }}c_\beta c_\gamma ^{}+\kappa _{\overline{\alpha }\overline{\beta }\overline{\gamma }}c_\beta ^{}c_\gamma ^{}\right],$$
(128)
where

$$\kappa _{\alpha \beta \gamma }=𝝃_\alpha ^{},𝐚^{(2)}[𝝃_\beta ,𝝃_\gamma ]$$
(129)
which is symmetric in $`\alpha `$, $`\beta `$, and $`\gamma `$. In Eq. (128) we have used a notational convention for the coupling coefficients $`\kappa `$ where a bar over an index means that the corresponding mode function is to be complex conjugated in the expression for the coupling coefficient. For example, we have

$$\kappa _{\alpha \beta \overline{\gamma }}=𝝃_\alpha ^{},𝐚^{(2)}[𝝃_\beta ,𝝃_\gamma ^{}].$$
(130)
Equation (128) is the final equation of motion. To model the saturation of $`r`$-modes, it will have to be supplemented by viscous damping and gravitational radiation reaction terms. We now turn to evaluating the coupling coefficients $`\kappa _{ABC}`$ to insert into this equation, with the understanding that each capital Roman index can be either a unbarred or a barred Greek lower case index, for example $`A=\alpha `$ or $`A=\overline{\alpha }`$. In the notation of Sec. II C 4, $`(\alpha ,+)`$ corresponds to $`\alpha `$ and $`(\alpha ,)`$ corresponds to $`\overline{\alpha }`$.

### C Explicit expression for the three-mode coupling coefficient

We can obtain an explicit expression for the coupling coefficient $`\kappa _{ABC}`$ in terms of the mode functions $`𝝃_A(𝐱)`$, $`𝝃_B(𝐱)`$ and $`𝝃_C(𝐱)`$ by substituting the formulae (113), (119) and (120) for the nonlinear acceleration into the definition (125) of $`\kappa _{ABC}`$, using the definition (21) of the inner product, and performing several integrations by parts. The computation can be simplified by (i) computing the functional $`𝝃^{},𝐚^{(2)}[𝝃,𝝃]`$ of $`𝝃`$, (ii) equating this functional to $`\kappa (𝝃,𝝃,𝝃)`$ where $`\kappa (𝝃,𝝃^{},𝝃^{\prime \prime })`$ is a symmetric trilinear functional of its arguments, in order to determine the functional $`\kappa (𝝃,𝝃^{},𝝃^{\prime \prime })`$, and (iii) evaluating $`\kappa (𝝃_A,𝝃_B,𝝃_C)`$ to obtain the coupling coefficient $`\kappa _{ABC}`$. For the gravitational terms, the domain of the spatial integral can be taken to be all of space, if one includes the $`\delta `$-function contributions in Eq. (I37). The gravitational terms can then be integrated by parts, and the boundary terms at infinity which are generated vanish. The domain of integration for the pressure terms can be taken to be the interior of the star; the boundary terms generated by the integration by parts vanish since $`p=0`$ on the stellar surface. The final result is

$`\kappa _{ABC}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d^3xp[(\mathrm{\Gamma }_11)(\mathrm{\Xi }_{AB}\mathbf{}𝝃_C+\mathrm{\Xi }_{BC}\mathbf{}𝝃_A+\mathrm{\Xi }_{CA}\mathbf{}𝝃_B)`$ (131)
$`+`$ $`\{(\mathrm{\Gamma }_11)^2+{\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{ln}\rho }}\}\mathbf{}𝝃_A\mathbf{}𝝃_B\mathbf{}𝝃_C+\chi _{ABC}+\chi _{ACB}]`$ (132)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d^3x\rho [\xi _A^i\xi _B^j\delta ^{(1)}\varphi _{C;ij}+\xi _B^i\xi _C^j\delta ^{(1)}\varphi _{A;ij}+\xi _C^i\xi _A^j\delta ^{(1)}\varphi _{B;ij}`$ (133)
$`+`$ $`\xi _A^i\xi _B^j\xi _C^k\varphi _{;ijk}].`$ (134)
The notations here are as follows. The quantity $`\delta ^{(1)}\varphi _A`$ is defined by Eq. (114) with $`𝝃`$ on the right hand side replaced by $`𝝃_A`$. We define the functions

$`\mathrm{\Theta }_j^i[𝝃,𝝃^{}]`$ $``$ $`\xi _{;j}^i\xi _{;k}^k,`$ (135)
$`\mathrm{\Xi }_j^i[𝝃,𝝃^{}]`$ $``$ $`\xi _{;k}^i\xi _{;j}^k,`$ (136)
and

$$\chi _j^i[𝝃,𝝃^{},𝝃^{\prime \prime }]\xi _{;l}^i\xi _{;k}^l\xi _{;j}^{\prime \prime k}.$$
(137)
The tensors $`\mathrm{\Theta }_j^i`$, $`\mathrm{\Xi }_j^i`$ and $`\chi _j^i`$ which we defined previously in Eq. (117) are obtained from these expressions evaluated at $`𝝃=𝝃^{}=𝝃^{\prime \prime }`$. Finally we define the traces

$$\mathrm{\Xi }_{AB}=\delta _i^j\mathrm{\Xi }_j^i[𝝃_A,𝝃_B]$$
(138)
and

$$\chi _{ABC}=\delta _i^j\chi _j^i[𝝃_A,𝝃_B,𝝃_C].$$
(139)
The domain of integration for the gravitational term in Eq. (134) is all of space. The domain of integration for the pressure terms can be taken to be either the interior of the star, or all of space.

Note that the expression (134) for the three mode coupling coefficient is valid for arbitrary angular velocities $`\mathrm{\Omega }`$ of the background star. The expression has no explicit dependence on $`\mathrm{\Omega }`$, however, and therefore is also valid for non-rotating stars. In fact the expression (134) agrees with previous results  for non-rotating stars, except for an overall factor of $`3`$ due to differing conventions for the form of the equation of motion. Some of the previous analyses omitted the pressure force term containing the factor $`\mathrm{\Gamma }_1/\mathrm{ln}\rho `$, and some used the Cowling approximation and thus omitted the gravitational terms†‡†‡†‡Ref.  makes a slight error in calculating the gravitational terms. However, since they use the Cowling approximation in which $`\delta \varphi =0`$, this error does not affect their results..

For the special case of constant $`\mathrm{\Gamma }_1`$, it is possible to obtain an alternative form of the pressure terms in the coupling coefficient by (i) taking the domain of integration in both the gravitational and pressure terms to be all of space, and (ii) performing further integrations by parts . The result can be expressed as

$$\kappa _{ABC}=\kappa (𝝃_A,𝝃_B,𝝃_C)$$
(140)
where

$`\kappa (𝝃,𝝃,𝝃)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle }d^3x[\xi ^i\xi ^j\xi ^kp_{;ijk}+3(\mathbf{}𝝃)\xi ^i\xi ^jp_{;ij}\mathrm{\Gamma }_1(\mathrm{\Gamma }_1+1)p(\mathbf{}𝝃)^3`$ (141)
$`+`$ $`3\xi ^i\xi ^j(\delta p)_{;ij}+6(\mathbf{}𝝃)\xi ^i(\delta p)_{;i}+3\xi ^i\xi ^j\rho (\delta \varphi )_{;ij}+\xi ^i\xi ^j\xi ^k\rho \varphi _{;ijk}],`$ (142)
and $`\delta p=\mathrm{\Gamma }_1p\left(\mathbf{}𝝃\right)𝝃\mathbf{}p`$ is the Eulerian pressure perturbation. The domain of integration here is all of space, and therefore one should include the $`\delta `$-functions at the stellar surface in the factors $`p_{;ijk}`$, $`p_{;ij}`$ and possibly $`\varphi _{;ijk}`$. However, one can show that all of the contributions from the $`\delta `$ functions cancel in general, and therefore one can take the domain of integration in the expression (142) to be the interior of the star.

### D Application to Rossby Modes

Although the expression for the three-mode coupling coefficient is the same for rotating stars as for non-rotating stars, it is much more difficult to evaluate the coupling coefficient for rotating stars, for two reasons. First, including centrifugal flattening of the background pressure and density profiles makes the integrals much more difficult to perform. Second, the mode functions $`𝝃_A(𝐱)`$ for rotating stars typically consist of a sum of terms with different values of the quantum number $`l`$, which complicates the calculation.

We shall avoid the first of these difficulties by using the slow rotation expansion. First, we will compute the mode functions to leading order in $`\mathrm{\Omega }`$, as explained in Sec. III above. Second, in evaluating the coupling coefficient expression (134), we will use the background density, pressure and potential profiles of a spherical star, i.e., we will use the $`\mathrm{\Omega }0`$ limits of the variables of the background star. This procedure will give the zeroth order piece of the coupling coefficient $`\kappa _{ABC}`$ in an expansion in powers of $`\mathrm{\Omega }`$. In cases where this zeroth order piece vanishes, the coupling coefficient might be nonzero at higher orders in $`\mathrm{\Omega }`$.

Before discussing in more detail how to compute the coefficients, we describe some useful selection rules.

## V Selection Rules for the Coupling Coefficients

### A Terminology

In this section we derive a number of selection rules that apply to the three-mode coupling coefficient. We start by reviewing some terminology to describe modes. Any vector field $`𝝃(𝐱)`$ can be expanded in vector spherical harmonics as

$$𝝃(r,\theta ,\phi )=\underset{lm}{}\left[A_{lm}(r)Y_{lm}(\theta ,\phi )𝐞_{\widehat{r}}+B_{lm}(r)\mathbf{}Y_{lm}(\theta ,\phi )+C_{lm}(r)𝐫\times Y_{lm}(\theta ,\phi )\right],$$
(143)
where $`(r,\theta ,\phi )`$ are spherical polar coordinates and $`𝐞_{\widehat{r}}`$, $`𝐞_{\widehat{\theta }}`$ and $`𝐞_{\widehat{\phi }}`$ is the associated orthonormal basis. Modes for which $`A_{lm}=B_{lm}=0`$ for all $`l,m`$ are called axial modes. We define $`h`$ to be the parity map $`(x,y,z)(x,y,z)`$ and $`f`$ be the map $`(x,y,z)(x,y,z)`$ or $`(r,\theta ,\phi )(r,\pi \theta ,\phi )`$ which we call $`z`$-parity. Rotating stars, both Newtonian and relativistic, are invariant under these maps, and therefore we can always choose to use a basis of modes for which all the modes have definite parity and $`z`$-parity transformation properties; see Refs. . Specifically, we define the pullback $`f_{}𝝃`$ of $`𝝃`$ under $`f`$ by

$`f_{}𝝃`$ $`=`$ $`f_{}\left[\xi ^{\widehat{r}}(r,\theta ,\phi )𝒆_{\widehat{r}}+\xi ^{\widehat{\theta }}(r,\theta ,\phi )𝒆_{\widehat{\theta }}+\xi ^{\widehat{\phi }}(r,\theta ,\phi )𝒆_{\widehat{\phi }}\right]`$ (144)
$`=`$ $`\xi ^{\widehat{r}}(r,\pi \theta ,\phi )𝒆_{\widehat{r}}\xi ^{\widehat{\theta }}(r,\pi \theta ,\phi )𝒆_{\widehat{\theta }}+\xi ^{\widehat{\phi }}(r,\pi \theta ,\phi )𝒆_{\widehat{\phi }}.`$ (145)
Then we can always find a basis for which all modes are either z-parity even, satisfying

$$f_{}𝝃=𝝃,$$
(146)
or z-parity odd, satisfying

$$f_{}𝝃=𝝃.$$
(147)
Note that $`f_{}Y_{lm}=(1)^{l+m}Y_{lm}`$. A similar statement is true for the regular parity map $`h`$; we can take all modes to be either parity even or parity odd. Finally we can take all modes to have definite values of $`m`$, so that if $`z`$ is the mapping $`(r,\theta ,\phi )(r,\theta ,\phi +\mathrm{\Delta }\phi )`$, then

$$z_{}𝝃=\mathrm{exp}[im\mathrm{\Delta }\phi ]𝝃.$$
(148)
Recall also that all modes can be classified as either rotational ($`\omega 0`$ as $`\mathrm{\Omega }0`$) or regular ($`\omega `$ finite as $`\mathrm{\Omega }0`$), as discussed in Sec. III A above.

### B Classes of modes

We next review the various different classes of modes to which we shall apply the selection rules. In zero-buoyancy stars, we distinguish three different classes of modes:

* The rotational modes which are axial to zeroth order in $`\mathrm{\Omega }`$. These modes are characterized by single values of $`l`$ and of $`m`$. If one is working with all the modes $`𝝃_A`$, then one can have $`m=l`$ or $`m=l`$. If one is working only with the set of modes $`\xi _\alpha `$ with positive rotating-frame frequency $`\omega `$ (the remaining modes being complex conjugates of these), then only $`m=l`$ is allowed. These modes are all $`z`$-parity odd. We shall call these modes the pure r-modes.
* The rotational modes whose $`\mathrm{\Omega }0`$ limits are not axial. These have fixed $`m`$ but do not have a fixed value of $`l`$. To zeroth order in $`\mathrm{\Omega }`$ they can be expanded as

$`𝝃(𝐱)={\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{W_{j+m+1}(r)}{r}}Y_{j+m+1}^m(r)𝒆_{\widehat{r}}+V_{j+m+1}(r)\mathbf{}Y_{j+m+1}^m{\displaystyle \frac{U_{j+m}(r)}{r}}𝐋Y_{j+m}^m\right],`$ (149)
in the $`z`$-parity odd case, and as

$`𝝃(𝐱)={\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{W_{j+m}(r)}{r}}Y_{j+m}^m𝒆_{\widehat{r}}+V_{j+m}(r)\mathbf{}Y_{j+m}^m{\displaystyle \frac{U_{j+m+1}(r)}{r}}𝐋Y_{j+m+1}^m\right]`$ (150)
in the $`z`$parity even case, where $`𝐋=i(𝐫\times \mathbf{})`$ and $`j`$ runs over $`0,2,4,\mathrm{}`$. Following Ref.  we shall call these the rotational hybrid modes, although they have also been called inertial modes or generalized $`r`$-modes.
* The regular modes, which can be $`z`$-parity even or $`z`$-parity odd, and consist of $`f`$ and $`p`$ modes.

Similarly, in stars with buoyancy we shall distinguish the following classes of modes:

* The rotational modes. All of these modes are axial to zeroth order in $`\mathrm{\Omega }`$ (Appendix E). They are also characterized by single values of $`l`$ and $`m`$, but any values of $`l`$, $`m`$ with $`l|m|`$ are allowed. They have $`z`$-parity $`(1)^{l+m+1}`$ and parity $`(1)^{l+1}`$. We shall call these the $`r`$-modes. \[Unlike the pure $`r`$-modes above, their radial eigenfunctions are not polynomial in $`r`$.\]
* The regular modes, which may be $`z`$-parity even or $`z`$-parity odd, and consist of $`f`$, $`p`$ and $`g`$ modes.

### C Selection rules

We can write the three mode coupling coefficient as a symmetric function of the three mode functions $`𝝃_A`$, $`𝝃_B`$, $`𝝃_C`$:

$$\kappa _{ABC}=\kappa (𝝃_A,𝝃_B,𝝃_C).$$
(151)
Now apply the pullback $`f_{}`$ to this equation. Since the left hand side is a pure number, it is invariant. On the right hand side, we can pull the $`f_{}`$ operator inside the function $`\kappa `$ since pullback commutes with geometrical operations like taking covariant derivatives. Thus

$$\kappa _{ABC}=\kappa (f_{}𝝃_A,f_{}𝝃_B,f_{}𝝃_C).$$
(152)
If we denote the $`z`$-parity of mode $`A`$ by $`ϵ_A=\pm 1`$, so that $`f_{}𝝃_A=ϵ_A𝝃_A`$, then it follows from Eqs. (151) and (152) that

$$\left[1ϵ_Aϵ_Bϵ_C\right]\kappa _{ABC}=0.$$
(153)
Therefore we get the selection rule

$$\text{z}parityselectionrule:\text{ odd number of z-parity odd modes }\kappa _{ABC}=0.$$
(154)
For example, if all three modes are $`z`$-parity odd, then the coupling coefficient vanishes. An identical argument gives the selection rule for regular parity

$$\text{ }parityselectionrule:\mathrm{odd}\mathrm{number}\mathrm{of}\mathrm{parity}\mathrm{odd}\mathrm{modes}\kappa _{ABC}=0.$$
(155)
Finally, applying Eq. (153) with $`f`$ replaced by the mapping $`\phi \phi +\mathrm{\Delta }\phi `$ and with $`ϵ_A`$ replaced by $`\mathrm{exp}[im_A\mathrm{\Delta }\phi ]`$ \[cf. Eq. (148) above\] yields

$$m\text{s}electionrule:m_A+m_B+m_C0\kappa _{ABC}=0.$$
(156)
These selection rules are valid in rapidly rotating stars as well as in the slow rotation limit, and also retain their validity for relativistic stars. The arguments can also be generalized to four-mode coupling coefficients $`\kappa _{ABCD}`$.

The selection rules (154), (155), (156) are based on symmetries of the background star and thus apply to all orders in $`\mathrm{\Omega }`$. Therefore, when they are applicable, they restrict not only the piece of $`\kappa _{ABC}`$ which is zeroth order in $`\mathrm{\Omega }`$ (which is all we compute in this paper), but also all correction terms that are higher order in $`\mathrm{\Omega }`$. By contrast, in Appendix G we derive two selection rules which are valid to zeroth order in $`\mathrm{\Omega }`$ only. First, we show the coupling coefficient for the coupling of three axial modes is $`O(\mathrm{\Omega })`$. Thus we have

$$\text{ }Axialselectionrule1:𝝃_A,𝝃_B,𝝃_C=\mathrm{axial}+O(\mathrm{\Omega })\kappa _{ABC}=O(\mathrm{\Omega }).$$
(157)
We also prove in Appendix G the following second selection rule. If (i) two of the modes are axial, (ii) the third has vanishing Eulerian density perturbation, and (iii) the background star and perturbations obey the same one-parameter equation of state (i.e., the zero-buoyancy or barotropic case), then $`\kappa _{ABC}=O(\mathrm{\Omega })`$. We can write this selection rule as

$$\text{ }Axialselectionrule2:\delta \rho _A=O(\mathrm{\Omega }),𝝃_B,𝝃_C=\mathrm{axial}+O(\mathrm{\Omega }),\mathrm{zero}\mathrm{buoyancy}\kappa _{ABC}=O(\mathrm{\Omega }).$$
(158)
Note that axial modes automatically have vanishing density perturbation to zeroth order in $`\mathrm{\Omega }`$, from Eq. (3) and using $`𝝃=\xi ^r=0`$ and the fact that the background density is a function of $`r`$ only to leading order in $`\mathrm{\Omega }`$. The selection rules (157) and (158) might fail for relativistic stars.

In incompressible stars, we prove in Appendix G that the coupling between any three rotational modes vanishes to zeroth order in $`\mathrm{\Omega }`$. A precise statement of this selection rule is

$$\text{ }Incompressibleselectionrule:\delta \rho _A,\delta p_A,\delta \rho _B,\delta p_B,\delta \rho _C,\delta p_C=O(1/\mathrm{\Gamma }_1)+O(\mathrm{\Omega })\kappa _{ABC}=O(1/\mathrm{\Gamma }_1)+O(\mathrm{\Omega }),$$
(159)
where $`\mathrm{\Gamma }_1`$ is the adiabatic index of the perturbations. Note that the regular modes ($`f`$-modes) of an incompressible star have $`\delta \rho \mathrm{const}`$ in the limit $`\mathrm{\Gamma }_1\mathrm{}`$, and do not satisfy the hypothesis of this selection rule, whereas the rotational modes have exactly vanishing $`\delta \rho `$ and $`\delta p`$ (to zeroth order in $`\mathrm{\Omega }`$) and so do satisfy the hypothesis.

The $`\mathrm{\Omega }`$ dependence of the coupling coefficients $`\kappa _{ABC}`$ arises in two ways: (i) through the dependence of the expression (134) on the background star, for which the leading order corrections arise at $`O(\mathrm{\Omega }^2)`$, and (ii) through the dependence of the mode functions themselves on $`\mathrm{\Omega }`$, for which the leading order corrections are $`O(\mathrm{\Omega }^2)`$ for pure $`r`$-modes , but can be $`O(\mathrm{\Omega })`$ for regular modes. Therefore, when the zeroth order coupling coefficient vanishes, the leading order corrections may be $`O(\mathrm{\Omega })`$ or $`O(\mathrm{\Omega }^2)`$, depending on the modes involved.

### D Applications of selection rules

We now discuss some applications of these rules. Consider first the case of zero-buoyancy stars. The most interesting coupling coefficients are those in which an unstable pure $`r`$-mode appears twice, since if the amplitudes of all the other modes are zero initially, it is via these coupling coefficients that energy can leak out of the pure $`r`$ mode. In particular we are interested in the most unstable mode, the $`l=m=2`$ mode. If $`𝝃_A=𝝃_B`$ is an unstable pure $`r`$-mode, then the third mode $`𝝃_C`$ must be a $`z`$-parity even mode from rule (154). In particular, the third mode cannot be a pure $`r`$-mode as pure $`r`$ modes have odd $`z`$ parity. Coupling coefficients between any three pure $`r`$-modes vanish. Nonzero coefficients can be obtained by taking the third mode to be a $`z`$-parity even regular mode. However, if the third mode is taken to be a $`z`$-parity even rotational hybrid, then the second axial selection rule (158) applies due to Eq. (108), and implies that the coupling coefficient vanishes to zeroth order in $`\mathrm{\Omega }`$. It may be nonvanishing to higher order in $`\mathrm{\Omega }`$ . If one considers one pure $`r`$-mode and two rotational hybrid modes, the corresponding coupling vanishes in an incompressible star by the incompressible selection rule, but may be nonvanishing at $`O(1/\mathrm{\Gamma }_1)`$.

In Table I we summarize which classes of modes can and cannot couple to zeroth order in $`\mathrm{\Omega }`$ in in zero-buoyancy stars. Stars indicate non-zero coupling coefficients. Since $`\kappa _{ABC}`$ is symmetric, Table I lists all possible types of coupling coefficient except one, the coupling coefficient between one pure $`r`$-mode, one regular mode and one hybrid rotational mode; this coupling coefficient is non-zero. The table shows the kinds of coupling coefficients that can be non-zero when $`z`$parity restrictions are met. For example, the table indicates that three hybrid rotational modes can couple together, but this can only happen when the coupling involves an even number of $`z`$parity odd modes.

For nonzero-buoyancy stars, consider the coupling coefficient between three $`r`$-modes. If the three modes have quantum numbers $`l_A`$, $`m_A`$, $`l_B`$, $`m_B`$ and $`l_C`$, $`m_C`$, then since each mode has $`z`$-parity $`(1)^{l+m+1}`$, Eq. (153) gives

$$\left[1(1)^{3+m_A+m_B+m_C+l_A+l_B+l_C}\right]\kappa _{ABC}=0.$$
(160)
Invoking the $`m`$ selection rule (156) implies that $`l_A+l_B+l_C`$ must be odd for $`\kappa _{ABC}`$ to be non-zero. If all three modes are of the $`|m|=l`$ variety, this contradicts $`m_A+m_B+m_C=0`$, so the coupling of three $`l=|m|`$ modes is zero. If some of the three modes do not satisfy $`l=|m|`$, then the $`z`$-parity and parity selection rules do not apply. However, the axial selection rule (157) does apply and therefore the coupling coefficient vanishes to zeroth order in $`\mathrm{\Omega }`$. To higher order in $`\mathrm{\Omega }`$ however, these coupling coefficients with $`l_A+l_B+l_C`$ odd and $`m_A+m_B+m_C=0`$ need not vanish. In fact, Morsink  has calculated coupling coefficients of this type to $`O(\mathrm{\Omega }^2)`$ and has found them to be nonzero.

## VI Computational method for Coupling Coefficients

In this section we describe an efficient computational method for computing coupling coefficients in rotating stars. Specific coupling coefficients relevant to $`r`$-mode saturation will be computed numerically in a subsequent paper.

### A Overview

Direct calculation of the expression (134) for the coupling coefficients is algebraically intensive even in the case where all three modes have no sums over $`l`$. If one proceeds in spherical coordinates to directly calculate all the covariant derivatives needed one encounters many terms that scale as $`1/\mathrm{sin}^n\theta `$ with $`n2`$. These terms are divergent as $`\theta 0`$ or $`\pi `$. Since the full coupling coefficient is obviously finite at the poles ($`\theta =0,\pi `$), these terms must cancel with other terms of the same order in $`1/\mathrm{sin}\theta `$. However, the amount of algebra necessary to cancel these terms by hand is unmanageable. Calculation of the coupling coefficients in Cartesian coordinates is also daunting. Even though all covariant derivatives are highly simplified in Cartesian coordinates, the total number of terms is greater than in spherical coordinates and one must still integrate over the sphere in the end, which is awkward. Cartesian coordinates also obscure selection rules that might lead to insight into which modes couple and which don’t. Methods based on integration by parts simplify the algebra for non-rotating stars , but, for rotating stars, modes acquire toroidal pieces that severely complicate the angular integrations even after integration by parts.

To circumvent these difficulties, we have found that using the spin-weighted spherical harmonics of Ref.  as angular basis functions is ideal. The use of these basis functions allows one to immediately reduce the angular integrals to Wigner $`3j`$ symbols. It also reduces dramatically the number of terms that must be calculated. For example, for modes that have only one value of $`l`$ each, (134) requires explicit calculation of $`100`$ terms, including many that require computing angular derivatives of the mode functions. When modes involve sums over $`n`$ values of $`l`$, the number of terms to calculate scales as $`n^3`$. Once this is all done, one must then integrate the resulting complicated expression over angles. However, in the spin-weighted formalism, the calculation involves just $`6n^3`$ terms, none of which requires computing angular derivatives. Moreover, those $`6n^3`$ terms are completely spherically symmetric, all angular dependence being already integrated out in terms of Wigner $`3j`$ symbols. The coupling coefficients are then built from these spherically symmetric terms, their corresponding $`3j`$ symbols, and their complex conjugates and symmetrizations. The method is easy to automate, including the various sums over the different $`l`$’s of the modes, leaving only $`6`$ terms to actually code up.

This section is organized as follows. In Sec. VI B we define the spin-weighted basis and express the covariant derivatives needed for the coupling coefficients in terms of components on that basis. Some details are relegated to Appendix H. Section VI C derives the explicit expression (213) for the coupling coefficient in terms of spherically symmetric functions and Wigner $`3j`$ symbols.

The method described here is specialized to slowly rotating stars and is valid only to zeroth order in $`\mathrm{\Omega }`$, since we take the background density and pressure to be functions of $`r`$ only. However, it is straightforward to extend the method to higher orders in $`\mathrm{\Omega }`$. For example, in calculating coupling coefficients to second order in $`\mathrm{\Omega }`$, contributions of order $`O(\mathrm{\Omega }^2)`$ will arise both from corrections to the background stellar model (centrifugal flattening), and from corrections to the modes. Including the additional angular dependence of the background density and pressure would require computing angular integrals of products of four spin-weighted spherical harmonics instead of three.

### B An appropriate choice of basis

Consider a given mode function $`𝝃(𝐱)`$. Since any mode has a fixed value of $`m`$, we can assume an expansion of the mode function on a basis of vector spherical harmonics of the form

$`𝝃={\displaystyle \underset{l=|m|}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{W_{lm}(r)}{r}}Y_{lm}𝒆_{\widehat{r}}+V_{lm}(r)\mathbf{}Y_{lm}{\displaystyle \frac{U_{lm}(r)}{r}}𝐋Y_{lm}\right],`$ (161)
where $`𝐋=i(𝐫\times \mathbf{})`$. We define a set of spin-weighted basis vectors consisting of two complex null 3-vectors, $`𝐦`$ and $`\overline{𝐦}`$ given by

$$𝐦=\frac{1}{\sqrt{2}}\left(𝐞_{\widehat{\theta }}+i𝐞_{\widehat{\phi }}\right)$$
(162)
and

$$\overline{𝐦}=\frac{1}{\sqrt{2}}\left(𝐞_{\widehat{\theta }}i𝐞_{\widehat{\phi }}\right)$$
(163)
and a vector orthogonal to those two,

$$𝐥=𝐞_{\widehat{r}}.$$
(164)
The components of the metric on this non-coordinate basis are

$`g^{ll}`$ $`=g_{ll}=1,`$ (165)
$`g^{m\overline{m}}`$ $`=g_{m\overline{m}}=1,`$ (166)
with all other components zero. The components of the mode function $`𝝃`$ on this basis are given by

$`\xi _l`$ $``$ $`\vartheta _0=𝐥𝝃`$ (167)
$`\xi _m`$ $``$ $`\vartheta _1=𝐦𝝃`$ (168)
$`\xi _{\overline{m}}`$ $``$ $`\vartheta _1=\overline{𝐦}𝝃,`$ (169)
where the subscript $`0,1`$ or $`1`$ on the scalar function $`\vartheta (r,\theta ,\phi )`$ represents the spin weight†§†§†§For a definition of spin-weight, see Appendix H. of the component. Note that

$`\xi _l`$ $`=`$ $`\xi ^l=\vartheta _0`$ (170)
$`\xi _m`$ $`=`$ $`\xi ^{\overline{m}}=\vartheta _1`$ (171)
$`\xi _{\overline{m}}`$ $`=`$ $`\xi ^m=\vartheta _1,`$ (172)
and that the mode function can be reconstructed as

$$𝝃=\xi ^l𝐥+\xi ^{\overline{m}}\overline{𝐦}+\xi ^m𝐦.$$
(173)
Spin-weighted spherical harmonics are defined in terms of the matrix representations of the rotation group, $`D_{sm}^l`$, by

$${}_{s}{}^{}Y_{lm}^{}(\theta ,\phi )=\sqrt{(2l+1)/4\pi }D_{sm}^l(\varphi ,\theta ,0),$$
(174)
where $`s`$ is the spin weight and $`{}_{0}{}^{}Y_{lm}^{}`$ are the ordinary spherical harmonics $`Y_{lm}`$. See Ref.  for a definition and a detailed list of the properties of the quantities $`D_{sm}^l`$. For each value of the spin weight $`s`$, these spin-weighted spherical harmonics form a complete orthonormal set, i.e.

$${}_{s}{}^{}Y_{lm}^{}{}_{}{}^{}{}_{s}{}^{}Y_{l^{}m^{}}^{}𝑑\mathrm{\Omega }=\delta _{l^{}l}\delta _{m^{}m}.$$
(175)
Thus we can expand $`\vartheta _0,\vartheta _1,\vartheta _1`$ on the spin-weighted spherical harmonics bases as

$`\vartheta _0(r,\theta ,\phi )`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}f_0^\mathrm{\Lambda }(r){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}(\theta ,\phi )`$ (176)
$`\vartheta _1(r,\theta ,\phi )`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}f_+^\mathrm{\Lambda }(r){}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}(\theta ,\phi )`$ (177)
$`\vartheta _1(r,\theta ,\phi )`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}f_{}^\mathrm{\Lambda }(r){}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}(\theta ,\phi ),`$ (178)
where the subscripts $`+`$ and $``$ are understood to mean $`+1`$ and $`1`$ respectively, and $`\mathrm{\Lambda }`$ denotes $`(lm)`$. Therefore the expansion of the mode function can be written as

$$𝝃=\underset{lm}{}\left[f_0^{lm}(r){}_{0}{}^{}Y_{lm}^{}𝐥+f_+^{lm}(r){}_{1}{}^{}Y_{lm}^{}\overline{𝐦}+f_{}^{lm}(r){}_{1}{}^{}Y_{lm}^{}𝐦\right].$$
(179)
In Appendix H we relate the coefficients $`f_s^{lm}(r)`$ appearing in this expansion to the coefficients $`U_{lm}(r)`$, $`V_{lm}(r)`$, and $`W_{lm}(r)`$ appearing in the expansion (161) above \[see Eqs. (H32) and (H37)\].

We will be using a nonstandard notational convention involving the symbol $`\mathrm{\Lambda }`$ which we now explain. In rotating stars, $`l`$ is generally not a good “quantum number” while $`m`$ is. In general there will be other good quantum numbers that will allow us to enumerate the modes. What is meant by the symbol $`\mathrm{\Lambda }`$ is a set of numbers representing both the good and bad quantum numbers of the object in question. For example, when $`\mathrm{\Lambda }`$ is used on a $`{}_{s}{}^{}Y_{lm}^{}`$ we mean the string of quantum numbers $`(l,m)`$. However, when $`\mathrm{\Lambda }`$ is used on the mode function $`f_+^\mathrm{\Lambda }`$ it may mean a string of numbers that includes, but is not limited to $`l`$ and $`m`$. When we sum over $`\mathrm{\Lambda }`$ it is implied that we sum only over bad quantum numbers. Thus we might label the above component $`\vartheta _0`$ with another index ($`A`$ say) representing only the good quantum numbers that are left after the bad ones have been summed over, e.g.

$$\vartheta _0^A(r,\theta ,\phi )=\underset{\mathrm{\Lambda }}{}f_0^\mathrm{\Lambda }(r){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}(\theta ,\phi ).$$
(180)
And finally, when we use the symbol $`\mathrm{\Lambda }`$ not in a sub- or super- script it is defined to be $`\mathrm{\Lambda }\sqrt{l(l+1)}`$.

We now express the covariant derivatives of the mode functions that appear in the coupling coefficient (134) in terms of our new expansion coefficients $`f_s^\mathrm{\Lambda }(r)`$. The details of the calculation are given in Appendix H. We define the functionals $`G_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$, $`F_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$, and $`H_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$ to be

$`G_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$ $``$ $`{\displaystyle \frac{1}{r}}\left(f_0^\mathrm{\Lambda }(r){\displaystyle \frac{s\mathrm{\Lambda }}{\sqrt{2}}}f_s^\mathrm{\Lambda }(r)\right)`$ (181)
$`F_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$ $``$ $`{\displaystyle \frac{1}{r}}\left(f_s^\mathrm{\Lambda }(r)+{\displaystyle \frac{s\mathrm{\Lambda }}{\sqrt{2}}}f_0^\mathrm{\Lambda }(r)\right)`$ (182)
$`H_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$ $``$ $`{\displaystyle \frac{s}{\sqrt{2}r}}f_s^\mathrm{\Lambda }(r)\sqrt{\mathrm{\Lambda }^22}.`$ (183)
Then the components in the basis $`(𝐥,𝐦,\overline{𝐦})`$ of the tensor $`\xi _{;j}^i`$ can be written as

$`\xi _{;l}^l`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}f_{0}^{\mathrm{\Lambda }}{}_{,r}{}^{}(r){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (184)
$`\xi _{;l}^{\overline{m}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}f_{+}^{\mathrm{\Lambda }}{}_{,r}{}^{}(r){}_{1}{}^{}Y_{\mathrm{\Lambda }}^{},\xi _{;l}^m={\displaystyle \underset{\mathrm{\Lambda }}{}}f_{}^{\mathrm{\Lambda }}{}_{,r}{}^{}(r){}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (185)
$`\xi _{;m}^m`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}G_+^\mathrm{\Lambda }[f]{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{},\xi _{;\overline{m}}^{\overline{m}}={\displaystyle \underset{\mathrm{\Lambda }}{}}G_{}^\mathrm{\Lambda }[f]{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{},`$ (186)
$`\xi _{;m}^{\overline{m}}`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}H_+^\mathrm{\Lambda }[f]{}_{2}{}^{}Y_{\mathrm{\Lambda }}^{},\xi _{;\overline{m}}^m={\displaystyle \underset{\mathrm{\Lambda }}{}}H_{}^\mathrm{\Lambda }[f]{}_{2}{}^{}Y_{\mathrm{\Lambda }}^{},`$ (187)
$`\xi _{;m}^l`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}F_+^\mathrm{\Lambda }[f]{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{},\xi _{;\overline{m}}^l={\displaystyle \underset{\mathrm{\Lambda }}{}}F_{}^\mathrm{\Lambda }[f]{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{},`$ (188)
where we abbreviate $`G_s^\mathrm{\Lambda }[f_s^\mathrm{\Lambda }(r)]`$ as $`G_s^\mathrm{\Lambda }[f]`$, etc. Note that this formalism allows us to write down forms for the covariant derivatives in terms of functions of $`r`$ only. The angular derivatives have all been computed. Thus the formalism avoids the large curvature terms that can arise from angular covariant derivatives in spherical coordinates. Another advantage is that the formalism allows us to integrate the coupling coefficients over angles immediately; no integrations by parts are needed. To see how this works we note that the integral over angles of three of the spin-weighted spherical harmonics can be written simply in terms of Wigner $`3j`$ symbols , which are easily calculated:

$`{\displaystyle {}_{s_1}{}^{}Y_{l_1m_1}^{}{}_{s_2}{}^{}Y_{l_2m_2}^{}{}_{s_3}{}^{}Y_{l_3m_3}^{}𝑑\mathrm{\Omega }}`$ $``$ $`{}_{s_1}{}^{}Y_{l_1m_1}^{}{}_{s_2}{}^{}Y_{l_2m_2}^{}{}_{s_3}{}^{}Y_{l_3m_3}^{}`$ (189)
$`=`$ $`4\pi \sqrt{(2l_1+1)/4\pi }\sqrt{(2l_2+1)/4\pi }\sqrt{(2l_3+1)/4\pi }`$ (191)
$`\times {\displaystyle }D_{s_1m_1}^{l_1}(\phi ,\theta ,0)D_{s_2m_2}^{l_2}(\phi ,\theta ,0)D_{s_3m_3}^{l_3}(\phi ,\theta ,0)d\mathrm{\Omega }`$
$`=`$ $`4\pi \sqrt{(2l_1+1)/4\pi }\sqrt{(2l_2+1)/4\pi }\sqrt{(2l_3+1)/4\pi }`$ (197)
$`\times \left(\begin{array}{ccc}l_1& l_2& l_3\\ s_1& s_2& s_3\end{array}\right)\left(\begin{array}{ccc}l_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right),`$
Note that this integral vanishes unless $`m_1+m_2+m_3=0`$, cf. the $`m`$ selection rule (156).

### C Integrating the coupling coefficients over the sphere

To begin integrating the expression (134) over angles we write

$`\kappa _{ABC}{\displaystyle r^2𝑑r\kappa _{ABC}(r)},`$ (198)
and note that in calculating $`\kappa _{ABC}`$ we can ignore mode indices as long as we symmetrize at the end of the computation. To this end define an operator $`𝒮`$ that symmetrizes over mode indices:

$`𝒮L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}L_{(\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3)}`$ $``$ $`{\displaystyle \frac{1}{6}}(L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}+L_{\mathrm{\Lambda }_2\mathrm{\Lambda }_1\mathrm{\Lambda }_3}`$ (200)
$`+L_{\mathrm{\Lambda }_3\mathrm{\Lambda }_2\mathrm{\Lambda }_1}+L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_3\mathrm{\Lambda }_2}+L_{\mathrm{\Lambda }_3\mathrm{\Lambda }_1\mathrm{\Lambda }_2}+L_{\mathrm{\Lambda }_2\mathrm{\Lambda }_3\mathrm{\Lambda }_1}).`$
Now expand the scalar $`\mathbf{}𝝃_A`$ as

$`\mathbf{}𝝃_A{\displaystyle \underset{\mathrm{\Lambda }}{}}g_0^\mathrm{\Lambda }(r){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{},`$ (201)
and the scalar $`\delta \varphi _A`$ as

$`a\delta \varphi _A{\displaystyle \underset{\mathrm{\Lambda }}{}}a_0^\mathrm{\Lambda }(r){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}.`$ (202)
The second covariant derivatives of this scalar are then given by

$`a_{;mm}`$ $`=`$ $`{\displaystyle \underset{\mathrm{\Lambda }}{}}E_0^\mathrm{\Lambda }[a]{}_{2}{}^{}Y_{\mathrm{\Lambda }}^{},a_{;\overline{m}\overline{m}}={\displaystyle \underset{\mathrm{\Lambda }}{}}E_0^\mathrm{\Lambda }[a]{}_{2}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (203)
$`a_{;lm}`$ $`=`$ $`a_{;ml}={\displaystyle \underset{\mathrm{\Lambda }}{}}C_+^\mathrm{\Lambda }[a]{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{},a_{;l\overline{m}}=a_{;\overline{m}l}={\displaystyle \underset{\mathrm{\Lambda }}{}}C_{}^\mathrm{\Lambda }[a]{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (204)
$`a_{;m\overline{m}}`$ $`=`$ $`a_{;\overline{m}m}={\displaystyle \underset{\mathrm{\Lambda }}{}}D_0^\mathrm{\Lambda }[a]{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{},a_{;ll}={\displaystyle \underset{\mathrm{\Lambda }}{}}a_{0}^{\mathrm{\Lambda }}{}_{,rr}{}^{}{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{},`$ (205)
where the functionals $`E_0^\mathrm{\Lambda }[a(r)],C_s^\mathrm{\Lambda }[a(r)]`$ and $`D_0^\mathrm{\Lambda }[a(r)]`$ are

$`E_0^\mathrm{\Lambda }[a(r)]`$ $``$ $`{\displaystyle \frac{\mathrm{\Lambda }}{2r^2}}\sqrt{\mathrm{\Lambda }^22}a_0^\mathrm{\Lambda }(r),`$ (206)
$`D_0^\mathrm{\Lambda }[a(r)]`$ $``$ $`{\displaystyle \frac{1}{2r^2}}\left(2ra_{0}^{\mathrm{\Lambda }}{}_{,r}{}^{}a_0^\mathrm{\Lambda }\mathrm{\Lambda }^2\right),`$ (207)
$`C_s^\mathrm{\Lambda }[a(r)]`$ $``$ $`{\displaystyle \frac{s\mathrm{\Lambda }}{\sqrt{2}r}}\left({\displaystyle \frac{1}{r}}a_0^\mathrm{\Lambda }a_{0}^{\mathrm{\Lambda }}{}_{,r}{}^{}\right).`$ (208)
We are now in a position to write down $`\kappa _{ABC}(r)`$:

$`\kappa _{ABC}(r)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}{}}\{L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(1)}(r){}_{0}{}^{}Y_{\mathrm{\Lambda }_1}^{}{}_{0}{}^{}Y_{\mathrm{\Lambda }_2}^{}{}_{0}{}^{}Y_{\mathrm{\Lambda }_3}^{}`$ (213)
$`+6𝒮\left[L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(2)}(r){}_{0}{}^{}Y_{\mathrm{\Lambda }_1}^{}{}_{1}{}^{}Y_{\mathrm{\Lambda }_2}^{}{}_{1}{}^{}Y_{\mathrm{\Lambda }_3}^{}\right]`$
$`+6𝒮\left[L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(3)}(r){}_{0}{}^{}Y_{\mathrm{\Lambda }_1}^{}{}_{2}{}^{}Y_{\mathrm{\Lambda }_2}^{}{}_{2}{}^{}Y_{\mathrm{\Lambda }_3}^{}\right]`$
$`+6𝒮\left[L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(4)}(r){}_{2}{}^{}Y_{\mathrm{\Lambda }_1}^{}{}_{1}{}^{}Y_{\mathrm{\Lambda }_2}^{}{}_{1}{}^{}Y_{\mathrm{\Lambda }_3}^{}\right]`$
$`+(+)\},`$
where the $`(+)`$ symbol means repeat all terms in the above equation replacing $`ss`$ in the integrals and interchanging $`+`$ and $``$ subscripts in the functions (detailed below). The functions $`L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(i)}(r)`$ are given by

$`L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(1)}(r)`$ $`=`$ $`{\displaystyle \frac{1}{2}}p\left\{(\mathrm{\Gamma }_11)^2+{\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{ln}\rho _0}}\right\}g_0^{\mathrm{\Lambda }_1}g_0^{\mathrm{\Lambda }_2}g_0^{\mathrm{\Lambda }_3}{\displaystyle \frac{1}{2}}\rho {\displaystyle \frac{d^3\varphi }{dr^3}}f_0^{\mathrm{\Lambda }_1}f_0^{\mathrm{\Lambda }_2}f_0^{\mathrm{\Lambda }_3}`$ (217)
$`+𝒮\{p(\mathrm{\Gamma }_11)3[{\displaystyle \frac{1}{2}}g_0^{\mathrm{\Lambda }_1}f_{0}^{\mathrm{\Lambda }_2}{}_{,r}{}^{}f_{0}^{\mathrm{\Lambda }_3}{}_{,r}{}^{}+g_0^{\mathrm{\Lambda }_1}G_+^{\mathrm{\Lambda }_2}[f]G_+^{\mathrm{\Lambda }_3}[f]]`$
$`+2p\left[{\displaystyle \frac{1}{2}}f_{0}^{\mathrm{\Lambda }_1}{}_{,r}{}^{}f_{0}^{\mathrm{\Lambda }_2}{}_{,r}{}^{}f_{0}^{\mathrm{\Lambda }_3}{}_{,r}{}^{}+G_+^{\mathrm{\Lambda }_1}[f]G_+^{\mathrm{\Lambda }_2}[f]G_+^{\mathrm{\Lambda }_3}[f]\right]`$
$`\rho {\displaystyle \frac{3}{2}}f_0^{\mathrm{\Lambda }_1}a_{0}^{\mathrm{\Lambda }_2}{}_{,rr}{}^{}f_0^{\mathrm{\Lambda }_3}\},`$
$`L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(2)}(r)`$ $`=`$ $`p(\mathrm{\Gamma }_11)g_0^{\mathrm{\Lambda }_1}F_+^{\mathrm{\Lambda }_2}[f]f_{}^{\mathrm{\Lambda }_3}{}_{,r}{}^{}`$ (221)
$`\rho \left({\displaystyle \frac{1}{2}}D_0^{\mathrm{\Lambda }_1}[a]f_+^{\mathrm{\Lambda }_2}f_{}^{\mathrm{\Lambda }_3}+f_0^{\mathrm{\Lambda }_1}C_+^{\mathrm{\Lambda }_2}[a]f_{}^{\mathrm{\Lambda }_3}\right)`$
$`+p\left[f_{0}^{\mathrm{\Lambda }_1}{}_{,r}{}^{}+G_+^{\mathrm{\Lambda }_1}[f]\right]F_+^{\mathrm{\Lambda }_2}[f]f_{}^{\mathrm{\Lambda }_3}{}_{,r}{}^{}`$
$`{\displaystyle \frac{1}{2}}\rho {\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{d\varphi }{dr}}\right)f_0^{\mathrm{\Lambda }_1}f_+^{\mathrm{\Lambda }_2}f_{}^{\mathrm{\Lambda }_3},`$
$`L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(3)}(r)`$ $`=`$ $`{\displaystyle \frac{1}{2}}p(\mathrm{\Gamma }_11)g_0^{\mathrm{\Lambda }_1}H_+^{\mathrm{\Lambda }_2}[f]H_{}^{\mathrm{\Lambda }_3}[f]+pG_+^{\mathrm{\Lambda }_1}[f]H_+^{\mathrm{\Lambda }_2}[f]H_{}^{\mathrm{\Lambda }_3}[f],`$ (222)
$`L_{\mathrm{\Lambda }_1\mathrm{\Lambda }_2\mathrm{\Lambda }_3}^{(4)}(r)`$ $`=`$ $`pH_{}^{\mathrm{\Lambda }_1}[f]f_{+}^{\mathrm{\Lambda }_2}{}_{,r}{}^{}F_+^{\mathrm{\Lambda }_3}[f]\rho {\displaystyle \frac{1}{2}}E_0^{\mathrm{\Lambda }_1}[a]f_+^{\mathrm{\Lambda }_2}f_+^{\mathrm{\Lambda }_3},`$ (223)
where a comma indicates ordinary differentiation.

To summarize, given a set of three modes $`𝝃_A`$, $`𝝃_B`$, $`𝝃_C`$, the computational method is to (i) Use Eqs. (H32) and (H37) to compute the expansion coefficients $`f_s^{lm}(r)`$ for each mode in terms of the more standard expansion coefficients $`U_{lm}(r)`$, $`V_{lm}(r)`$, and $`W_{lm}(r)`$ appearing in the expansion (161), and (ii) Use the expressions (197), (198) and (213) – (223) to compute the coupling coefficient $`\kappa _{ABC}`$ from the functions $`f_s^{lm}(r)`$.

## VII Summary

In this paper, we have formulated a perturbative approach to the nonlinear interactions of unstable $`r`$modes in a neutron star. Our formalism presumes that mode growth saturates at moderately low amplitudes, so that we can model the modal interactions via three-mode couplings. By developing further a previous perturbation theory for rotating stars, we have found equations of motion for the mode amplitudes of rotating stars that are uncoupled at linear order when acted upon by an external force. This feature of our formalism is essential for following the cascade of energy from one mode to another when lowest order nonlinear couplings are included.

The important astrophysical question is what determines the saturation amplitude for unstable $`r`$modes. Fundamentally, if mode-mode coupling is the dominant damping, we expect the saturation amplitude to be set by a competition between the growth rate of the instability, and the amplitude-dependent rate of drainage of energy from an unstable mode to other stellar modes. Our formalism gives explicit formulae for the lowest-order coupling coefficients among stellar modes. The numerical simulations by Stergioulas and Font , and Lindblom, Tohline and Vallisneri , have found the coupling of the unstable $`l=|m|=2`$ $`r`$mode to other stellar modes to be surprisingly weak, resulting in a timescale for energy transfer to other modes that is at least $`20`$ rotation periods even when the modal amplitude is substantial.

Although we shall present detailed numerical solutions of our equations of motion for modal amplitudes in a subsequent publication, there are already hints, in the formal developments presented here, of an explanation for the apparent weakness of the coupling of the $`l=|m|=2`$ $`r`$mode to other stellar modes. Strong coupling requires either near-resonance or a large coupling coefficient, or both. Although the rotational modes have frequencies comparable to the $`r`$-mode, we have found that the possible couplings to such modes may be limited. For example, in zero-buoyancy stars parity arguments prevent the coupling of three $`r`$modes, and the couplings involving two $`r`$ modes and one hybrid rotational mode vanish to zeroth order in the stellar angular velocity. In nonzero-buoyancy stars, the coupling between $`r`$-modes is again vanishing to zeroth order in $`\mathrm{\Omega }`$, but is nonzero when correction terms of order $`O(\mathrm{\Omega }^2)`$ are included . Thus, the couplings among $`r`$modes are small, generically. The coupling of two $`r`$modes to a $`f`$mode need not be small. However, since the eigenfrequencies of $`f`$modes are $`(GM/R^3)^{1/2}`$, which is much larger than $`\mathrm{\Omega }`$ for a slowly rotating star, the excitation of an $`f`$mode from small amplitude by an unstable $`r`$mode is likely to be suppressed, by factors $`\mathrm{\Omega }^2R^3/GM`$.

Finally, as explained in the introduction, the possibility of the $`r`$-modes saturating at small amplitude is not necessarily incompatible with the numerical simulations to date, or with the apparent weakness of the nonlinear couplings of the $`r`$-modes. If the modes do saturate at small or moderate amplitude, then the formalism developed in this paper may be sufficient to explore the saturation process. On the other hand, it is possible that the formalism developed here may have to be developed further in order to understand the saturation. Some possibilities are

* It may be that the $`r`$-modes stop growing only in the very nonlinear regime, where shocks develop, as in the simulations of Lindblom, Tohline and Vallisneri . If this is the case then the perturbative formalism developed here will not be useful. However, as explained in the introduction, the simulations to date do not show that the strongly nonlinear regime is reached.
* It may be that an analysis with only three-mode couplings will be insufficient, as cubic potentials generically have instabilities at large amplitude , but that including four-mode couplings will be sufficient to allow an exploration of the saturation process.
* This paper has only studied the coupling coefficients to zeroth order in the stellar angular velocity. However, it may be that the dominant energy-transfer channels involve coupling coefficients that are nonzero only at $`O(\mathrm{\Omega }^2)`$.
* We have focused attention in this paper on zero-buoyancy stars, for simplicity. However the quantitative details of the saturation process in real neutron stars may be altered by the presence of buoyancy forces.
* Finally, we have used Newtonian gravity throughout. It is conceivable that important energy transfer channels could come about via coupling coefficients which vanish at Newtonian order (to zeroth order in $`\mathrm{\Omega }`$), but which are non-vanishing when post-1-Newtonian corrections are included. This will depend on the relative sizes of the two dimensionless parameters $`GM/(c^2R)`$ and $`R^3\mathrm{\Omega }^2/(GM)`$. The fact that the gravitational radiation reaction force on the $`l=m=2`$ $`r`$-mode is dominated not by the “Newtonian” quadrupole coupling but instead by the “post-Newtonian” gravitomagnetic coupling is a hint in this direction.

###### Acknowledgements.

This work was supported in part by NSF grants PHY-9900672 and PHY-9722189 and NASA grants NAG5-7264 and NAG5-8356 to Cornell University. PA received support from the Jeffrey L. Bishop Fellowship and EF from the Alfred P. Sloan foundation. PA wishes to thank Yanqin Wu and Chris Matzner for many useful conversations. We thank Larry Kidder, Dong Lai, Sharon Morsink and Mark Scheel for useful conversations, and Sharon Morsink for detailed and helpful suggestions on the manuscript.

## A Hamiltonian Analysis of Linearized Perturbations

The purpose of this appendix is to derive the mode decomposition formalism described in Sec. II of the body of the paper. As discussed in Sec. II, most of that formalism is contained in a series of papers by Schutz and collaborators . Schutz introduces the phase space mode expansion, shows that the set of vectors (34) form a basis for $`\times `$, and gives an extensive discussion of Jordan chains. The main new feature that we introduce is the explicit computation of left eigenvectors \[defined in Eq. (A43) below\] in terms of right eigenvectors, and the subsequent derivation of the equations of motion in the explicit form (46) above. For completeness, we sketch in this Appendix derivations of all the building blocks of the formalism. Secs. A 1A 4 detail the formalism for non-Jordan chain modes. However, evolving the differential rotation and/or total spin of the star requires an extension of the formalism to compute the equations of motion for Jordan chain modes. That extension is given in Sec. A 5.

### 1 Phase space equations of motion

We start by reformulating the equation of motion (19) as a pair of first order equations rather than a single second order equation. Our treatment closely follows , except that we use canonically conjugate variables instead of $`𝝃`$ and $`\dot{𝝃}`$.

The equation of motion (19) can be derived from the Lagrangian density

$$=\frac{1}{2}\dot{𝝃}\dot{𝝃}+\frac{1}{2}\dot{𝝃}𝐁𝝃\frac{1}{2}𝝃𝐂𝝃+𝐚_{\mathrm{ext}}(t)𝝃.$$
(A1)
The momentum canonically conjugate to $`𝝃`$ is

$$𝝅=\frac{}{\dot{𝝃}}=\dot{𝝃}+\frac{1}{2}𝐁𝝃,$$
(A2)
and the associated Hamiltonian density is

$$=\frac{1}{2}\left(𝝅\frac{1}{2}𝐁𝝃\right)^2+\frac{1}{2}𝝃𝐂𝝃𝐚_{\mathrm{ext}}(t)𝝃.$$
(A3)
The Hamiltonian equations of motion can be written as

$`\dot{𝜻}=𝐓𝜻+𝐅(t),`$ (A4)
where

$`𝜻(t,𝐱)\left[\begin{array}{c}𝝃(t,𝐱)\\ 𝝅(t,𝐱)\end{array}\right],`$ (A7)
the operator $`𝐓`$ is

$`𝐓=\left[\begin{array}{cc}\frac{1}{2}𝐁& \mathrm{𝟏}\\ 𝐂+\frac{1}{4}𝐁^2& \frac{1}{2}𝐁\end{array}\right],`$ (A10)
and where

$`𝐅(t)=\left[\begin{array}{c}0\\ 𝐚_{\mathrm{ext}}(t)\end{array}\right].`$ (A13)
If we now specialize to the case of no forcing term, $`𝐚_{\mathrm{ext}}(t)=0`$, and assume a solution of the form

$`𝜻(t,𝐱)=e^{i\omega t}𝜻(𝐱),`$ (A14)
then we get from Eq. (A4) the eigenvalue equation

$`\left[𝐓+i\omega \right]𝜻(𝐱)=0.`$ (A15)
It is straightforward to show using Eqs. (A2) and (A7) that Eq. (A4) is equivalent to the configuration space equation of motion (19), and that the eigenvalue equation (A15) is equivalent to the quadratic eigenvalue equation (19).

### 2 Right and left eigenvectors and Jordan chains

Label the distinct right eigenvectors of $`𝐓`$ as $`𝜻_A`$, and the associated eigenfrequencies as $`\omega _A`$, so that

$`\left[𝐓+i\omega _A\right]𝜻_A=0.`$ (A16)
Since the operator $`𝐓`$ is not Hermitian, we will have also left eigenvectors $`𝝌_A`$ (distinct from the right eigenvectors) that satisfy

$`\left[𝐓^{}i\omega _A^{}\right]𝝌_A=0.`$ (A17)
Here

$`𝐓^{}=\left[\begin{array}{cc}\frac{1}{2}𝐁& 𝐂+\frac{1}{4}𝐁^2\\ \mathrm{𝟏}& \frac{1}{2}𝐁\end{array}\right],`$ (A20)
is the Hermitian conjugate of $`𝐓`$, since $`𝐂^{}=𝐂`$ and $`𝐁^{}=𝐁`$.

Since the operator $`𝐓`$ is not Hermitian, the set of its right eigenvectors will not in general be a complete basis. However, one can obtain a complete basis if one includes all Jordan chains . What this means is as follows. For a given eigenvalue $`i\omega _A`$, let $`𝒱_A`$ be the subspace of $`_2`$ consisting of vectors $`𝜻`$ that satisfy

$`\left[𝐓+i\omega _A\right]^m𝜻=0`$ (A21)
for some integer $`m1`$. Clearly the space $`𝒱_A`$ contains all the right eigenvectors associated with $`\omega _A`$. Now it can be shown that the direct sum of all the subspaces $`𝒱_A`$ (one for each distinct eigenvalue) gives the entire space $`_2`$. Hence, to obtain a basis for $`_2`$, it suffices to find a basis for each space $`𝒱_A`$. One can always find such a basis of the form $`\{𝜻_{A,\sigma }\}`$ where $`0\sigma p_A`$, $`𝜻_{A,0}=𝜻_A`$ is the right eigenvector, $`𝜻_{A,1},\mathrm{},𝜻_{A,p_A}`$ are a set vectors in $`𝒱_A`$ that satisfy

$`\left[𝐓+i\omega _A\right]𝜻_{A,\sigma }=𝜻_{A,\sigma 1},`$ (A22)
for $`1\sigma p_A`$, and that form the Jordan chain of length $`p_A`$ associated with $`𝜻_A`$. For each eigenvector $`𝜻_A`$, either $`p_A=0`$ and there is no associated Jordan chain, or there is a chain of length $`p_A1`$. The set $`\left\{𝜻_{A,\sigma }\right|A=1,2,3,\mathrm{},0\sigma p_A\}`$ of Jordan chains, including the right eigenvectors, forms a basis of $`_2`$.

These Jordan chains are right Jordan chains. We also have for each $`A`$ a left eigenvector $`𝝌_A=𝝌_{A,0}`$ as discussed above, and left Jordan chains consisting of vectors $`𝝌_{A,\sigma }`$ for $`1\sigma p_A`$ satisfying

$`\left[𝐓^{}i\omega _A^{}\right]𝝌_{A,\sigma }=𝝌_{A,\sigma 1},`$ (A23)
for $`1Ap_A`$. Note that the length $`p_A`$ of the left Jordan chain must be the same as that of the right Jordan chain, for each $`A`$. The basis $`\{𝝌_{A,\sigma }\}`$ of left Jordan chains can be chosen to be dual to the basis of right Jordan chains in the sense that

$`𝝌_{A,\sigma },𝜻_{B,\lambda }=\delta _{AB}\delta _{p_A,\sigma +\lambda },`$ (A24)
for all $`A,B`$ and for $`0\sigma p_A`$ and $`0\lambda p_B`$ . The inner product here is defined in the obvious way as

$$[𝝃,𝝅],[𝝃^{},𝝅^{}]𝝃,𝝃^{}+𝝅,𝝅^{},$$
(A25)
where the inner products on the right hand side are given by the definition (21).

### 3 General mode expansion and equations of motion for mode coefficients

We expand $`𝜻(t,𝐱)`$ on the basis of right Jordan chains as

$`𝜻(t,𝐱)={\displaystyle \underset{A}{}}{\displaystyle \underset{\sigma =0}{\overset{p_A}{}}}c_{A,\sigma }(t)𝜻_{A,\sigma }(𝐱),`$ (A26)
where from the orthogonality relation (A24) we have

$$c_{A,\sigma }=𝝌_{A,p_A\sigma },𝜻.$$
(A27)
Substituting the expansion (A26) into the dynamical equation (A4), and using the defining properties (A16) and (A22) of right Jordan chains yields

$`{\displaystyle \underset{A}{}}{\displaystyle \underset{\sigma =0}{\overset{p_A}{}}}\dot{c}_{A,\sigma }(t)𝜻_{A,\sigma }`$ $`=`$ $`{\displaystyle \underset{A}{}}{\displaystyle \underset{\sigma =0}{\overset{p_A}{}}}c_{A,\sigma }(t)\left[i\omega _A𝜻_{A,\sigma }+𝜻_{A,\sigma 1}\right]+𝐅(t),`$ (A28)
where we define $`𝜻_{A,1}0`$. Now multiplying on the left by $`𝝌_{B,\lambda }`$, using the orthogonality relation (A24), and then relabeling the indices yields the equations

$`\dot{c}_{A,\sigma }+i\omega _Ac_{A,\sigma }c_{A,\sigma +1}=𝝌_{A,p_A\sigma },𝐅(t),`$ (A29)
for $`0\sigma p_A1`$, and

$`\dot{c}_{A,p_A}+i\omega _Ac_{A,p_A}=𝝌_{A,0},𝐅(t).`$ (A30)
Thus, the equations for the mode coefficients $`c_{A,0}(t)`$, $`c_{A,1}(t)`$, $`\mathrm{}`$, $`c_{A,p_A}(t)`$ are coupled together for each fixed $`A`$, but the different $`A`$’s are uncoupled. The general solution of Eqs. (A29) and (A30) in the case of no forcing terms is a polynomial in time multiplied by the usual complex exponential:

$`c_{A,\sigma }(t)=e^{i\omega _At}{\displaystyle \underset{\lambda =0}{\overset{p_A\sigma }{}}}\gamma _{A,\sigma +\lambda }t^\lambda ,`$ (A31)
where $`\gamma _{A,0},\gamma _{A,1},\mathrm{},\gamma _{A,p_A}`$ are constants of integration.

### 4 Specialization to no Jordan chains

We now specialize to the situation where there are no Jordan chains. In this case the orthonormality relation (A24) reduces to

$`𝝌_A,𝜻_B=\delta _{AB},`$ (A32)
and the mode expansion (A26) becomes

$`𝜻(t,𝐱)={\displaystyle \underset{A}{}}c_A(t)𝜻_A(𝐱).`$ (A33)
The inverse of this mode expansion is

$$c_A(t)=𝝌_A,𝜻(t),$$
(A34)
and the equation of motions (A29) and (A30) for the mode coefficients reduce to

$`\dot{c}_A+i\omega _Ac_A=𝝌_A,𝐅(t).`$ (A35)
#### a Translation of results to configuration space variables

Our goal now is to write these results for no Jordan chains entirely in terms of the Lagrangian displacement $`𝝃`$ and its time derivative $`\dot{𝝃}`$, and the modes $`(𝝃_A,\omega _A)`$ of the configuration space formalism. To do this we proceed as follows. We write the right eigenvector $`𝜻_A`$ as

$`𝜻_A=\left[\begin{array}{c}𝝃_A\\ 𝝅_A\end{array}\right],`$ (A38)
and similarly write the left eigenvector $`𝝌_A`$ as

$`𝝌_A=\left[\begin{array}{c}𝝈_A\\ 𝝉_A\end{array}\right].`$ (A41)
Then the definitions (A16) –(A17) of right and left eigenvectors together with the formula (A10) for the operator $`𝐓`$ shows that $`𝝃_A`$ and $`𝝉_A`$ are right and left eigenvectors of the operator $`𝐋(\omega _A)=\omega _A^2i\omega _A𝐁+𝐂`$ \[cf. Eqs. (29) and (31) above\]

$`𝐋(\omega _A)𝝃_A`$ $`=`$ $`0`$ (A42)
$`𝐋(\omega _A)^{}𝝉_A`$ $`=`$ $`0.`$ (A43)
In addition we obtain the formulae

$$𝝅_A=i\omega _A𝝃_A+𝐁𝝃_A/2$$
(A44)
and

$$𝝈_A=i\omega _A^{}𝝉_A𝐁𝝉_A/2,$$
(A45)
which using Eqs. (A48) and (A51) allow us to write the phase space right and left eigenvectors $`𝜻_A`$ and $`𝝌_A`$ entirely in terms of the configuration space left and right eigenvectors $`𝝃_A`$ and $`𝝉_A`$:

$`𝜻_A=\left[\begin{array}{c}𝝃_A\\ i\omega _A𝝃_A+\frac{1}{2}𝐁𝝃_A\end{array}\right]`$ (A48)
and

$`𝝌_A=\left[\begin{array}{c}i\omega _A^{}𝝉_A\frac{1}{2}𝐁𝝉_A\\ 𝝉_A\end{array}\right].`$ (A51)
Now using the definitions (A2) and (A7) and the formula (A48) we can rewrite the mode expansion (A33) as

$`\left[\begin{array}{c}𝝃(t)\\ \dot{𝝃}(t)+\frac{1}{2}𝐁𝝃(t)\end{array}\right]={\displaystyle \underset{A}{}}c_A(t)\left[\begin{array}{c}𝝃_A\\ i\omega _A𝝃_A+\frac{1}{2}𝐁𝝃_A\end{array}\right],`$ (A56)
which is equivalent to the expansion (43) quoted in the body of the paper. Note that there is a one-to-one correspondence between right eigenvectors solutions $`𝜻_A`$ of Eq. (A16) and solutions $`(𝝃_A,\omega _A)`$ of Eq. (19), which allows us to identify the sums over $`A`$ that appear in the expansions (43) and (A56).

Next, we can use the formulae (A48) and (A51) to rewrite the orthogonality relation (A32), the inverse mode expansion (A34) and the equation of motion (A35). The results are

$$\delta _{AB}=𝝉_A,𝐁𝝃_Bi(\omega _A+\omega _B)𝝉_A,𝝃_B,$$
(A57)
$$c_A(t)=i𝝉_A,\omega _A𝝃(t)+i𝐁𝝃(t)+i\dot{𝝃}(t).$$
(A58)
and

$`\dot{c}_A+i\omega _Ac_A=𝝉_A,𝐚_{\mathrm{ext}}(t),`$ (A59)
where we have also used Eqs. (23), (A2), (A7) and (A13).

#### b Computation of the left eigenmodes

The last step in the construction is the explicit computation of the left eigenvectors $`𝝉_A`$ in terms of the right eigenvectors $`𝝃_A`$. For a given eigenfrequency $`\omega `$, let $`𝒱_R(\omega )`$ denote the space of associated right eigenvectors $`𝝃`$ satisfying $`𝐋(\omega )𝝃=0`$ \[cf. Eq. (29)\], and let $`𝒱_L(\omega )`$ denote the space of associated left eigenvectors $`𝝉`$ satisfying $`𝐋(\omega )^{}𝝉=0`$. The spaces $`𝒱_R(\omega _A)`$ and $`𝒱_L(\omega _A)`$ have the same dimension, which is the degeneracy associated with the eigenfrequency $`\omega _A`$. However, $`𝒱_R(\omega _A)`$ and $`𝒱_L(\omega _A)`$ will not in general coincide for complex $`\omega _A`$. We now show that for real $`\omega _A`$, the spaces $`𝒱_L(\omega _A)`$ and $`𝒱_R(\omega _A)`$ do coincide, which means that the left eigenvectors associated with a given eigenfrequency can be expressed as linear combinations of the corresponding right eigenvectors. This is the key result that we use to compute the left eigenvectors.

We shall need the following properties of the spaces $`𝒱_R(\omega )`$ and $`𝒱_L(\omega )`$ of right and left eigenvectors discussed by Schutz . Suppose that $`𝝃`$ lies in $`𝒱_R(\omega )`$. Then

* The function $`𝝃`$ also lies in $`𝒱_L(\omega ^{})`$, ie, is a left eigenvector with eigenfrequency $`\omega ^{}`$. This can be derived by taking the Hermitian conjugate of Eq. (19), and by using the fact that $`𝐁`$ is anti-Hermitian and $`𝐂`$ is Hermitian, which yields $`𝐋(\omega ^{})^{}𝝃=0`$.
* The function $`𝝃(𝐱)^{}`$ lies in $`𝒱_R(\omega ^{})`$. This follows directly from the fact that the equation of motion (7) is real.
* Let $`g`$ be the mapping that takes $`(r,\theta ,\phi )`$ to $`(r,\theta ,\phi )`$ in spherical polar coordinates, or equivalently $`(x,y,z)(x,y,z)`$ in Cartesian coordinates, i.e., reflection in the $`xz`$ plane. We define $`g_{}𝝃`$ be the pullback of $`𝝃`$ under $`g`$ , so that if

$`𝝃=\xi ^r(r,\theta ,\phi ){\displaystyle \frac{}{r}}+\xi ^\theta (r,\theta ,\phi ){\displaystyle \frac{}{\theta }}+\xi ^\phi (r,\theta ,\phi ){\displaystyle \frac{}{\phi }},`$ (A60)
then

$`g_{}𝝃`$ $`=`$ $`\xi ^r(r,\theta ,\phi ){\displaystyle \frac{}{r}}+\xi ^\theta (r,\theta ,\phi ){\displaystyle \frac{}{\theta }}\xi ^\phi (r,\theta ,\phi ){\displaystyle \frac{}{\phi }}.`$ (A61)
Since the background star is invariant under the transformation $`\phi \phi `$, $`tt`$, it follows that $`g_{}𝝃`$ lies in $`𝒱_R(\omega )`$ .
* By combining the previous properties, it follows that $`(g_{}𝝃)^{}`$ belongs to both $`𝒱_R(\omega ^{})`$ and $`𝒱_L(\omega )`$. The mapping $`𝝃(g_{}𝝃)^{}`$ was first written down by Schutz in Ref.  where it was denoted $`S`$.

We next recall the notation used in Sec. II C 3 above. We write the distinct eigenfrequencies as $`\omega _a`$, and the right eigenvectors as $`𝝃_A=𝝃_{a,k}`$ where $`1kn_a`$ and $`n_a`$ is the degeneracy associated with the eigenfrequency $`\omega _a`$.

Consider first the case of complex $`\omega _a`$. In this case it follows from the above that the mapping $`𝝃(g_{}𝝃)^{}`$ maps $`𝒱_R(\omega _a)`$ onto $`𝒱_L(\omega _a)`$ . Hence we can write

$$𝝉_{a,k}=\underset{k^{}=1}{\overset{n_a}{}}T_{kk^{}}(g_{}𝝃_{a,k^{}})^{},$$
(A62)
for some $`n_a\times n_a`$ matrix $`T_{kk^{}}`$. By substituting Eq. (A62) into the orthogonality relation (A57) we can solve for the matrix $`T_{kk^{}}`$ and thus obtain the left eigenvectors $`𝝉_{a,k}`$. We note that Schutz uses a different basis of left eigenmodes for which $`𝝉_{a,k}=(g_{}𝝃_{a,k})^{}`$ \[Eq. (3.14) of \], but for which the orthonormality condition (A57) does not hold. This method of obtaining the left eigenmodes for complex frequencies is discussed in more detail in Sec. A 5 below.

For real frequencies, the situation is much simpler. Every $`𝝃`$ in $`𝒱_R(\omega )`$ also lies in $`𝒱_L(\omega ^{})=𝒱_L(\omega )`$, and hence the spaces $`𝒱_R(\omega )`$ and $`𝒱_L(\omega )`$ coincide †¶†¶†¶For a mode $`(𝝃_A,\omega _A)`$ with no degeneracy, it follows from this that we can choose the phase of $`𝝃_A`$ so that $`(g_{}𝝃_A)^{}=𝝃_A`$. Hence, if $`m`$ is the azimuthal quantum number and we write $`𝝃_A=\mathrm{exp}[im\phi ]\widehat{𝝃}_A`$, the $`r`$ and $`\theta `$ components of $`\widehat{𝝃}_A`$ are purely real and the $`\phi `$ component is purely imaginary.. Hence we can write

$$𝝉_{a,k}=\underset{k^{}=1}{\overset{n_a}{}}i_{kk^{}}^{(a)}𝝃_{a,k^{}},$$
(A63)
for some $`n_a\times n_a`$ matrix $`_{kk^{}}^{(a)}`$. (The factor of $`i`$ and the complex conjugation are included for convenience.) By substituting Eq. (A63) with $`𝝉_A=𝝉_{a,k}`$ and $`𝝃_B=𝝃_{b,l}`$ into the orthogonality relation (A57) we obtain

$$\delta _{ab}\delta _{kl}=_{kk^{}}^{(a)}_{ak^{},bl},$$
(A64)
where the matrix $`_{ak^{},bl}`$ is defined in Eq. (50). It follows from Eq. (A64) that, first, the matrix $`_{ak,bl}`$ is block diagonal in the sense claimed in Sec. II C 3, and second, for a given value of $`a`$ the matrix $`_{kk^{}}^{(a)}`$ is just the inverse of the the diagonal block $`_{ak,al}`$.

Now as explained in Sec. II C 3, we can always choose the basis $`𝝃_{a,k}`$ of $`𝒱_R(\omega _a)`$ to diagonalize the matrix $`_{ak^{},al}`$. In particular, this will be automatically true when $`n_a=1`$ and there is no degeneracy. For such bases, the matrix $`B_{kk^{}}`$ will also be diagonal. It then follows from Eq. (A63) that each left eigenvector can be written as

$$𝝉_A=\frac{i}{b_A}𝝃_A$$
(A65)
for some constant $`b_A`$. Substituting this into the orthonormality relation (A57) yields the formula (45) for $`b_A`$. Finally, the formula (A65) for left eigenvectors can be combined with the orthogonality relation (A57), the inverse mode expansion (A58) and equation of motion (A59) to yield the versions (35), (44) and (46) of these relations quoted in the body of the paper †∥†∥†∥The validity of the orthogonality relation (35), for those values of $`A`$ and $`B`$ for which $`\omega _A\omega _B`$, can also be derived directly by contracting Eq. (31) on the left with $`𝝃_B`$, subtracting from this the same equation complex conjugated with $`A`$ and $`B`$ interchanged, and dividing by $`\omega _A\omega _B`$..

Finally, we note that Eq. (A64) shows that the matrix $`_{ak,al}`$ is non-degenerate. Thus, if that matrix is degenerate, then the assumption underlying the derivation of Eq. (A64) — that there is no Jordan chain associated with the eigenfrequency $`\omega _a`$ — must fail. In particular, for non-degenerate modes, if the constant $`b_A`$ defined by Eq. (45) vanishes, then the mode must be a Jordan chain mode.

### 5 Computation of the left eigenmodes for Jordan chain modes

In this subsection we show how to generalize the above analysis to obtain the left eigenvectors $`𝝌_{A,\sigma }`$ for Jordan chain modes. The construction is useful, for example, in computing the evolution of stellar differential rotation, which is described by an infinite set of Jordan chains of length 1 (see Appendix D).

The essential idea is to use the phase space variables rather than the more complicated configuration space variables, and to use an operator which maps right eigenvectors onto left eigenvectors. We define the operator

$`𝐌=i\left[\begin{array}{cc}0& \mathrm{𝟏}\\ \mathrm{𝟏}& 0\end{array}\right],`$ (A68)
which is essentially the symplectic structure defined by Friedman and Schutz . It satisfies $`𝐌^{}=𝐌=𝐌^1`$ and

$$𝐌^{}𝐓𝐌=𝐓^{}.$$
(A69)
It follows from Eqs. (A69), (A16) and (A17) that if $`𝜻`$ is a right eigenvector of $`𝐓`$ with eigenfrequency $`\omega `$, then $`𝐌𝜻`$ is a left eigenvector with eigenfrequency $`\omega ^{}`$.

Now the set of vectors $`𝐌𝜻_{A,\sigma }`$ forms a basis, since $`𝐌`$ is invertible and the set of vectors $`𝜻_{A,\sigma }`$ is a basis. Therefore we can write the left eigenvectors and left Jordan chain vectors as

$$𝝌_{A,\sigma }=\underset{B,\lambda }{}_{A\sigma ,B\lambda }^{}𝐌𝜻_{B,\lambda },$$
(A70)
for some matrix $`_{A\sigma ,B\tau }`$, where the complex conjugation is included for later convenience. Inserting this into the orthogonality relation (A24) gives

$$_{A\sigma ,C\tau }_{C\tau ,B\lambda }=\delta _{AB}\delta _{\sigma +\lambda ,p_A},$$
(A71)
where

$$_{C\tau ,B\lambda }=𝜻_{C,\tau },𝐌𝜻_{B,\lambda }$$
(A72)
are the matrix elements of the operator $`𝐌`$. Equation (A71) says that the matrices $``$ and $``$ are inverses of each other (up to the index permutation $`\sigma p_A\sigma `$).

The reason the ansatz (A70) is useful is that the matrix $``$ is almost diagonal, in the sense that it satisfies the identity

$$(\omega _A^{}\omega _B)_{A\sigma ,B\lambda }=0.$$
(A73)
It follows from Eq. (A73) that $``$ is block diagonal, with one block for each real frequency and one block for each pair $`(\omega ,\omega ^{})`$ of complex frequencies. Therefore, we can obtain the left eigenvectors by inverting each diagonal block of $`_{A\sigma ,B\lambda }`$ to obtain $`_{A\sigma ,B\lambda }`$, and by using Eq. (A70).

To derive the identity (A73) we use Eqs. (A69) and (A22) to show that the vectors $`𝐌𝜻_{A,\sigma }`$ form a left Jordan chain with frequency $`\omega _A^{}`$:

$$\left[𝐓^{}i\omega _A\right](𝐌𝜻_{A,\sigma })=𝐌𝜻_{A,\sigma 1}.$$
(A74)
Next, we have

$`\omega _A^{}_{A\sigma ,B\lambda }`$ $`=`$ $`\omega _A^{}𝜻_{A,\sigma },𝐌𝜻_{B,\lambda }`$ (A75)
$`=`$ $`\omega _A𝐌𝜻_{A,\sigma },𝜻_{B,\lambda }`$ (A76)
$`=`$ $`i𝐌𝜻_{A,\sigma 1}i𝐓^{}𝐌𝜻_{A,\sigma },𝜻_{B,\lambda }`$ (A77)
$`=`$ $`i𝜻_{A,\sigma 1},𝐌𝜻_{B,\lambda }+i𝐌𝜻_{A,\sigma },𝐓𝜻_{B,\lambda }`$ (A78)
$`=`$ $`i𝜻_{A,\sigma 1},𝐌𝜻_{B,\lambda }+i𝐌𝜻_{A,\sigma },i\omega _B𝜻_{B,\lambda }+𝜻_{B,\lambda 1},`$ (A79)
where we have used Eqs. (A22) and (A74). The result (A79) can be written as

$$(\omega _A^{}\omega _B)_{A\sigma ,B\lambda }=i\left(_{A(\sigma 1),B\lambda }+_{A\sigma ,B(\lambda 1)}\right),$$
(A80)
and the identity (A73) follows from iterating the identity (A80).

Schutz  suggested a different general method of constructing left Jordan chains. That method is based on the anti-linear operator $`𝐒_2`$ defined by

$`𝐒_2\left[\begin{array}{c}𝝃\\ 𝝅\end{array}\right]=\left[\begin{array}{c}(g_{}𝝃)^{}\\ (g_{}𝝅)^{}\end{array}\right],`$ (A85)
which maps right Jordan chains $`𝜻_{A,\sigma }`$ of $`𝐓`$ with frequency $`\omega _A`$ onto right Jordan chains $`𝐒_2𝜻_{A,\sigma }`$ with frequency $`\omega _A^{}`$. Therefore composing this map with the operator $`𝐌`$ yields a mapping $`𝜻𝐌𝐒_2𝜻`$ which takes right Jordan chains of frequency $`\omega `$ to left Jordan chains of frequency $`\omega `$. The approach here omits the mapping $`𝐒_2`$ and is simpler to use when the eigenfrequencies are real.

For Jordan chains of length zero (ordinary eigenvectors), the method (A70) reduces to the method (A63) of Sec. A 4 above. Dropping the Jordan chain indices $`\sigma `$ and $`\lambda `$, the ansatz (A70) reduces to $`𝝌_A=_B_{AB}^{}𝐌𝜻_B`$ with $`_{AC}_{CB}=\delta _{AB}`$, which using Eqs. (A48), (A51) and (A68) implies

$$𝝉_A=i\underset{B}{}_{AB}^{}𝝃_B,$$
(A86)
cf. Eq. (A63) above. Similarly, the definition (A72) can be simplified using Eqs. (A48) and (A68) to give

$`_{AB}`$ $`=`$ $`𝜻_A,𝐌𝜻_B`$ (A87)
$`=`$ $`𝝃_A,i𝐁𝝃_B+(\omega _A^{}+\omega _B)𝝃_A,𝝃_B,`$ (A88)
which is a generalization of the definition (50) valid for complex frequencies.

As an example, consider now the case of a mode $`(𝝃_1,\omega _1)`$ with complex eigenfrequency which is non-degenerate. Then there is an associated mode $`(𝝃_2,\omega _2)`$ with $`\omega _2=\omega _1^{}`$, and by choosing the normalization of $`𝝃_2`$ we can without loss of generality take

$$𝝃_2=(g_{}𝝃_1)^{},$$
(A89)
cf. Sec. A 4 b above. The corresponding $`2\times 2`$ block of the matrix $``$ is from Eqs. (A73) and (A88) of the form

$`_{AB}=\left[\begin{array}{cc}0& 𝒟\\ 𝒟^{}& 0\end{array}\right],`$ (A92)
where $`𝒟=_{12}=N[𝝃_1,\omega _1]`$, where the functional $`N`$ is

$$N[𝝃,\omega ]𝝃,i𝐁(g_{}𝝃)^{}+2\omega ^{}𝝃,(g_{}𝝃)^{}.$$
(A93)
The corresponding left eigenmodes are therefore

$$𝝉_1=\frac{i}{𝒟}𝝃_2=\frac{i}{𝒟}(g_{}𝝃_1)^{}$$
(A94)
and

$$𝝉_2=\frac{i}{𝒟^{}}𝝃_1.$$
(A95)
#### a Zero frequency Jordan chains of length one

We now carry through the explicit computation of the left Jordan chains for the case that arises in practice in stable rotating stars, that of zero-frequency Jordan chains of length one. Consider first the case of no degeneracy. Dropping the index $`A`$, we have from Eq. (A22) two Jordan chain vectors $`𝜻_0`$ and $`𝜻_1`$ which satisfy

$$𝐓𝜻_0=0,𝐓𝜻_1=𝜻_0.$$
(A96)
The $`2\times 2`$ matrix $`_{\sigma \tau }=𝜻_\sigma ,𝐌𝜻_\tau `$ satisfies the identity

$$_{(\sigma 1)\tau }+_{\sigma (\tau 1)}=0$$
(A97)
from Eq. (A80). Using this identity together with the fact that $``$ is Hermitian gives

$`_{\sigma \tau }=\left[\begin{array}{cc}0& i\beta \\ i\beta & \gamma \end{array}\right],`$ (A100)
where $`\beta `$ and $`\gamma `$ are real. Under the transformation $`𝜻_1𝜻_1+\chi 𝜻_0`$, which preserves the defining relations (A96), we have

$$\beta \beta ,\gamma \gamma +2\beta \mathrm{Im}(\chi ).$$
(A101)
It follows that we can choose $`𝜻_1`$ to make $`\gamma =0`$. Now combining Eqs. (A70), (A71) and (A100) gives for the left Jordan chain vectors

$`𝝌_0`$ $`=`$ $`{\displaystyle \frac{i}{\beta }}𝐌𝜻_0`$ (A102)
$`𝝌_1`$ $`=`$ $`{\displaystyle \frac{i}{\beta }}𝐌𝜻_1.`$ (A103)
We now write these results in terms of configuration space variables, using the notation

$`𝜻_\sigma =\left[\begin{array}{c}𝝃_\sigma \\ 𝝅_\sigma \end{array}\right],`$ (A106)
$`𝝌_\sigma =\left[\begin{array}{c}𝝈_\sigma \\ 𝝉_\sigma \end{array}\right],`$ (A109)
for $`\sigma =0,1`$. First, the relations (A96) can be written as

$$𝐂𝝃_0=0,𝐂𝝃_1=𝐁𝝃_0,𝝅_0=𝐁𝝃_0/2,𝝅_1=𝝃_0+𝐁𝝃_1/2.$$
(A110)
Second, the formula for $`\beta `$ is, from Eqs. (A68), (A106) and (A110)

$`\beta `$ $`=`$ $`i𝜻_0,𝐌𝜻_1`$ (A111)
$`=`$ $`𝝃_0,𝝃_0+𝝃_0,𝐁𝝃_1.`$ (A112)
Third, the relations (A102) and (A103) together with Eqs. (A68), (A106) and (A109) imply that

$$𝝉_0=\frac{1}{\beta }𝝃_0,𝝉_1=\frac{1}{\beta }𝝃_1.$$
(A113)
The corresponding equations of motion are, from Eqs. (A29) and (A30),

$`\dot{c}_1`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}𝝃_0,𝐚_{\mathrm{ext}}`$ (A114)
$`\dot{c}_0`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}𝝃_1,𝐚_{\mathrm{ext}}+c_1.`$ (A115)
We now generalize these results to allow for degeneracy, since the space of zero-frequency Jordan chain modes in rotating stars is highly degenerate (see Appendix D). We use the notation

$$A=(a,k)$$
(A116)
of Sec. II C 3 above, where $`a`$ labels the distinct eigenfrequencies and $`k`$ the eigenvectors associated with each eigenfrequency. Therefore the Jordan chain vectors can be written as

$$𝜻_{A\sigma }=𝜻_{ak\sigma }$$
(A117)
where $`\sigma `$ labels the Jordan chain vectors associated with each eigenvector. For simplicity we drop the label $`a`$ in what follows. The defining relations for Jordan chains of length one are

$$𝐓𝜻_{k0}=0,𝐓𝜻_{k1}=𝜻_{k0}.$$
(A118)
The matrix $`_{k\sigma ,l\tau }=𝜻_{k\sigma },𝐌𝜻_{l\tau }`$ satisfies the identity

$$_{k(\sigma 1),l\tau }+_{k\sigma ,l(\tau 1)}=0$$
(A119)
from Eq. (A80). Using this identity together with the fact that $``$ is Hermitian gives

$`_{k0,l0}`$ $`=`$ $`0`$ (A120)
$`_{k0,l1}`$ $`=`$ $`_{k1,l0}=i\beta _{kl}`$ (A121)
$`_{k1,l1}`$ $`=`$ $`\gamma _{kl},`$ (A122)
where the matrices $`\beta _{kl}`$ and $`\gamma _{kl}`$ are Hermitian.

Next, the transformation

$`𝜻_{k0}`$ $``$ $`𝜻_{k0}^{}=F_{kl}^{}𝜻_{l0}`$ (A123)
$`𝜻_{k1}`$ $``$ $`𝜻_{k1}^{}=F_{kl}^{}𝜻_{l1}`$ (A124)
preserves the relations (A118). Under this transformation the matrix $`\beta _{kl}`$ transforms as $`𝜷𝜷^{}=𝐅𝜷𝐅^{}`$, and it follows that we can choose the basis $`𝜻_{k0}`$ to diagonalize $`\beta _{kl}`$, so that

$$\beta _{kl}=\beta _k\delta _{kl}.$$
(A125)
Similarly the transformation

$`𝜻_{k0}`$ $``$ $`𝜻_{k0}^{}=𝜻_{k0}`$ (A126)
$`𝜻_{k1}`$ $``$ $`𝜻_{k1}^{}=𝜻_{k1}+F_{kl}^{}𝜻_{l0}`$ (A127)
preserves the relations (A118). Under this transformation the matrices transform as

$`𝜷`$ $``$ $`𝜷`$ (A128)
$`𝜸`$ $``$ $`𝜸^{}=𝜸+i𝐅𝜷i𝜷𝐅^{}.`$ (A129)
Now it follows from Eqs. (A120)–(A122) and the fact that the matrix $`_{k\sigma ,l\tau }`$ is invertible that $`\beta _{kl}`$ is invertible. Therefore we can find a transformation matrix $`F_{kl}`$ which achieves $`𝜸^{}=0`$ in Eq. (A129), namely

$$𝐅=i𝜸𝜷^1/2.$$
(A130)
Therefore we can choose the basis $`𝜻_{k1}`$ to make $`\gamma _{kl}=0`$. The rest of the analysis now proceeds exactly as for the non-degenerate case above, and we obtain the equations of motion

$`\dot{c}_{k1}`$ $`=`$ $`{\displaystyle \frac{1}{\beta _k}}𝝃_{k0},𝐚_{\mathrm{ext}}`$ (A131)
$`\dot{c}_{k0}`$ $`=`$ $`{\displaystyle \frac{1}{\beta _k}}𝝃_{k1},𝐚_{\mathrm{ext}}+c_{k1}.`$ (A132)
## B Linear Independence of a Subset of the Right Eigenvectors

In this appendix we show that the of non-Jordan-chain right eigenmodes $`𝝃_A`$ for which the constant $`b_A`$ is positive are linearly independent. We start by recalling the notation used in Sec. II C 4 above: the index $`\alpha `$ labels the distinct right eigenvectors $`𝝃_A`$ for which $`b_A>0`$, and we write $`A=(\alpha ,ϵ)`$ where $`ϵ`$ takes on the values $`ϵ=+`$ and $`ϵ=`$. Then $`𝝃_\alpha 𝝃_{\alpha ,+}`$ and $`𝝃_{\alpha ,}=𝝃_\alpha ^{}`$, while $`\omega _{\alpha ,\pm }=\pm \omega _\alpha `$ and $`b_{\alpha ,\pm }=\pm b_\alpha `$.

By combining the orthogonality relation (A57) with the formula (A65) for left eigenvectors $`𝝉_A`$, and choosing $`A=(\alpha ,+)`$ and $`B=(\beta ,+)`$ we obtain the following form of the orthogonality relation:

$$b_\alpha \delta _{\alpha \beta }=𝝃_\alpha ,i𝐁𝝃_\beta +(\omega _\alpha +\omega _\beta )𝝃_\alpha ,𝝃_\beta .$$
(B1)
Suppose now that the vectors $`𝝃_\alpha `$ are not linearly independent. Then, one can find coefficients $`c_\alpha `$, not all zero, so the vector

$$𝐪\underset{\alpha }{}c_\alpha 𝝃_\alpha $$
(B2)
is zero. Now contract Eq. (B1) with $`c_\alpha ^{}c_\beta `$ and sum over $`\alpha `$ and $`\beta `$. The result is

$`{\displaystyle \underset{\alpha }{}}|c_\alpha |^2b_\alpha `$ $`=`$ $`𝐪,i𝐁𝐪+2𝐪,{\displaystyle \underset{\alpha }{}}c_\alpha \omega _\alpha 𝝃_\alpha `$ (B3)
$`=`$ $`0,`$ (B4)
where the second line follows from $`𝐪=0`$. Since on the left hand side $`b_\alpha >0`$ for all $`\alpha `$, this forces $`c_\alpha =0`$ for all $`\alpha `$, which contradicts our assumption above.

## C Alternative Form of Equations of Motion

In this appendix we present a version of the general equations of motion which is second order in time and which is more similar in form to the standard equations of motion for non-rotating stars. We start with the equations derived in Sec. II C 4 for the coefficients $`c_{\alpha ,+}(t)`$ and $`c_{\alpha ,}(t)`$:

$`\dot{c}_{\alpha ,+}+i\omega _\alpha c_{\alpha ,+}`$ $`=`$ $`{\displaystyle \frac{i}{b_\alpha }}f_{\alpha ,+}(t),`$ (C1)
$`\dot{c}_{\alpha ,}i\omega _\alpha c_{\alpha ,+}`$ $`=`$ $`{\displaystyle \frac{i}{b_\alpha }}f_{\alpha ,}(t),`$ (C2)
where $`f_{\alpha ,\pm }(t)𝝃_{\alpha ,\pm },𝐚_{\mathrm{ext}}(t)`$. As noted in Sec. II C 4, the Lagrangian displacement $`𝝃(𝐱,t)`$ will be real if and only if $`c_{\alpha ,+}=c_{\alpha ,}^{}`$ and $`f_{\alpha ,+}=f_{\alpha ,}^{}`$. Here, however, we allow arbitrary complex $`𝝃(𝐱,t)`$. We define

$`e_\alpha (t)`$ $``$ $`c_{\alpha ,+}(t)+c_{\alpha ,}(t)`$ (C3)
$`=`$ $`{\displaystyle \frac{2\omega _\alpha }{b_\alpha }}\mathrm{Re}(𝝃_\alpha ),𝝃(t)+{\displaystyle \frac{2}{b_\alpha }}\mathrm{Im}(𝝃_\alpha ),\dot{𝝃}(t)+𝐁𝝃(t),`$ (C4)
where we have used Eqs. (66) and (67). It is now straightforward to show using the definition (C3) that the equations of motion (C1)–(C2) are equivalent to

$$\ddot{e}_\alpha (t)+\omega _\alpha ^2e_\alpha (t)=_\alpha (t),$$
(C5)
which is the standard forced harmonic oscillator equation of motion. The forcing term here is

$`_\alpha `$ $`=`$ $`{\displaystyle \frac{i}{b_\alpha }}(\dot{f}_{\alpha ,}+i\omega _\alpha f_{\alpha ,})+{\displaystyle \frac{i}{b_\alpha }}(\dot{f}_{\alpha ,+}i\omega _\alpha f_{\alpha ,+})`$ (C6)
$`=`$ $`{\displaystyle \frac{2\omega _\alpha }{b_\alpha }}\mathrm{Re}(𝝃_\alpha ),𝐚_{\mathrm{ext}}+{\displaystyle \frac{2}{b_\alpha }}\mathrm{Im}(𝝃_\alpha ),\dot{𝐚}_{\mathrm{ext}}.`$ (C7)
It is clear that $`e_\alpha (t)`$ and $`_\alpha (t)`$ will be real if and only if $`𝐚_{\mathrm{ext}}(t)`$ and $`𝝃(𝐱,t)`$ are real. We can express the Lagrangian displacement $`𝝃(𝐱,t)`$ in terms of $`e_\alpha (t)`$ and $`\dot{e}_\alpha (t)`$ at any time $`t`$ by using the first component of the mode expansion (65), and by using the relations

$`c_{\alpha ,+}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(e_\alpha +{\displaystyle \frac{i}{\omega _\alpha }}\dot{e}_\alpha \right)+{\displaystyle \frac{1}{2\omega _\alpha b_\alpha }}\left(f_{\alpha ,+}f_{\alpha ,}\right),`$ (C8)
$`c_{\alpha ,}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(e_\alpha {\displaystyle \frac{i}{\omega _\alpha }}\dot{e}_\alpha \right){\displaystyle \frac{1}{2\omega _\alpha b_\alpha }}\left(f_{\alpha ,+}f_{\alpha ,}\right).`$ (C9)
In this formalism, the number of “modes” $`e_\alpha (t)`$ is the same for rotating stars as for non-rotating stars. In particular, for non-rotating stars the expansion coefficients $`e_\alpha `$ coincide with the standard expansion coefficients $`q_\alpha `$ given by Eq. (26), when the mode functions $`𝝃_\alpha (𝐱)`$ are chosen to be real.

Note that the equation of motion (C5) for a given mode has a characteristic feature which is peculiar to rotating stars, namely the forcing term depends on the time derivative $`\dot{𝐚}_{\mathrm{ext}}(𝐱,t)`$ of the externally applied force per unit mass, as well as on $`𝐚_{\mathrm{ext}}(𝐱,t)`$ itself. For non-rotating stars one can choose a real mode basis $`𝝃_\alpha `$ which removes the dependence on $`\dot{𝐚}_{\mathrm{ext}}`$ \[cf. the second term in Eq. (C7)\], but for rotating stars it is not possible to find real mode bases.

## D A Set of Zero-Frequency, Jordan-Chain Modes

In this appendix we justify the claim made in Sec. II C 2 above that there are always Jordan-chain modes present in rotating zero-buoyancy stars, even for stable stars. We show that the well known class of purely axial, zero frequency modes that move the star to nearby equilibrium states with different angular velocity are Jordan-chain modes. Note that in this appendix we do not need to assume that the angular velocity $`\mathrm{\Omega }`$ of the star is small. For non-rotating stars, the existence of these Jordan chains was previously noted in a footnote by Ref. .

It is most convenient to work directly with the linearized Eulerian equations of motion (5) and (6). Substitute into these equations the ansatz $`\delta \rho (𝐱,t)=\delta \rho (𝐱)`$ and

$$\delta 𝐮(𝐱,t)=\delta 𝐮(𝐱)=r_{}\delta \mathrm{\Omega }(r_{})𝐞_\phi ,$$
(D1)
where $`(r_{},\phi ,z)`$ are cylindrical coordinates with $`r_{}`$ the distance from the rotation axis, and $`𝐞_\phi `$ is the unit vector in the $`\phi `$ direction. One finds that this yields a time-independent solution, and that one can specify the perturbation $`\delta \mathrm{\Omega }(r_{})`$ arbitrarily and solve for the density perturbation $`\delta \rho (𝐱)`$ †\**†\**†\**The density perturbation always vanishes if the background star is non-rotating but not in general..

Now switch back to the language of Lagrangian perturbation theory. The equations for a Jordan chain of length one with $`\omega _A=0`$ are as follows. \[We drop the index $`A`$ and keep only the index $`\sigma `$ for convenience.\] The phase space right eigenvector is, from Eq. (A48),

$`𝜻_0=\left[\begin{array}{c}𝝃_0\\ \frac{1}{2}𝐁𝝃_0\end{array}\right],`$ (D4)
and the associated right Jordan chain vector is

$`𝜻_1=\left[\begin{array}{c}𝝃_1\\ \frac{1}{2}𝐁𝝃_1+𝝃_0\end{array}\right].`$ (D7)
Here $`𝝃_0`$ and $`𝝃_1`$ satisfy, from Eq. (A22),

$$𝐋(0)𝝃_0=𝐂𝝃_0=0,$$
(D8)
and

$$𝐋(0)𝝃_1=𝐂𝝃_1=𝐁𝝃_0.$$
(D9)
Now let $`\xi _0(𝐱)`$ be a mode function of the form (D1). Then one finds that $`𝝃_0=0`$ and $`(\rho 𝝃_0)=0`$, and hence from the definition (14)–(16) of the operator $`𝐂`$, Eq. (D8) is satisfied. One can then always solve Eq. (D9) to determine $`𝝃_1`$, since the left hand side is a pure gradient for zero-buoyancy stars and the curl of the right hand side vanishes by Eq. (D1). Thus the solutions (D1) correspond to a Jordan chains of length one. The vector $`𝝃_1`$ encodes the density perturbation. The general solution $`𝝃(𝐱,t)`$ associated with $`𝝃_0`$ is then, from Eq. (A31),

$$𝝃(𝐱,t)=\gamma _0𝝃_0(𝐱)+\gamma _1\left[𝝃_1(𝐱)+𝝃_0(𝐱)t\right],$$
(D10)
where $`\gamma _0`$ and $`\gamma _1`$ are constants of integration.

The piece of the solution proportional to $`\gamma _0`$ has vanishing Eulerian density and velocity perturbations, and is pure gauge in the sense of Ref. . By contrast, the piece proportional to $`\gamma _1`$ is a physical, non-gauge mode.

For non-rotating stars, the situation is a little different: all zero frequency modes are associated with Jordan chains of length one. This can be seen from the fact that Eq. (D9) always has solutions for non-rotating stars, since $`𝐁=0`$, whereas for rotating stars a solution only exists if the curl of $`𝐁𝝃_0`$ vanishes. The associated solutions of Eq. (7) are of the form $`𝝃(t)=𝝃_1+t𝝃_0`$ where $`𝝃_0`$ and $`𝝃_1`$ satisfy $`𝐂𝝃_0=0`$ and $`𝐂𝝃_1=0`$.

Finally, Ref.  shows that in situations where there are Jordan-chain modes, if one adds any small perturbation to the system then unstable modes are generically created. Hence, Jordan-chain modes should generically not occur in stable stars. However, Schutz notes that a necessary condition for this argument to apply is that a certain matrix element of the perturbation not vanish. Specifically, if $`\mathrm{\Delta }𝐓`$ is the perturbation to the operator (A10) and $`𝜻_{A,\sigma }`$ is a right Jordan chain with associated left Jordan chain $`𝝌_{A,\sigma }`$, then the argument requires

$$𝝌_{A,0},\mathrm{\Delta }𝐓𝜻_{A,0}0.$$
(D11)
In the current context, the argument is evaded because of the fact that the matrix element (D11) vanishes for arbitrary linear perturbations $`\mathrm{\Delta }𝐁`$ and $`\mathrm{\Delta }𝐂`$ to the operators $`𝐁`$ and $`𝐂`$, from Eqs. (A10) and (D4). Generic physical perturbations to the system do not correspond to generic mathematical perturbations to the operator $`𝐓`$, because of the structure of Eq. (A10).

## E The Zero Frequency Subspace for a Non-Rotating Star

This appendix characterizes the space $`_0`$ of zero frequency modes for non-rotating stars. The proofs in this appendix are Newtonian versions of the fully relativistic proofs given in Sec. III of Ref. .

The Lagrangian displacements $`𝝃(𝐱)`$ for zero frequency modes satisfy

$`𝐂𝝃=0,`$ (E1)
where the operator $`𝐂`$ is defined by Eq. (12) specialized to a non-rotating background star. From Eqs. (12)–(16), Eq. (E1) can be written as

$$\left[g_1(r)\delta \rho (𝐱)+\delta \varphi (𝐱)+g_2(r)\xi ^r(𝐱)\right]\beta (𝐱)g_2(r)𝐞_{\widehat{r}}=0,$$
(E2)
where

$$\beta (𝐱)𝝃,$$
(E3)
the radial component of $`𝝃`$ is $`\xi ^r𝐞_r`$, and the functions $`g_1(r)=\mathrm{\Gamma }_1p/\rho ^2`$ and $`g_2(r)=\mathrm{\Gamma }_1pA^r/\rho `$ are fixed functions determined by the background stellar model and by the adiabatic index $`\mathrm{\Gamma }_1`$ of the perturbations. The last two terms involving the function $`g_2`$ vanish for zero-buoyancy stars.

The argument of Lockitch and Friedman , simplified for Newtonian stars, is as follows. The perturbed pressure $`\overline{p}p+\delta p`$ and perturbed density $`\overline{\rho }=\rho +\delta \rho `$ obey the equation of hydrostatic balance and Poisson’s equation, to linear order in the perturbation. Hence, the pair of functions $`(\overline{\rho },\overline{p})`$ describe a static self-gravitating perfect fluid configuration, if one assumes that to each solution of the linearized fluid equations there corresponds a solution of the exact equations. However, it is known that all static, self-gravitating perfect fluid configurations are spherically symmetric†††††††††Note that the spherical-symmetry theorem does not require the barotropic assumption $`\overline{\rho }\times \overline{p}=0`$. It requires only that the perturbed pressure $`\overline{p}`$ and perturbed density $`\overline{\rho }`$ satisfy the equation of hydrostatic balance and Poisson’s equation.. Hence, all perturbations $`(\delta \rho ,\delta p)`$ must correspond to displacements to nearby static spherical configurations. Solutions of this type include (i) perturbations taking the star to a nearby static, spherically symmetric configuration with a different total mass $`M`$. If one assumes that the change $`\delta M`$ to the total mass of the star vanishes, such solutions are disallowed. (ii) Solutions with $`l=1`$ of the form $`𝝃=\mathrm{const}`$ corresponding to the displacement of the center of mass of the star. If one demands that the perturbation leaves fixed the center of mass of the star, these solutions are disallowed. (iii) Solutions with $`l=0`$ where the background star is marginally stable to radial collapse. Henceforth we will assume the background star is such that there are no $`l=0`$ zero-frequency solutions of this type.

If one makes sufficient assumptions to outlaw the cases (i), (ii) and (iii) above, it follows that all solutions to Eq. (E2) must satisfy

$$\delta \rho =\delta p=0.$$
(E4)
We now consider the cases of zero-buoyancy and nonzero-buoyancy stars.

### 1 Zero-buoyancy stars

We showed above that all solutions in $`_0`$ have $`\delta \rho =0`$. However, for $`𝐀=0`$, Eq. (E2) is always satisfied when $`\delta \rho =0`$, since the last two terms vanish. Hence, the space $`_0`$ consists precisely of the perturbations with

$`\delta \rho =(\rho 𝝃)=0.`$ (E5)
Note that this characterization does not hold for rotating stars, where there are zero frequency modes with $`\delta \rho 0`$ (see Appendix D).

We next derive a property of $`_0`$ that is used in Sec. III above. Let $`𝐏_0`$ to be the operator that projects orthogonally onto the subspace $`_0`$. We now show that

$`𝐏_0𝝃=0\text{if and only if}\times 𝝃=0.`$ (E6)
Let $`T^{ij}(𝐱)`$ be any antisymmetric tensor that vanishes on the boundary of the star. Then an integration by parts shows that

$$d^3x_{[i}\xi _{j]}^{}T^{ij}=𝝃,𝐯$$
(E7)
where the vector $`𝐯`$ is given by

$$v^i=\frac{1}{\rho }_jT^{ij}.$$
(E8)
Now if $`𝐏_0𝝃=0`$, then $`𝝃`$ is orthogonal to all vectors $`𝐯`$ in $`_0`$. But vectors of the form (E8) are automatically in $`_0`$, by Eq. (E5). Hence both sides of Eq. (E7) vanish. Since this is true for all $`T^{ij}(𝐱)`$ that vanish on the boundary, and since $`_i\xi _j`$ is continuous, it follows that $`\times 𝝃=0`$. It can be checked that the argument also carries through in the other direction.

### 2 Non-zero buoyancy stars

Turn now to stars with buoyancy forces. For these stars, the Lockitch-Friedman argument shows that the space of zero frequency modes is precisely the set of axial vectors (except possibly in the $`l=0`$ sector).

We use the expansion (143) of Lagrangian perturbation on a basis of vector spherical harmonics. Here the first two terms are the polar part of $`𝝃`$ and the last term is the axial part. It is easy to see that the axial part gives a vanishing contribution to both $`\delta \rho `$ and $`\beta =𝝃`$, and hence all axial vectors lie in $`_0`$, from Eq. (E2). Thus the $`C_{lm}`$ terms do not contribute at all to Eq. (E2) and so we can drop them and focus on the $`A_{lm}`$ and $`B_{lm}`$ terms. The vanishing of the density perturbation $`\delta \rho =(\rho 𝝃)`$ from Eq. (E4) can be written as

$$\frac{1}{\rho }(\rho A_{lm})^{}+2A_{lm}/rl(l+1)B_{lm}/r^2=0.$$
(E9)
Also, the relation between the Eulerian density and pressure perturbations for $`𝐀0`$ is, from Eq. (4),

$$\delta p=\frac{\mathrm{\Gamma }_1p}{\rho }\left(\delta \rho +\rho 𝝃𝐀\right).$$
(E10)
Substituting from Eq. (E4) now shows that the radial component of $`𝝃`$ vanishes. Hence $`A_{lm}=0`$ for all $`l,m`$, and it follows from Eq. (E9) that $`B_{lm}=0`$ for all $`l,m`$ with $`l>0`$ too. Since we are restricting attention to stars for which there are no $`l=0`$ zero-frequency modes (see above), it follows that the space $`_0`$ of zero frequency modes is therefore precisely the set of axial vectors.

Lastly, let $`𝐏_0`$ be the orthogonal projection operator that projects onto $`_0`$. In this case it follows from the definition of axial \[Eq. (143) with $`A_{lm}=B_{lm}=0`$\] that

$`𝐏_0𝝃=0\text{if and only if}𝐫\left(\times 𝝃\right)=0.`$ (E11)
## F Non-Existence of Modes Whose Frequencies Scale as the Square Root of the Rotational Frequency

In this appendix, we show that there are no modes whose frequencies and mode functions can be expressed as power series in the square root $`\sqrt{\mathrm{\Omega }}`$ of the star’s angular velocity and whose frequency in the non-rotating limit is zero. As discussed in Sec. III A above, one might suspect the existence of such modes from the general perturbation theory analysis of  and the fact that non-rotating stars have Jordan chains of length one. The result is valid for both zero-buoyancy and nonzero-buoyancy stars. The non-existence of such modes is related to the vanishing of the matrix element (D11) discussed in Appendix D above.

We assume a one parameter family of modes $`𝝃(\mathrm{\Omega }),\omega (\mathrm{\Omega })`$ with expansions of the form

$`\omega (\mathrm{\Omega })=\sqrt{\mathrm{\Omega }}\omega ^{(1/2)}+\mathrm{\Omega }\omega ^{(1)}+O(\mathrm{\Omega }^{3/2}),`$ (F1)
and

$`𝝃(\mathrm{\Omega })=𝝃^{(0)}+\sqrt{\mathrm{\Omega }}𝝃^{(1/2)}+\mathrm{\Omega }𝝃^{(1)}+O(\mathrm{\Omega }^{3/2}).`$ (F2)
If we substitute the expansions (F1)–(F2) and also the expansions (79)–(80) for the operators $`𝐁`$ and $`𝐂`$ into the quadratic eigenvalue equation (19) we get at orders $`O(\mathrm{\Omega }^0)`$, $`O(\sqrt{\mathrm{\Omega }})`$, and $`O(\mathrm{\Omega })`$ the equations $`𝐂^{(0)}𝝃^{(0)}=0`$, $`𝐂^{(0)}𝝃^{(1)}=0`$, and

$$𝐂^{(0)}𝝃^{(1)}\omega ^{(1/2)\mathrm{\hspace{0.17em}2}}𝝃^{(0)}=0.$$
(F3)
We now multiply Eq. (F3) by the projection operator $`𝐏_0`$ into the space $`_0`$ of zero-frequency modes, and use the fact that $`𝐏_0𝝃^{(0)}=𝝃^{(0)}`$ and $`𝐏_0𝐂^{(0)}=0`$. This gives

$$w^{(1/2)\mathrm{\hspace{0.17em}2}}𝝃^{(0)}=0,$$
(F4)
which shows that $`\omega ^{(1/2)}=0`$.

## G Derivation of Selection Rules

In this appendix we derive the two axial selection rules (157) and (158) and the selection rule (159) for an incompressible star. Throughout this appendix we will drop all terms that are not of zeroth order in $`\mathrm{\Omega }`$. Therefore we will use the pressure, density and potential profiles of a spherical background star.

### 1 Axial selection rules

We assume that the modes $`𝝃_B`$ and $`𝝃_C`$ are axial, and therefore

$$\mathbf{}𝝃_B=\mathbf{}𝝃_C=\delta \rho _B=\delta \rho _C=\xi _B^r=\xi _C^r=0.$$
(G1)
In addition we assume that the Eulerian density perturbation $`\delta \rho _A`$ of mode $`𝝃_A`$ vanishes, and so from Eq. (3) we have

$$𝝃_A=\frac{1}{\rho }\frac{d\rho }{dr}𝝃_A^r.$$
(G2)
Inserting these results into the expression (134) for the coupling coefficient we obtain

$`2\kappa _{ABC}`$ $`=`$ $`{\displaystyle \mathrm{d}^3xp\left[(\mathrm{\Gamma }_11)(\mathbf{}𝝃_A)\mathrm{\Xi }_{BC}+\chi _{ABC}+\chi _{ACB}\right]}{\displaystyle d^3x\rho \xi _A^i\xi _B^j\xi _C^k\varphi _{;ijk}}.`$ (G3)
Now the background potential $`\varphi `$ is a function of $`r`$ only, so we have

$$\varphi _{;jki}=\frac{1}{r}\frac{d}{dr}\left(\frac{1}{r}\frac{d\varphi }{dr}\right)(x_i\delta _{jk}+x_k\delta _{ij}+x_j\delta _{ik})+\frac{1}{r}\frac{d}{dr}\left[\frac{1}{r}\frac{d}{dr}\left(\frac{1}{r}\frac{d\varphi }{dr}\right)\right]x_ix_jx_k.$$
(G4)
Substituting Eq. (G4) into Eq. (G3) and using Eq. (G1) gives

$`2\kappa _{ABC}`$ $`=`$ $`{\displaystyle \mathrm{d}^3xp\left(\chi _{ABC}+\chi _{ACB}\right)}+{\displaystyle \mathrm{d}^3xp(\mathrm{\Gamma }_11)(\mathbf{}𝝃_A)\mathrm{\Xi }_{BC}}`$ (G6)
$`{\displaystyle \mathrm{d}^3x\rho \frac{d}{dr}\left(\frac{1}{r}\frac{d\varphi }{dr}\right)\xi _A^r\left(𝝃_B𝝃_C\right)}.`$
We now evaluate the first term in Eq. (G6) using the definition (139) of $`\chi _{ABC}`$. This gives

$$\mathrm{d}^3xp\chi _{ABC}=\mathrm{d}^3xp\xi _{A;j}^i\xi _{B;k}^j\xi _{C;i}^k=\mathrm{d}^3xp(\xi _A^i\xi _{B;k}^j)_{;j}\xi _{C;i}^k,$$
(G7)
where the second form of the integral was obtained using $`\xi _{B;kj}^j=\xi _{B;jk}^j=0`$. Upon integrating by parts†‡‡†‡‡†‡‡ Throughout this appendix we discard surface terms that arise from integrations by parts which involve the radial component $`\xi _A^r`$ of the mode function $`𝝃_A`$ evaluated on the star’s surface, because the quantity
$$\frac{\mathrm{\Delta }p}{\rho }=\frac{\mathrm{\Gamma }_1p}{\rho ^2}\frac{d\rho }{dr}\xi _A^r=\frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}\varphi ^{}\xi _A^r$$
must vanish on the surface . we find

$$\mathrm{d}^3xp\chi _{ABC}=\mathrm{d}^3x(p_{,j}\xi _A^i\xi _{B;k}^j\xi _{C;i}^k+p\xi _A^i\xi _{B;k}^j\xi _{C;ij}^k).$$
(G8)
We can rewrite this using the relationship that follows from $`p=p(r)`$

$$(p_{,j}\xi _B^j)_{;k}=0=p_{,j}\xi _{B;k}^j+p_{;jk}\xi _B^j;$$
(G9)
the result is

$`{\displaystyle \mathrm{d}^3xp\chi _{ABC}}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x\left[p_{;jk}\xi _A^i\xi _B^j\xi _{C;i}^kp\xi _A^i\xi _{B;k}^j\xi _{C;ij}^k\right]}`$ (G10)
$`=`$ $`{\displaystyle \mathrm{d}^3x\left[p_{;jki}\xi _A^i\xi _B^j\xi _C^k+p_{;jk}\xi _{A;i}^i\xi _B^j\xi _C^k+p_{;jk}\xi _A^i\xi _{B;i}^j\xi _C^k+p\xi _A^i\xi _{B;k}^j\xi _{C;ij}^k\right]},`$ (G11)
where we have again integrated by parts. We can evaluate the term involving $`p_{;jki}`$ using an equation analogous to Eq. (G4), which gives

$`{\displaystyle \mathrm{d}^3xp\chi _{ABC}}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x\left[\frac{d}{dr}\left(\frac{1}{r}\frac{dp}{dr}\right)\xi _A^r\left(𝝃_A𝝃_B\right)+p_{;jk}\xi _{A;i}^i\xi _B^j\xi _C^k+p_{;jk}\xi _A^i\xi _{B;i}^j\xi _C^k+p\xi _A^i\xi _{B;k}^j\xi _{C;ij}^k\right]},`$ (G12)
Next, we symmetrize over the indices $`B`$ and $`C`$. This gives

$`{\displaystyle \mathrm{d}^3xp\left(\chi _{ABC}+\chi _{ACB}\right)}`$ $`=`$ $`2{\displaystyle \mathrm{d}^3x\left[\frac{d}{dr}\left(\frac{1}{r}\frac{dp}{dr}\right)\xi _A^r\left(𝝃_A𝝃_B\right)+p_{;jk}\xi _{A;i}^i\xi _B^j\xi _C^k\right]}`$ (G14)
$`{\displaystyle \mathrm{d}^3x\left[p_{;jk}\xi _A^i\xi _{B;i}^j\xi _C^k+p\xi _A^i\xi _{B;k}^j\xi _{C;ij}^k+p_{;jk}\xi _A^i\xi _{C;i}^k\xi _B^j+p\xi _A^i\xi _{C;j}^k\xi _{B;ik}^j\right]},`$
where on the right hand side we have relabeled some indices. We can rewrite the second line of Eq. (G14) as

$`{\displaystyle \mathrm{d}^3x\left[p_{;jk}\xi _A^i\left(\xi _B^j\xi _C^k\right)_{;i}+p\xi _A^i\left(\xi _{B;k}^j\xi _{C;j}^k\right)_{;i}\right]},`$ (G15)
$`=`$ $`{\displaystyle \mathrm{d}^3x[p_{;jki}\xi _A^i\xi _B^j\xi _C^k+p_{;jk}\xi _{A;i}^i\xi _B^j\xi _C^k+p,i\xi _A^i\xi _{B;k}^j\xi _{C;j}^k+p\xi _{A;}^i\xi _{B;k}^j\xi _{C;j}^k]},`$ (G16)
where we have again integrated by parts. Rewriting the terms involving $`p_{;jki}`$ and $`p_{;jk}`$ the same way as before, and using the definition (138) of $`\mathrm{\Xi }_{BC}`$, we can write the expression (G14) as

$`{\displaystyle \mathrm{d}^3xp\left(\chi _{ABC}+\chi _{ACB}\right)}`$ $`=`$ $`{\displaystyle \mathrm{d}^3x\left[\frac{d}{dr}\left(\frac{1}{r}\frac{dp}{dr}\right)\xi _A^r\left(𝝃_B𝝃_C\right)+\frac{1}{r}\frac{dp}{dr}\left(𝝃_A\right)\left(𝝃_B𝝃_C\right)\right]}`$ (G18)
$`+{\displaystyle \mathrm{d}^3x\left[p_{,i}\xi _A^i\mathrm{\Xi }_{BC}+p\xi _{A;i}^i\mathrm{\Xi }_{BC}\right]}.`$
Now since the Eulerian density perturbation $`(\rho \xi _A^i)_{;i}`$ vanishes, we have

$$(p\xi _A^i)_{;i}=\frac{dp}{d\rho }\rho _{;i}\xi _A^i+p\xi _{A;i}^i=p\xi _{A;i}^i(\frac{d\mathrm{ln}p}{d\mathrm{ln}\rho }+1)=p(\mathrm{\Gamma }1)\xi _{A;i}^i,$$
(G19)
where $`\mathrm{\Gamma }`$ is the adiabatic index of the background spherical star. Combining Eqs. (G6), (G16), and (G18) now yields

$`2\kappa _{ABC}`$ $`=`$ $`{\displaystyle d^3x\left(𝝃_B𝝃_C\right)\left[\rho \xi _A^r\frac{d}{dr}\left(\frac{1}{r}\frac{d\varphi }{dr}\right)\xi _A^r\frac{d}{dr}\left(\frac{1}{r}\frac{dp}{dr}\right)\frac{1}{r}\frac{dp}{dr}\left(𝝃_A\right)\right]}`$ (G21)
$`+{\displaystyle d^3xp(\mathrm{\Gamma }_1\mathrm{\Gamma })\left(𝝃_A\right)\mathrm{\Xi }_{BC}}.`$
Using Eq. (G2) and the equation of hydrostatic equilibrium in the background star one can show that the first line of the right hand side of Eq. (G21) vanishes. Therefore we find

$$\kappa _{ABC}=\frac{1}{2}d^3xp(\mathrm{\Gamma }_1\mathrm{\Gamma })\left(𝝃_A\right)\mathrm{\Xi }_{BC}.$$
(G22)
We can derive two different selection rules from Eq. (G22). First, if $`𝝃_A`$ is axial, then $`𝝃_A=0`$ and therefore $`\kappa _{ABC}=0`$; this proves the selection rule (157) for three axial modes. Second, if the perturbations obey the same equation of state as the background star, then the adiabatic index $`\mathrm{\Gamma }_1`$ of the perturbations coincides with the adiabatic index $`\mathrm{\Gamma }`$ of the background star, and $`\kappa _{ABC}`$ again vanishes. This proves the second axial selection rule (158).

### 2 Incompressible selection rule

We can consider an incompressible star to be the limit $`\mathrm{\Gamma }_1\mathrm{}`$ of a zero-buoyancy polytrope with $`\mathrm{\Gamma }=\mathrm{\Gamma }_1`$. The expressions (134) and (142) for the coupling coefficient $`\kappa _{ABC}`$ are valid for incompressible stars, since they are valid for finite $`\mathrm{\Gamma }_1`$ and they have a finite limit as $`\mathrm{\Gamma }_1\mathrm{}`$. Note however that one must be careful to include all the distributional contributions in the integrands at the stellar surface that arise in the $`\mathrm{\Gamma }_1\mathrm{}`$ limit.

Consider now a linear mode $`𝝃(𝐱)`$ of the large-$`\mathrm{\Gamma }_1`$ star, whose $`\mathrm{\Gamma }_1\mathrm{}`$ limit is a mode of the incompressible star. Clearly the Lagrangian perturbation of the density must scale as

$$\frac{\mathrm{\Delta }\rho }{\rho }\frac{1}{\mathrm{\Gamma }_1},$$
(G23)
which implies from Eq. (I18) that

$$\mathbf{}𝝃\frac{1}{\mathrm{\Gamma }_1}.$$
(G24)
It follows that to zeroth order in $`1/\mathrm{\Gamma }_1`$, the Eulerian density perturbation is, from Eq. (3)

$`\delta \rho `$ $`=`$ $`(\rho 𝝃)=𝝃\rho \rho \left(\mathbf{}𝝃\right)`$ (G25)
$`=`$ $`\xi ^r\stackrel{~}{\rho }\delta (rR)`$ (G26)
where $`\stackrel{~}{\rho }=3M/4\pi R^3`$ and $`M`$ and $`R`$ are the stellar mass and radius.

Consider now three modes $`𝝃_A`$, $`𝝃_B`$, $`𝝃_C`$ that satisfy the hypothesis of the selection rule (159). That is, they all have vanishing Eulerian density and pressure perturbations to zeroth order in $`1/\mathrm{\Gamma }_1`$. It follows from Eq. (G26) that they all have vanishing radial components at the stellar surface,

$$\xi ^r(r=R)=0.$$
(G27)
Consider now the expression (142) for the coupling coefficient. In this expression, the second, third and fifth terms vanish since $`\mathbf{}𝝃=0`$ for each mode, and the fourth and sixth terms vanish since each mode has vanishing $`\delta \rho `$, $`\delta \varphi `$ and $`\delta p`$. The remaining terms yield the following expression for the coupling coefficient:

$$\kappa _{ABC}=\frac{1}{2}d^3x\xi _A^i\xi _B^j\xi _C^k\left[p_{;ijk}+\rho \varphi _{;ijk}\right].$$
(G28)
The integrand in Eq. (G28) vanishes in the interior of the star, since $`p`$ and $`\varphi `$ are quadratic functions of position for a constant density star. However, there are distributional contributions to the integral at the stellar surface, since $`p_{;i}`$ and $`\varphi _{;ij}`$ are discontinuous across the surface. To evaluate these contributions, substitute in the equation of hydrostatic equilibrium in the background star, $`p_{;i}+\rho \varphi _{;i}=0`$, which yields

$$\kappa _{ABC}=\frac{1}{2}d^3x\xi _{(A}^i\xi _B^j\xi _{C)}^k\left[2\rho _{;i}\varphi _{;jk}+\varphi _{;i}\rho _{;jk}\right].$$
(G29)
The first term vanishes since $`\varphi _{;ji}`$ is finite on the stellar surface, and for any of the modes we have

$$\xi ^k\rho _{;k}=\stackrel{~}{\rho }\xi ^r(r=R)=0$$
(G30)
from Eq. (G27). Integrating the second term by parts gives

$$\kappa _{ABC}=\frac{1}{2}d^3x\left[\xi _{(A}^i\xi _B^j\xi _{C)}^k\varphi _{;i}\right]_{;j}\rho _{;k}.$$
(G31)
When the derivative of the square bracket is expanded out, each term will contain either a factor $`\xi ^i\rho _{;i}`$ or a factor $`\xi ^i\varphi _{;i}`$, both of which vanish on the stellar surface. Therefore we obtain

$$\kappa _{ABC}=0$$
(G32)
to zeroth order in $`1/\mathrm{\Gamma }_1`$ and in $`\mathrm{\Omega }`$. Note that the cancellation of the distributional components in the expression (G28) would not be obtained if one uses the Cowling approximation of neglecting the gravitational terms in the coupling coefficient, as in Ref. .

## H Calculating Covariant Derivatives Using the Spin Weighted Spherical Harmonics and the $`ð`$ Operator

We define the “spin-weight” $`s`$ of any tensor as in  (also see  for a very readable reference), in terms of how it transforms under rotations around the radial direction. For example, the vector $`𝐦`$ of Sec. VI B transforms under these types of rotations as $`𝐦𝐦^{}`$, where

$`𝐦^{}=e^{i\psi }𝐦`$ (H1)
where $`\psi `$ is the angle of rotation. We say that this vector has spin-weight $`s=1`$. In general, a quantity $`\vartheta `$ is said to have spin-weight $`s`$ if it transforms as $`\vartheta \vartheta ^{}`$ where

$`\vartheta ^{}=e^{is\psi }\vartheta .`$ (H2)
In particular, the basis vectors $`𝐥`$,$`𝐦`$ and $`\overline{𝐦}`$ defined in Sec. VI B have spin-weights $`0,1`$ and $`1`$ respectively. Thus we can decompose any given vector function, $`𝐁(r,\theta ,\phi )`$, into three components with definite spin-weight by dotting it with one of the three basis vectors $`(𝐥,𝐦,\overline{𝐦})`$. Each scalar component can then be expanded in the complete set of basis functions $`{}_{s}{}^{}Y_{lm}^{}`$, e.g.

$`𝐁𝐥`$ $``$ $`\zeta _0(r,\theta ,\phi )={\displaystyle \underset{\mathrm{\Lambda }}{}}b_0^\mathrm{\Lambda }(r){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (H3)
$`𝐁𝐦`$ $``$ $`\zeta _1(r,\theta ,\phi )={\displaystyle \underset{\mathrm{\Lambda }}{}}b_+^\mathrm{\Lambda }(r){}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (H4)
$`𝐁\overline{𝐦}`$ $``$ $`\zeta _1(r,\theta ,\phi )={\displaystyle \underset{\mathrm{\Lambda }}{}}b_{}^\mathrm{\Lambda }(r){}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}.`$ (H5)
For notational brevity we note that $`𝐦`$ and $`\overline{𝐦}`$, and consequently $`\zeta _1`$ and $`\zeta _1`$, are complex conjugates of each other. Thus, in subsequent calculations we will write $`+c.c.`$ when convenient.‡\*‡\*‡\*Note that in Section VI B we choose to implement this by simply adding terms with $`s`$ replaced by $`s`$. The connection coefficients on this basis are

$`\mathrm{\Gamma }_{m\overline{m}}^l`$ $`=\mathrm{\Gamma }_{lm}^m={\displaystyle \frac{1}{r}}`$ (H6)
$`\mathrm{\Gamma }_{m\overline{m}}^m`$ $`=\mathrm{\Gamma }_{mm}^m={\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \frac{\mathrm{cot}\theta }{r}},`$ (H7)
where all the other components are zero except for complex conjugates of the above, e.g.

$`(\mathrm{\Gamma }_{mm}^m)^{}=\mathrm{\Gamma }_{\overline{m}\overline{m}}^{\overline{m}}={\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \frac{\mathrm{cot}\theta }{r}}.`$ (H8)
It follows that the components of the covariant derivative of $`𝐁`$ are

$`B_{;l}^l`$ $`=`$ $`B_{,l}^l`$ (H9)
$`B_{;m}^m`$ $`=`$ $`B_{,m}^m+{\displaystyle \frac{1}{r}}B^l+{\displaystyle \frac{\mathrm{cot}\theta }{\sqrt{2}r}}B^m`$ (H10)
$`B_{;\overline{m}}^m`$ $`=`$ $`B_{,\overline{m}}^m{\displaystyle \frac{\mathrm{cot}\theta }{\sqrt{2}r}}B^m`$ (H11)
$`B_{;l}^m`$ $`=`$ $`B_{,l}^m`$ (H12)
$`B_{;m}^l`$ $`=`$ $`B_{,m}^l{\displaystyle \frac{1}{r}}B^{\overline{m}},`$ (H13)
where all others are zero except for complex conjugates of the above. Now consider the $`ð`$ (“edth”) and $`\overline{ð}`$ operators of Ref. , defined by their action on a function, $`\zeta _s`$ of spin-weight $`s`$,

$`ð\zeta _s`$ $`=`$ $`(\mathrm{sin}\theta )^s\sqrt{2}(𝐦𝐋)(\mathrm{sin}\theta )^s\zeta _s=\left(_\theta +i\mathrm{csc}\theta _\phi s\mathrm{cot}\theta \right)\zeta _s`$ (H14)
$`\overline{ð}\zeta _s`$ $`=`$ $`(\mathrm{sin}\theta )^s\sqrt{2}(\overline{𝐦}𝐋)(\mathrm{sin}\theta )^s\zeta _s=\left(_\theta +i\mathrm{csc}\theta _\phi s\mathrm{cot}\theta \right)\zeta _s.`$ (H15)
Using these definitions in Eq. (H13) we get

$`B_{;m}^m`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}r}}ð\zeta _1+{\displaystyle \frac{1}{r}}\zeta _0`$ (H16)
$`B_{;m}^l`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}r}}ð\zeta _0{\displaystyle \frac{1}{r}}\zeta _1`$ (H17)
$`B_{;\overline{m}}^m`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}r}}\overline{ð}\zeta _1`$ (H18)
$`B_{;l}^m`$ $`=`$ $`_r\zeta _1`$ (H19)
$`B_{;l}^l`$ $`=`$ $`_r\zeta _0`$ (H20)
and complex conjugates of the above. The $`ð`$ ($`\overline{ð}`$) operator, when acting on the spin-weighted spherical harmonics, raises (lowers) the spin-weight by one,

$`ð{}_{s}{}^{}Y_{lm}^{}`$ $`=`$ $`\sqrt{(ls)(l+s+1)}{}_{s+1}{}^{}Y_{lm}^{}`$ (H21)
$`\overline{ð}{}_{s}{}^{}Y_{lm}^{}`$ $`=`$ $`\sqrt{(l+s)(ls+1)}{}_{s1}{}^{}Y_{lm}^{}.`$ (H22)
Thus when calculating a covariant derivative such as $`B_{;m}^m`$ we can write

$`B_{;m}^m={\displaystyle \frac{1}{r}}{\displaystyle \underset{\mathrm{\Lambda }}{}}\left(b_0^\mathrm{\Lambda }(r){\displaystyle \frac{\mathrm{\Lambda }}{\sqrt{2}}}b_{}^\mathrm{\Lambda }(r)\right){}_{0}{}^{}Y_{\mathrm{\Lambda }}^{},`$ (H23)
and similarly for the other covariant derivatives. This allows us to define the functionals as in Eq. (183). For decomposing vectors of the form (161), note that

$`𝐦{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{\sqrt{2}r}}{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (H24)
$`\overline{𝐦}{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{\sqrt{2}r}}{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (H25)
$`𝐦𝐋{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{\sqrt{2}}}{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}`$ (H26)
$`\overline{𝐦}𝐋{}_{0}{}^{}Y_{\mathrm{\Lambda }}^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Lambda }}{\sqrt{2}}}{}_{1}{}^{}Y_{\mathrm{\Lambda }}^{}.`$ (H27)
Using these relations together with Eqs. (149) and (150) we find the functions $`f_0^\mathrm{\Lambda }`$, $`f_+^\mathrm{\Lambda }`$ and $`f_{}^\mathrm{\Lambda }`$ for $`z`$parity odd modes

$`f_0^{lm}(r)`$ $`=`$ $`{\displaystyle \frac{W_{j+m+1}}{r}}\delta _{j+m+1,l}`$ (H28)
$`f_+^{lm}(r)`$ $`=`$ $`{\displaystyle \frac{V_{j+m+1}}{r}}{\displaystyle \frac{\sqrt{(j+m+1)(j+m+2)}}{\sqrt{2}}}\delta _{j+m+1,l}`$ (H30)
$`+{\displaystyle \frac{U_{j+m}}{r}}{\displaystyle \frac{\sqrt{(j+m)(j+m+1)}}{\sqrt{2}}}\delta _{j+m,l}`$
$`f_{}^{lm}(r)`$ $`=`$ $`{\displaystyle \frac{V_{j+m+1}}{r}}{\displaystyle \frac{\sqrt{(j+m+1)(j+m+2)}}{\sqrt{2}}}\delta _{j+m+1,l}`$ (H32)
$`+{\displaystyle \frac{U_{j+m}}{r}}{\displaystyle \frac{\sqrt{(j+m)(j+m+1)}}{\sqrt{2}}}\delta _{j+m,l},`$
and for $`z`$parity even modes

$`f_0^{lm}(r)`$ $`=`$ $`{\displaystyle \frac{W_{j+m}}{r}}\delta _{j+m,l}`$ (H33)
$`f_+^{lm}(r)`$ $`=`$ $`{\displaystyle \frac{V_{j+m}}{r}}{\displaystyle \frac{\sqrt{(j+m)(j+m+1)}}{\sqrt{2}}}\delta _{j+m,l}`$ (H35)
$`+{\displaystyle \frac{U_{j+m+1}}{r}}{\displaystyle \frac{\sqrt{(j+m+1)(j+m+2)}}{\sqrt{2}}}\delta _{j+m+1,l}`$
$`f_{}^{lm}(r)`$ $`=`$ $`{\displaystyle \frac{V_{j+m}}{r}}{\displaystyle \frac{\sqrt{(j+m)(j+m+1)}}{\sqrt{2}}}\delta _{j+m,l}`$ (H37)
$`+{\displaystyle \frac{U_{j+m+1}}{r}}{\displaystyle \frac{\sqrt{(j+m+1)(j+m+2)}}{\sqrt{2}}}\delta _{j+m+1,l}.`$
## I Direct Derivation of Second Order Equations of Motion

There are at least three ways to derive the equation of motion for the second order Lagrangian perturbations: (1) Expand the variational principle to third order in the Lagrangian displacement $`𝝃`$, and then vary with respect to $`𝝃`$; (2) Vary with respect to $`𝝃`$ to get the equation of motion, then expand to second order; (3) Directly perturb the equation of motion to second order in $`𝝃`$:

$`\mathrm{\Delta }\left({\displaystyle \frac{du_i}{dt}}+{\displaystyle \frac{1}{\rho }}_ip+_i\varphi \right)=0.`$ (I1)
The second method is carried out in Appendix J. Here we carry out the third method.

Recall that the Eulerian and Lagrangian perturbations of a quantity $`Q`$ are defined as

$`\delta Q`$ $`=`$ $`Q(𝐱,t)Q_0(𝐱,t),`$ (I2)
$`\mathrm{\Delta }Q`$ $`=`$ $`Q(𝐱+𝝃(𝐱,t),t)Q_0(𝐱,t),`$ (I3)
where $`Q_0`$ is the value of $`Q`$ in the unperturbed solution. \[In this appendix we shall sometimes use $`\rho _0`$, $`p_0`$ and $`\varphi _0`$ to denote the background density, pressure and potential; these quantities were denoted simply $`\rho `$, $`p`$ and $`\varphi `$ in the body of the paper.\] Suppressing the time dependence of the quantities, we can write to second order

$`\mathrm{\Delta }Q`$ $`=`$ $`Q(𝐱+𝝃)Q_0(𝐱)`$ (I4)
$`=`$ $`Q(𝐱+𝝃)Q_0(𝐱+𝝃)+Q_0(𝐱+𝝃)Q_0(𝐱)`$ (I5)
$`=`$ $`\delta Q(𝐱+𝝃)+Q_0(𝐱+𝝃)Q_0(𝐱)`$ (I6)
$`=`$ $`\delta Q+\xi ^i_i\delta Q+\xi ^i_iQ+{\displaystyle \frac{1}{2}}\xi ^i\xi ^j_i_jQ.`$ (I7)
In the last line we have dropped the subscript 0 on the last two terms. Note that $`\delta Q`$ in the second term on the right-hand side need only be evaluated to first order. Substituting the first order relation

$`\delta Q=\mathrm{\Delta }Q\xi ^i_iQ`$ (I8)
in this term gives the alternative relation

$`\mathrm{\Delta }Q=\delta Q+\xi ^i_i\mathrm{\Delta }Q+\xi ^i_iQ\xi ^i_i\xi ^j_jQ{\displaystyle \frac{1}{2}}\xi ^i\xi ^j_i_jQ.`$ (I9)
By definition of the Lagrangian displacement vector the unperturbed and perturbed fluid trajectories are related by

$`x_0^i=c^i(t),x^i=c^i(t)+\xi ^i(c^j(t),t).`$ (I10)
Thus

$`u_0^i={\displaystyle \frac{dc^i}{dt}},u^i={\displaystyle \frac{dc^i}{dt}}+{\displaystyle \frac{dc^j}{dt}}{\displaystyle \frac{\xi ^i}{x^j}}+{\displaystyle \frac{\xi ^i}{t}}.`$ (I11)
Hence,

$`\mathrm{\Delta }u^i={\displaystyle \frac{d\xi ^i}{dt}},`$ (I12)
an equation that is exact to all orders in $`𝝃`$ ‡†‡†‡† Note that we have defined the Lagrangian perturbation of the vector $`v^i`$ by subtracting the components at two locations separated by $`𝝃`$. This procedure is valid for Cartesian vector components in flat spacetime. We can derive formulas valid in arbitrary coordinate systems by working first in Cartesian coordinates, and then rewriting all partial derivatives as covariant derivatives at the end. Friedman and Schutz  have suggested defining the Lagrangian perturbation by mapping the perturbed vector back to the original location so that vectors are subtracted only at the same point (Lie derivative). This gives different intermediate values for the Lagrangian perturbations of various quantities than the procedure followed here. However, when the final equations of motion are written in terms of $`𝝃`$, the results are the same since one is using the same definition of $`𝝃`$. While the Lie derivative is essential in relativistic applications, the present approach is simpler here. (see also Ref. ).

An expression for $`\mathrm{\Delta }\rho `$ follows from conservation of mass. The variation in the total mass due to the perturbation is

$`\delta M`$ $`=`$ $`{\displaystyle _{V+\mathrm{\Delta }V}}\rho (𝐱)d^3x{\displaystyle _V}\rho _0(𝐱)d^3x`$ (I13)
$`=`$ $`{\displaystyle _V}[\rho (𝐱+𝝃)J\rho _0(𝐱)]d^3x.`$ (I14)
Here we have made the transformation $`𝐱^{}=𝐱𝝃(𝐱)`$ in the first integral to make the volumes the same, and then dropped the primes. The Jacobian of the transformation is

$`J={\displaystyle \frac{(𝐱+𝝃)}{(𝐱)}}=\mathrm{det}(\delta _j^i+L^i{}_{j}{}^{}),`$ (I15)
where

$`L^i{}_{j}{}^{}={\displaystyle \frac{\xi ^i}{x^j}}.`$ (I16)
Setting $`\delta M=0`$ and substituting $`\rho (𝐱+𝝃)=\rho _0(𝐱)+\mathrm{\Delta }\rho `$ in Eq. (I14), we get

$`{\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}={\displaystyle \frac{1J}{J}},`$ (I17)
exact to all orders in $`𝝃`$. Keeping terms through second order from Eq. (J20), we find

$`{\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}=_i\xi ^i+{\displaystyle \frac{1}{2}}[(_i\xi ^i)^2+_i\xi ^k_k\xi ^i].`$ (I18)
Consider an equation of state $`p=p(\rho ,\mu )`$ where $`\mu `$ denotes quantities such as entropy or composition. For adiabatic perturbations in which the composition does not have time to change $`\mathrm{\Delta }\mu =0`$, so to second order

$`\mathrm{\Delta }p={\displaystyle \frac{p}{\rho }}|_\mu \mathrm{\Delta }\rho +{\displaystyle \frac{1}{2}}{\displaystyle \frac{^2p}{\rho ^2}}|_\mu (\mathrm{\Delta }\rho )^2.`$ (I19)
In terms of the adiabatic index governing the perturbations,

$`\mathrm{\Gamma }_1={\displaystyle \frac{\mathrm{ln}p}{\mathrm{ln}\rho }}|_\mu ,`$ (I20)
we get

$`{\displaystyle \frac{\mathrm{\Delta }p}{p}}=\mathrm{\Gamma }_1{\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}+{\displaystyle \frac{1}{2}}\left[\mathrm{\Gamma }_1(\mathrm{\Gamma }_11)+\rho {\displaystyle \frac{\mathrm{\Gamma }_1}{\rho }}\right]\left({\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}\right)^2.`$ (I21)
Here $`\mathrm{\Delta }\rho `$ need be evaluated only to first order in the last term.

To derive the equation of motion for the perturbations, we need to perturb various derivatives of quantities. Consider for example $`Q/t`$. By Eq. (I7),

$`\mathrm{\Delta }{\displaystyle \frac{Q}{t}}=\delta {\displaystyle \frac{Q}{t}}+\xi ^i_i\delta {\displaystyle \frac{Q}{t}}+\xi ^i_i{\displaystyle \frac{Q}{t}}+{\displaystyle \frac{1}{2}}\xi ^i\xi ^j_i_j{\displaystyle \frac{Q}{t}}.`$ (I22)
The Eulerian perturbation operator commutes with partial differentiation:

$`\delta {\displaystyle \frac{}{t}}={\displaystyle \frac{}{t}}\delta ,`$ (I23)
so the first term on the right-hand side of Eq. (I22) is

$`{\displaystyle \frac{}{t}}\delta Q={\displaystyle \frac{}{t}}\left(\mathrm{\Delta }Q\xi ^i_i\delta Q\xi ^i_iQ{\displaystyle \frac{1}{2}}\xi ^i\xi ^j_i_jQ\right),`$ (I24)
where we have used Eq. (I7) again. Thus

$`\mathrm{\Delta }{\displaystyle \frac{Q}{t}}={\displaystyle \frac{}{t}}\mathrm{\Delta }Q{\displaystyle \frac{\xi ^i}{t}}\left(_i\delta Q+_iQ+\xi ^j_i_jQ\right).`$ (I25)
Here we have used the symmetry of $`_i_j`$ to rewrite the last term without the factor of $`1/2`$. An alternative form of Eq. (I25) follows from substituting Eq. (I8):

$`\mathrm{\Delta }{\displaystyle \frac{Q}{t}}={\displaystyle \frac{}{t}}\mathrm{\Delta }Q{\displaystyle \frac{\xi ^i}{t}}\left(_i\mathrm{\Delta }Q+_iQ_i\xi ^j_jQ\right).`$ (I26)
A similar equation holds if $`/t`$ is replaced by $`/x^i`$, and then we can let $`/x^i_i`$ to get a tensor equation valid in any coordinate system.

For the perturbation of the acceleration, we have

$`\mathrm{\Delta }{\displaystyle \frac{du_i}{dt}}`$ $`=`$ $`\mathrm{\Delta }\left({\displaystyle \frac{}{t}}+u^k_k\right)u_i`$ (I27)
$`=`$ $`\mathrm{\Delta }{\displaystyle \frac{}{t}}u_i+u^k\mathrm{\Delta }_ku_i+\mathrm{\Delta }u^k_ku_i+\mathrm{\Delta }u^k\mathrm{\Delta }_ku_i.`$ (I28)
Note that our formulas hold in non-Cartesian coordinates whether we perturb the covariant components $`u_i`$ or the contravariant components $`u^i`$, for the reason explained in footnote ‡†I. Commute the operators $`\mathrm{\Delta }/t`$ and $`\mathrm{\Delta }_k`$ in Eq. (I28) using Eq. (I25) and the corresponding equation for $`_k`$. Many terms cancel, and we are left with

$`\mathrm{\Delta }{\displaystyle \frac{du_i}{dt}}={\displaystyle \frac{d^2\xi _i}{dt^2}}.`$ (I29)
For a time-independent unperturbed star $`u^i/t=0`$, and so

$`\mathrm{\Delta }{\displaystyle \frac{du_i}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{}{t}}+u^j_j\right)\left({\displaystyle \frac{}{t}}+u^k_k\right)\xi _i`$ (I30)
$`=`$ $`{\displaystyle \frac{^2\xi _i}{t^2}}+2u^j_j{\displaystyle \frac{\xi _i}{t}}+u^j_j(u^k_k\xi _i).`$ (I31)
Using Eq. (I7) and commuting $`\delta `$ and $`_i`$, we get

$`\mathrm{\Delta }_i\varphi =_i\delta \varphi +\xi ^j_j_i\delta \varphi +\xi ^j_j_i\varphi +{\displaystyle \frac{1}{2}}\xi ^j\xi ^k_j_k_i\varphi .`$ (I32)
The Eulerian perturbation is found from Poisson’s equation:

$`^2\varphi =4\pi G\rho ^2\delta \varphi =4\pi G\delta \rho ,`$ (I33)
which implies

$`\delta \varphi =G{\displaystyle \frac{\delta \rho ^{}}{|𝐱𝐱^{}|}d^3x^{}}.`$ (I34)
We get an expression for $`\delta \rho `$ by substituting Eq. (I18) in Eq. (I9):

$`\delta \rho `$ $`=`$ $`\rho _i\xi ^i+{\displaystyle \frac{1}{2}}\rho [(_i\xi ^i)^2+_i\xi ^k_k\xi ^i]`$ (I36)
$`+\xi ^i_i(\rho _j\xi ^j)\xi ^i_i\rho +\xi ^i_i\xi ^j_j\rho +{\displaystyle \frac{1}{2}}\xi ^i\xi ^j_i_j\rho .`$
After some algebra, we can simplify this to ‡‡‡‡‡‡Note that the second term in Eq. (I37) typically includes a $`\delta `$-function at the stellar surface coming from the second derivative of the density. This $`\delta `$-function is physical and should be included; it corrects for the fact that the expression (I37) vanishes outside the surface of the unperturbed star whereas in reality $`\delta \rho `$ can be nonvanishing there.

$`\delta \rho =_i(\rho \xi ^i)+{\displaystyle \frac{1}{2}}_i_j(\rho \xi ^i\xi ^j).`$ (I37)
We have

$`\mathrm{\Delta }\left({\displaystyle \frac{1}{\rho }}_ip\right)`$ $`=`$ $`\mathrm{\Delta }\left({\displaystyle \frac{1}{\rho }}\right)_ip+{\displaystyle \frac{1}{\rho }}\mathrm{\Delta }_ip+\mathrm{\Delta }\left({\displaystyle \frac{1}{\rho }}\right)\mathrm{\Delta }_ip`$ (I38)
$`=`$ $`\left[{\displaystyle \frac{1}{\rho ^2}}\mathrm{\Delta }\rho +{\displaystyle \frac{1}{\rho ^3}}(\mathrm{\Delta }\rho )^2\right]_ip+{\displaystyle \frac{1}{\rho }}\mathrm{\Delta }_ip{\displaystyle \frac{1}{\rho ^2}}\mathrm{\Delta }\rho \mathrm{\Delta }_ip.`$ (I39)
Substitute Eq. (I26) with $`_ip`$ replacing $`Q/t`$:

$`\rho \mathrm{\Delta }\left({\displaystyle \frac{1}{\rho }}_ip\right)`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}+\left({\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}\right)^2\right]_ip+_i\mathrm{\Delta }p`$ (I40)
$``$ $`_i\xi ^j(_j\mathrm{\Delta }p_j\xi ^k_kp+_jp)`$ (I41)
$`{\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}(_i\mathrm{\Delta }p_i\xi ^j_jp)`$ (I42)
Recall the definitions of the tensors $`\mathrm{\Theta }^j_i`$ and $`\mathrm{\Xi }^j_i`$ and their traces from Eq. (117). Their divergences satisfy the identity

$`_j\mathrm{\Xi }^j{}_{i}{}^{}_j\mathrm{\Theta }^j{}_{i}{}^{}={\displaystyle \frac{1}{2}}_i(\mathrm{\Xi }\mathrm{\Theta }).`$ (I43)
Substitute Eqs. (I18) and (I21) in Eq. (I40) to get

$`\rho \mathrm{\Delta }\left({\displaystyle \frac{1}{\rho }}_ip\right)`$ $`=`$ $`_j\xi ^j_ip_i\xi ^j_jp_i(p\mathrm{\Gamma }_1_j\xi ^j)`$ (I44)
$`+{\displaystyle \frac{1}{2}}(\mathrm{\Theta }\mathrm{\Xi })_ip`$ (I45)
$`+`$ $`_i\left[p\mathrm{\Gamma }_1{\displaystyle \frac{1}{2}}(\mathrm{\Theta }+\mathrm{\Xi })+{\displaystyle \frac{1}{2}}p\mathrm{\Theta }\left(\mathrm{\Gamma }_1^2\mathrm{\Gamma }_1+\rho {\displaystyle \frac{\mathrm{\Gamma }_1}{\rho }}\right)\right]`$ (I47)
$`+(\mathrm{\Xi }^k{}_{i}{}^{}\mathrm{\Theta }^k{}_{i}{}^{})_kp+_i\xi ^j_j(p\mathrm{\Gamma }_1_k\xi ^k)_j\xi ^j_i(p\mathrm{\Gamma }_1_k\xi ^k)`$
It will be convenient to rewrite the second order terms on the right-hand side of this equation as a divergence. (The first order terms can also be written as a divergence.) The last two terms are

$`_i\xi ^j_j(p\mathrm{\Gamma }_1_k\xi ^k)_j\xi ^j_i(p\mathrm{\Gamma }_1_k\xi ^k)`$ $`=`$ $`_j(_i\xi ^jp\mathrm{\Gamma }_1_k\xi ^k)`$ (I49)
$`_i(_j\xi ^jp\mathrm{\Gamma }_1_k\xi ^k)`$
$`=`$ $`_j(p\mathrm{\Gamma }_1\mathrm{\Theta }^j{}_{i}{}^{})_i(p\mathrm{\Gamma }_1\mathrm{\Theta }).`$ (I50)
Also, by Eq. (I43) we have

$`(\mathrm{\Xi }^k{}_{i}{}^{}\mathrm{\Theta }^k{}_{i}{}^{})_kp`$ $`=`$ $`_k[p(\mathrm{\Xi }^k{}_{i}{}^{}\mathrm{\Theta }^k{}_{i}{}^{})]p_k(\mathrm{\Xi }^k{}_{i}{}^{}\mathrm{\Theta }^k{}_{i}{}^{})`$ (I51)
$`=`$ $`_k[p(\mathrm{\Xi }^k{}_{i}{}^{}\mathrm{\Theta }^k{}_{i}{}^{})]{\displaystyle \frac{1}{2}}p_i(\mathrm{\Xi }\mathrm{\Theta }).`$ (I52)
Thus Eq. (I44) becomes

$`\rho \mathrm{\Delta }\left({\displaystyle \frac{1}{\rho }}_ip\right)`$ $`=`$ $`_j\xi ^j_ip_i\xi ^j_jp_i(p\mathrm{\Gamma }_1_j\xi ^j)`$ (I53)
$`+_j\{p(\mathrm{\Gamma }_11)\mathrm{\Theta }^j{}_{i}{}^{}+{\displaystyle \frac{1}{2}}\delta ^j{}_{i}{}^{}p[(\mathrm{\Gamma }_11)^2+\rho {\displaystyle \frac{\mathrm{\Gamma }_1}{\rho }}]\mathrm{\Theta }`$ (I55)
$`+p\mathrm{\Xi }^j{}_{i}{}^{}+{\displaystyle \frac{1}{2}}\delta ^j{}_{i}{}^{}p(\mathrm{\Gamma }_11)\mathrm{\Xi }\}.`$
Multiplying Eq. (I1) by $`\rho `$, and using eqs. (I29), (I53), and (I37), we obtain the equation of motion for the perturbation to second order:

$`\rho {\displaystyle \frac{d^2\xi _i}{dt^2}}`$ $`=`$ $`_j\xi ^j_ip+_i\xi ^j_jp+_i(p\mathrm{\Gamma }_1_j\xi ^j)\rho _i\delta ^{(1)}\varphi \rho \xi ^j_j_i\varphi `$ (I59)
$`_j\{p(\mathrm{\Gamma }_11)\mathrm{\Theta }^j{}_{i}{}^{}+{\displaystyle \frac{1}{2}}\delta ^j{}_{i}{}^{}p[(\mathrm{\Gamma }_11)^2+\rho {\displaystyle \frac{\mathrm{\Gamma }_1}{\rho }}]\mathrm{\Theta }`$
$`+p\mathrm{\Xi }^j{}_{i}{}^{}+{\displaystyle \frac{1}{2}}\delta ^j{}_{i}{}^{}p(\mathrm{\Gamma }_11)\mathrm{\Xi }\}\rho _i\delta ^{(2)}\varphi `$
$`\rho \xi ^j_j_i\delta ^{(1)}\varphi {\displaystyle \frac{1}{2}}\rho \xi ^j\xi ^k_j_k_i\varphi .`$
Here we have split the quantity $`\delta \varphi `$ defined by eqs. (I34) and (I37) into first and second order pieces $`\delta ^{(1)}\varphi `$ and $`\delta ^{(2)}\varphi `$. Accordingly the first line on the right-hand side of Eq. (I59) gives the first order terms in the equation of motion, while the remaining lines give the second order terms.

Finally, we note that the above derivation was carried out in the inertial frame. To obtain the equation of motion in the rotating frame, we replace the inertial-frame Euler equation (I1) at the start of the calculation by the rotating-frame Euler equation (2). This involves adding to the left hand side the term

$$2𝛀\times 𝐮,$$
(I60)
and replacing the background potential $`\varphi `$ by $`\varphi +\varphi _{\mathrm{rot}}`$, where the centrifugal potential $`\varphi _{\mathrm{rot}}`$ is given by Eq. (11). Taking the variation of the term (I60) using Eq. (I12) yields

$$2𝛀\times \frac{d𝝃}{dt},$$
(I61)
which should be added to the left hand side of Eq. (I59), and substituting $`\varphi \varphi +\varphi _{\mathrm{rot}}`$ on the right hand side results in the additional term $`\rho \xi ^j_j_i\varphi _{\mathrm{rot}}`$ on the right hand side. In the application of this paper, the background velocity in the rotating frame vanishes, and so we can replace $`d𝝃/dt`$ with $`𝝃/t`$ and $`d^2𝝃/dt^2`$ with $`^2𝝃/t^2`$. With these modifications, Eq. (I59) reduces to the result quoted in the body of the paper, given by Eqs. (8), (10), (112), (113), (119) and (120).

## J Variational Derivation of Second Order Equations of Motion

Derivations of the hydrodynamical equations directly from a variational principle have a long history , and have been applied to perturbation theory for nonrotating stars (e.g. Ref. ). Hydrodynamical equations can be derived from a Lagrangian viewpoint by extremizing the action

$$S=𝑑t(t)=𝑑t[T(t)V(t)],$$
(J1)
where the Lagrangian is the difference between kinetic and potential energies for the fluid. Included in the potential energy are both the gravitational potential (both the self potential and external potentials) and the thermal energy. The specific variational principle we use here is the following . Label fluid elements by their positions $`𝐲_0`$ at some initial time $`t=t_0`$, and let the position of fluid element $`𝐲_0`$ at time $`t`$ be $`𝐲(𝐲_0,t).`$ Let $`\rho _0(𝐲_0)`$ be the mass density at time $`t_0`$, then the element of mass in the fluid is $`dM=\rho _0d^3y_0=\rho d^3y`$ where $`\rho =\rho _0/J`$ and $`J`$ is the Jacobian

$$J=\mathrm{det}\left(\frac{y_i}{y_{0j}}\right)$$
(J2)
Then the kinetic energy is

$$T=\frac{1}{2}𝑑M\left(\frac{𝐲}{t}\right)_{𝐲_0}^2,$$
(J3)
and the potential is

$$V=𝑑M(\rho )\frac{G}{2}d^3y_0d^3y_0^{}\rho _0(𝐲_0)\rho _0(𝐲_0^{})\frac{1}{\left|𝐲(𝐲_0,t)𝐲(𝐲_0^{},t)\right|}.$$
(J4)
Here $`\rho `$ is understood to be the functional of $`𝐲(𝐲_0,t)`$ given by $`\rho =\rho _0/J`$, and $``$ is the internal energy per unit mass of the fluid, related to the equation of state by $`d/d\rho =p(\rho )/\rho ^2`$. With these definitions, varying the action (J1) with respect to $`𝐲(𝐲_0,t)`$ gives the Euler equation .

We can generalize the above variational principle slightly by allowing the energy per unit mass to depend on additional hydrodynamic variables such as specific entropy or composition. We denote the set of such variables by the vector $`\stackrel{}{\mu }`$, so that

$$=(\rho ,\stackrel{}{\mu }).$$
(J5)
In this case, the action principle (J1) does not determine the evolution of the variables $`\stackrel{}{\mu }`$. However, if we restrict attention to situations in which the variables $`\stackrel{}{\mu }`$ are determined by the function $`𝐲(𝐲_0,t)`$, then the action (J1) does determine the motion of the fluid. An example of such a situation is where each fluid element conserves its values of $`\stackrel{}{\mu }`$; see below.

We now apply this principle to perturbations of a uniformly rotating star. We denote the background quantities by $`\rho _0`$, $`p_0`$ and $`\varphi _0`$; these were denote by $`\rho `$, $`p`$ and $`\varphi `$ in the body of the paper. We can reformulate the variational principle (J1) by transforming to the co-rotating frame. Define $`𝐲_0=𝐔(t)𝐱_0`$ and $`𝐲=𝐔(t)𝐱`$, where $`𝐔(t)`$ is the appropriate orthogonal rotation matrix. The form of the piece (J4) of the action does not change under this transformation, so we merely have to transform the kinetic energy (J3). In the rotating frame, the position at time $`t`$ of a mass element with comoving coordinate $`𝐱_0`$ is

$$𝐱=𝐱(𝐱_0,t)𝐱_0+𝝃(𝐱_0,t).$$
(J6)
\[Note that the $`𝐱_0`$ of this appendix was denoted $`𝐱`$ in the body of the paper.\] The position as seen in the inertial frame in which the star is at rest is

$$y_i=U_{ij}(t)[x_{0j}+\xi _j(𝐱_0,t)].$$
(J7)
The inertial-frame fluid velocity is found by differentiating this equation with respect to time:

$$u_i=\frac{dU_{ij}}{dt}[x_{0j}+\xi _j(𝐱_0,t)]+U_{ij}(t)\frac{\xi _j(𝐱_0,t)}{t}.$$
(J8)
Using the identity $`dU_{ij}/dt=ϵ_{ipq}\mathrm{\Omega }_pU_{qj}(t)`$, we find

$$𝐔^1(t)𝐮=𝛀\times \left[𝐱_0+𝝃(𝐱_0,t)\right]+\frac{𝝃(𝐱_0,t)}{t},$$
(J9)
where we have used $`\mathrm{\Omega }_m=U_{pm}(t)\mathrm{\Omega }_p`$. This yields

$$|𝐮|^2=|𝛀\times [𝐱_0+𝝃(𝐱_0,t)]|^22[𝐱_0+𝝃(𝐱_0,t)]\left(𝛀\times \frac{𝝃(𝐱_0,t)}{t}\right)+\left|\frac{𝝃}{t}\right|^2.$$
(J10)
To find the equations of motion for the perturbations of the star, we can dispense with all terms in the Lagrangian that are either zeroth or first order in $`𝝃`$ and its time derivatives‡§‡§‡§Note that there is a term proportional to $`\dot{𝝃}`$ in $`|𝐮|^2`$, but varying the action with respect to this term does not contribute to the equations of motion. This is because the result of varying this term is a total time derivative. We will denote with subscripts 2 quantities from which the zeroth and first order pieces in $`𝝃`$ have been subtracted. Thus, the relevant piece of the kinetic energy is

$$T_2=𝑑M\left[\frac{1}{2}\left|\frac{𝝃}{t}\right|^2+\frac{1}{2}[\mathrm{\Omega }^2|𝝃|^2(𝛀𝝃)^2]𝝃\left(𝛀\times \frac{𝝃(𝐱_0,t)}{t}\right)\right].$$
(J11)
We also note that we can define the canonical momentum

$$𝝅=\frac{𝒯_2}{(𝝃/t)}=\frac{𝝃}{t}+𝛀\times 𝝃,$$
(J12)
where $`𝒯_2`$ is the integrand in Eq. (J11), and that

$$T_2𝝅\frac{𝝃}{t}=\frac{|𝝅|^2}{2}+𝝅(𝛀\times 𝝃)=\frac{|𝝅𝛀\times 𝝃|^2}{2}+\frac{|𝛀\times 𝝃|^2}{2}.$$
(J13)
A Hamiltonian for the system is thus given by

$$H_2=\frac{1}{2}𝑑M\left[|𝝅𝛀\times 𝝃|^2|𝛀\times 𝝃|^2\right]+V_2(t),$$
(J14)
where $`V_2`$ is the potential (J4) with zeroth order and first order terms dropped, which only depends on $`𝝃`$, $`\mathbf{}_0𝝃`$ and the background star, but not on $`𝝅`$. Hamilton’s equations then yield the equations of motion

$`\dot{𝝃}`$ $`=`$ $`𝝅𝛀\times 𝝃`$ (J15)
$`\dot{𝝅}`$ $`=`$ $`{\displaystyle \frac{𝐟}{\rho _0}}𝛀\times 𝝅,`$ (J16)
where the accelerations $`𝐟`$ is determined from the variation of the potential $`V_2`$. In the linear approximation, these results are equivalent to the first six equations in Appendix A, with no external force. We can combine Eqs. (J16) into the single equation

$$\frac{^2𝝃}{t^2}+2𝛀\times \frac{𝝃}{t}=\frac{𝐟}{\rho _0}𝛀\times (𝛀\times 𝝃).$$
(J17)
It remains to compute the variations of the potential, $`V_2`$, which contains two parts, internal energy and gravitational energy. The internal energy is

$$U=𝑑M(v,\stackrel{}{\mu })$$
(J18)
where $`v`$ is the volume per unit mass. (We express $``$ here as a function of $`v`$ instead of $`\rho =1/v`$). As explained above, $`\stackrel{}{\mu }`$ is a set of variables such as the specific entropy and composition of the fluid. The volume per unit mass in the fluid is related to its value in the unperturbed state, $`v_0`$, by

$$v=J(𝐱,𝐱_0)v_0,$$
(J19)
where $`J(𝐱,𝐱_0)=\mathrm{det}(1+L)`$, with $`L_j^i=\delta _{ij}+\xi ^i/x_{0,j}`$ \[cf. Eqs. (I15) and (I16)\]. The Jacobian $`J`$ is given by the exact formula

$$J(𝐱,𝐱_0)=1+𝝃+\frac{1}{2}\left(\mathrm{\Theta }\mathrm{\Xi }\right)+\frac{1}{6}(𝝃)\left(\mathrm{\Theta }3\mathrm{\Xi }\right)+\frac{1}{3}\chi $$
(J20)
where $`\mathrm{\Theta }`$, $`\mathrm{\Xi }`$ and $`\chi `$ are the traces of the tensors defined in Eq. (117). This formula is essentially just the characteristic equation for the $`3\times 3`$ matrix $`\xi _{;j}^i`$.

To go further, we need to specify how $`\stackrel{}{\mu }`$ changes in the presence of fluid displacements. We assume that for each fluid element $`\stackrel{}{\mu }`$ remains frozen at its background value $`\stackrel{}{\mu }_0`$ ‡¶‡¶‡¶More generally, we could assume that $`\stackrel{}{\mu }`$ changes only in response to gradients in $`𝝃`$, and, since the various components of $`\stackrel{}{\mu }`$ are scalars, presumably would depend only on scalar quantities that can be constructed from $`\mathbf{}_0𝝃`$.. This means that we focus on perturbations that preserve for each fluid element the composition, specific entropy, etc. of the unperturbed star. We are essentially assuming that the timescales for the fluid to relax to local composition are relatively long.‡∥‡∥‡∥To do better, we could introduce equations of heat transport, composition and so on.

Although it is easy to subtract the contributions to $`U`$ that are zeroth and first order in $`𝝃`$, for computing its variation it is easier to work with $`U`$ rather than $`U_2`$, and subtract off the displacement field through $`J(𝐱,𝐱_0)`$ only. The result is

$`\delta U`$ $`=`$ $`{\displaystyle 𝑑Mv_0p[v_0J(𝐱,𝐱_0),\stackrel{}{\mu }_0]\delta J(𝐱,𝐱_0)}`$ (J21)
$`=`$ $`{\displaystyle \mathrm{d}^3𝐱_0p[v_0J(𝐱,𝐱_0),\stackrel{}{\mu }_0]\delta J(𝐱,𝐱_0)},`$ (J22)
where $`p(v)=(v,\stackrel{}{\mu })/v`$ is the pressure. To simplify further, we substitute $`\delta J(𝐱,𝐱_0)=J(𝐱,𝐱_0)(1+L)_{ji}^1\delta \xi _i/x_{0j}`$, and then integrate by parts, and impose the surface boundary condition $`p=0`$; after subtracting the pressure gradient force in the background star, the result is

$`\delta U_2`$ $`=`$ $`{\displaystyle \mathrm{d}^3𝐱_0\delta \xi _i\frac{}{x_{0j}}\left[p[v_0J(𝐱,𝐱_0),\stackrel{}{\mu }_0]J(𝐱,𝐱_0)(1+L)_{ji}^1\delta _{ij}p(v_0,\stackrel{}{\mu }_0)\right]}`$ (J23)
Note the exact relation

$$J(1+L)_{ji}^1=[1+𝝃+\frac{1}{2}(\mathrm{\Theta }\mathrm{\Xi })]\delta _{ji}\xi _{j;i}+\mathrm{\Xi }_{ji}\mathrm{\Theta }_{ji}.$$
(J24)
(For a derivation, see e.g. Ref. .) Expanding the quantity in square brackets in Eq. (J23) to second order, we find

$`p[v_0J(𝐱,𝐱_0),\stackrel{}{\mu }_0]J(𝐱,𝐱_0)(1+L)_{ji}^1\delta _{ij}p(v_0,\stackrel{}{\mu }_0)=p_0(v_0,\stackrel{}{\mu }_0)\{\delta _{ij}(1\mathrm{\Gamma }_1)\mathbf{}_0𝝃{\displaystyle \frac{\xi _j}{x_{0i}}}`$ (J25)
$`+{\displaystyle \frac{\delta _{ij}}{2}}[\mathrm{\Theta }((\mathrm{\Gamma }_11)^2+\rho _0{\displaystyle \frac{\mathrm{\Gamma }_2}{\rho _0}})+(\mathrm{\Gamma }_11)\mathrm{\Xi }]+(\mathrm{\Gamma }_11)\mathrm{\Theta }_i^j+\mathrm{\Xi }_i^j+\mathrm{}\},`$ (J26)
which gives the same result, to second order, as Eq. (I59). Here, we have defined a generalized version of the adiabatic index of the perturbation,

$$\mathrm{\Gamma }_1=\left(\frac{\mathrm{ln}p_0(v_0,\stackrel{}{\mu }_0)}{\mathrm{ln}v_0}\right)_{\stackrel{}{\mu }_0}=\left(\frac{\mathrm{ln}p_0(\rho _0,\stackrel{}{\mu }_0)}{\mathrm{ln}\rho _0}\right)_{\stackrel{}{\mu }_0},$$
(J27)
which also applies in a zero entropy star, for perturbations whose composition, and possibly other characteristics, remain fixed at their unperturbed values. It is straightforward to extend the expansion to any desired order, and also to expand the internal energy per unit mass, $`(v_0J,\stackrel{}{\mu }_0)`$, to any desired order.

Finally, we need to consider the gravitational potential energy,

$$V_G=\frac{G}{2}\frac{\mathrm{d}M\mathrm{d}M^{}}{|𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)|}.$$
(J28)
Varying with respect to the displacement field gives

$$\delta V_G=\frac{G}{2}dMdM^{}[\delta 𝝃(𝐱_0,t)\delta 𝝃(𝐱_0^{},t)]_i\frac{[𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)]_i}{|(𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)|^3};$$
(J29)
by switching primed and unprimed variables, and subtracting off the zeroth order gravitational force, we find

$$\delta V_{G}^{}{}_{2}{}^{}=GdM\delta 𝝃(𝐱_0,t)dM^{}\left[\frac{𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)}{|(𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)|^3}\frac{𝐱_0𝐱_0^{}}{|𝐱_0𝐱_0^{}|^3}\right].$$
(J30)
Gathering results, we see that the acceleration in Eq.(J16), computed to any desired order in the amplitude of the Lagrangian perturbation field, $`𝝃(𝐱_0,t)`$, is $`𝐟/\rho _0=𝐟_P/\rho _0+𝐟_G/\rho _0`$, where

$`𝐟_P^i`$ $`=`$ $`{\displaystyle \frac{}{x_{0,j}}}\left[p[v_0J(𝐱,𝐱_0),\stackrel{}{\mu }_0]J(𝐱,𝐱_0)(1+L)_{ji}^1\delta _{ij}p(v_0,\stackrel{}{\mu }_0)\right]`$ (J31)
$`{\displaystyle \frac{𝐟_G}{\rho _0}}`$ $`=`$ $`G{\displaystyle dM^{}\left[\frac{𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)}{|(𝐱_0𝐱_0^{}+𝝃(𝐱_0,t)𝝃(𝐱_0^{},t)|^3}\frac{𝐱_0𝐱_0^{}}{|𝐱_0𝐱_0^{}|^3}\right]}.`$ (J32)
Equations (J26) and (J31) allow us to calculate $`𝐟_P`$ to second order in the perturbation; the result is Eq. (120). Similarly, let us expand $`𝐟_G=𝐟_G^{(1)}+𝐟_G^{(2)}+\mathrm{}`$; the contributions from various orders can be found by expanding Eq. (J32). To first order, the result of the calculation is that ‡\**‡\**‡\**An integration by parts is required to get this form of the gravitational force. The outer boundary has been taken to be outside the star, so if the mass density goes to zero discontinuously there, a surface term must be included.

$$𝐅_i^{(1)}=\rho _0𝝃_j(𝐱_0,t)\mathbf{}_{0,i}\mathbf{}_{0,j}\varphi _0(𝐱_0)\rho _0\mathbf{}_{0,i}\delta ^{(1)}\varphi (𝐱_0,t),$$
(J33)
where

$$\delta ^{(1)}\varphi (𝐱_0,t)=G\mathrm{d}^3𝐱_0^{}\frac{\mathbf{}_{0,j}^{}[\rho _0(𝐱_0^{})𝝃_j(𝐱_0^{},t)]}{|𝐱_0𝐱_0^{}|}.$$
(J34)
The second order contribution to the gravitational force is (compare Eq.\[I59\])

$$𝐅_i^{(2)}=\frac{\rho _0𝝃_j𝝃_k}{2}\mathbf{}_{0,j}\mathbf{}_{0,k}\mathbf{}_{0,i}\varphi _0\rho _0𝝃_j\mathbf{}_{0,j}\mathbf{}_{0,i}\delta \varphi \rho _0\mathbf{}_{0,i}\delta ^{(2)}\varphi ,$$
(J35)
where

$$\delta ^{(2)}\varphi \frac{G}{2}\mathrm{d}^3𝐱_0^{}\rho _0(𝐱_0)𝝃_j(𝐱_0^{},t)𝝃_k(𝐱_0^{},t)\mathbf{}_{0,j}^{}\mathbf{}_{0,k}^{}|𝐱_0𝐱_0^{}|^1.$$
(J36)
It is possible to integrate by parts to change the form of $`\delta ^{(2)}\varphi `$ to show it coincides with the expression (115).

## K Energy and Angular Momentum of stellar perturbations

In Sec. K 1 of this appendix we give general formulae for the physical energy $`E_{\mathrm{phys}}`$ and physical ($`z`$-component of) angular momentum $`J_{\mathrm{phys}}`$ of a perturbation, and also the rotating-frame canonical energy $`E_{\mathrm{can},\mathrm{rot}}`$. General expressions for these quantities have been derived by Friedman and Schutz , in a context more general than that considered here (uniform rotation of the background star). However, the results of Friedman and Schutz are expressed in terms of an inertial-frame Lagrangian displacement $`𝝃_{\mathrm{in}}`$. For a uniformly rotating background star, it is more natural and simpler to use instead the rotating-frame Lagrangian displacement $`𝝃`$ as we do in this paper. The expressions in Sec. K 1 below can be obtained by specializing and translating the results of Friedman and Schutz, or, as we outline below, can also be easily obtained directly.

We also derive, in Sec. K 2, expressions for these quantities in terms of the mode coefficients $`c_{A\sigma }`$, and in Sec. K 3 compute the physical energy deposited in a star by an externally applied acceleration $`𝐚_{\mathrm{ext}}`$ to second order in $`𝐚_{\mathrm{ext}}`$.

Friedman and Schutz  derived several general properties of the quantities $`E_{\mathrm{phys}}`$, $`J_{\mathrm{phys}}`$ and $`E_{\mathrm{can},\mathrm{rot}}`$. First, the physical energy $`E_{\mathrm{phys}}[𝜻]`$ of a perturbation $`𝜻`$ ‡††‡††‡††This is defined to be the difference between the energy of the perturbed configuration and the energy of the background, unperturbed configuration. has a piece that is linear in $`𝜻`$ as well as pieces that are quadratic and higher order in $`𝜻`$:

$$E_{\mathrm{phys}}[𝜻]=E_{\mathrm{phys}}^{(1)}[𝜻]+E_{\mathrm{phys}}^{(2)}[𝜻,𝜻]+O(𝜻^3),$$
(K1)
and similarly for the physical angular momentum:

$$J_{\mathrm{phys}}[𝜻]=J_{\mathrm{phys}}^{(1)}[𝜻]+J_{\mathrm{phys}}^{(2)}[𝜻,𝜻].$$
(K2)
Second, one has the identity

$$E_{\mathrm{phys}}^{(1)}[𝜻]=\mathrm{\Omega }J_{\mathrm{phys}}^{(1)}[𝜻],$$
(K3)
and so the physical energy in the rotating frame, defined as

$$E_{\mathrm{phys},\mathrm{rot}}E_{\mathrm{phys}}\mathrm{\Omega }J_{\mathrm{phys}}$$
(K4)
has no linear term in $`𝜻`$:

$$E_{\mathrm{phys},\mathrm{rot}}[𝜻]=E_{\mathrm{phys},\mathrm{rot}}^{(2)}[𝜻,𝜻]+O(𝜻^3)$$
(K5)
where

$$E_{\mathrm{phys},\mathrm{rot}}^{(2)}[𝜻,𝜻]=E_{\mathrm{phys}}^{(2)}[𝜻,𝜻]\mathrm{\Omega }J_{\mathrm{phys}}^{(2)}[𝜻,𝜻].$$
(K6)
Third, this quantity coincides with the so-called rotating frame canonical energy $`E_{\mathrm{can},\mathrm{rot}}[𝜻]`$ defined in Ref. , that is, the conserved quantity related to the time translation symmetry in the rotating frame by Noether’s theorem.

### 1 Explicit expressions in terms of $`𝜻`$

The physical energy of a perturbation can be computed using the procedure of Appendix J above, and adding the kinetic and potential energies $`T(t)`$ and $`V(t)`$. To linear order, the result is

$`E_{\mathrm{phys}}^{(1)}[𝜻]`$ $`=`$ $`{\displaystyle d^3x\rho _0(𝐱)(𝛀\times 𝐱)\left[\dot{𝝃}+2𝛀\times 𝝃\right]}`$ (K7)
$`=`$ $`𝝃_r,\dot{𝝃}+𝐁𝝃=𝝃_r,𝝅+𝐁𝝃/2,`$ (K8)
where

$$𝝃_r(𝐱)=𝛀\times 𝐱$$
(K9)
is the mode function corresponding to uniform rotation and we have used Eq. (A2). The second order piece is

$`E_{\mathrm{phys}}^{(2)}[𝜻,𝜻]`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{𝝃},\dot{𝝃}+{\displaystyle \frac{1}{2}}𝝃,𝐂𝝃+{\displaystyle \frac{1}{2}}\dot{𝝃},𝐁𝝃+{\displaystyle \frac{1}{4}}𝐁𝝃,𝐁𝝃`$ (K10)
$`=`$ $`{\displaystyle \frac{1}{2}}𝝅,𝝅+{\displaystyle \frac{1}{2}}𝝃,𝐂𝝃+{\displaystyle \frac{1}{8}}𝐁𝝃,𝐁𝝃.`$ (K11)
Here we have used the fact that $`𝝃`$ and $`\dot{𝝃}`$ are real to eliminate terms such as $`\dot{𝝃},𝐁\dot{𝝃}`$. The physical angular momentum is given by

$$𝐉_{\mathrm{phys}}=𝑑M𝐲\times 𝐮(𝐲).$$
(K12)
Here we are using the notation of Appendix J, where $`𝐲`$ and $`𝐮`$ are the inertial-frame location and velocity of a fluid element. Using Eqs. (J7) and (J9) in Eq. (K12) gives

$$𝐔(t)^1𝐉_{\mathrm{phys}}=d^3x\rho _0(𝐱)(𝐱+𝝃)\times \left[𝛀\times 𝐱+𝛀\times 𝝃+\dot{𝝃}\right],$$
(K13)
and the $`z`$-component is given by

$`\mathrm{\Omega }J_{\mathrm{phys}}[𝝃]`$ $`=`$ $`𝛀𝐉_{\mathrm{phys}}=𝛀𝐔(t)^1𝐉_{\mathrm{phys}}`$ (K14)
$`=`$ $`{\displaystyle d^3x\rho _0(𝐱)\left[𝛀\times 𝐱+𝛀\times 𝝃\right]\left[𝛀\times 𝐱+𝛀\times 𝝃+\dot{𝝃}\right]}.`$ (K15)
Expanding this in powers of $`𝝃`$ one obtains at linear order using Eq. (K8) the identity (K3). At second order we obtain

$`\mathrm{\Omega }J_{\mathrm{phys}}^{(2)}[𝜻,𝜻]`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{𝝃},𝐁𝝃+{\displaystyle \frac{1}{4}}𝐁𝝃,𝐁𝝃`$ (K16)
$`=`$ $`{\displaystyle \frac{1}{2}}𝝅,𝐁𝝃.`$ (K17)
The physical energy in the rotating frame is therefore, from Eqs. (K6), (K11), and (K17)

$`E_{\mathrm{phys},\mathrm{rot}}[𝜻]`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{𝝃},\dot{𝝃}+{\displaystyle \frac{1}{2}}𝝃,𝐂𝝃+(𝜻^3)`$ (K18)
$`=`$ $`{\displaystyle \frac{1}{2}}𝝅,𝝅+{\displaystyle \frac{1}{2}}𝝃,𝐂𝝃+{\displaystyle \frac{1}{8}}𝐁𝝃,𝐁𝝃{\displaystyle \frac{1}{2}}𝝅,𝐁𝝃+O(𝜻^3).`$ (K19)
The expression (K18) agrees with the expression Eq. (43) of Ref.  for rotating-frame canonical energy, as expected ‡‡‡‡‡‡‡‡‡Equation (43) of Ref.  is actually an expression for inertial-frame canonical energy. The operators which appear in that expression, which we will denote by $`𝐀_{\mathrm{in}}`$, $`𝐁_{\mathrm{in}}`$ and $`𝐂_{\mathrm{in}}`$, are defined such that the linearized equation of motion is $`𝐀_{\mathrm{in}}\ddot{𝝃}_{\mathrm{in}}+𝐁_{\mathrm{in}}\dot{𝝃}+𝐂_{\mathrm{in}}𝝃=0`$, where $`𝝃_{\mathrm{in}}`$ is the inertial-frame Lagrangian displacement. However, one can obtain an expression for the rotating-frame canonical energy, as defined by Eq. (64) of Ref. , by replacing in Eq. (43) the inertial frame quantities $`𝐀_{\mathrm{in}}`$, $`𝐁_{\mathrm{in}}`$, $`𝐂_{\mathrm{in}}`$ and $`𝝃_{\mathrm{in}}`$ by their rotating-frame counterparts used in this paper, $`𝐀`$, $`𝐁`$, $`𝐂`$, and $`𝝃`$. Note that in our notational convention the operator $`𝐀`$ is the identity operator..

### 2 Expressions in terms of mode coefficients

From the above discussion we can write the total energy of the perturbation, to quadratic order, as the sum of three terms, a quadratic rotating-frame energy (physical or canonical), a linear angular momentum term, and a quadratic angular momentum term:

$$E_{\mathrm{phys}}[𝜻]=E_{\mathrm{phys},\mathrm{rot}}^{(2)}[𝜻,𝜻]+\mathrm{\Omega }J_{\mathrm{phys}}^{(1)}[𝜻]+\mathrm{\Omega }J_{\mathrm{phys}}^{(2)}[𝜻,𝜻]+O(𝜻^3).$$
(K20)
Of these three terms, the first is diagonalized by the mode basis $`𝜻_{A\sigma }`$, the second can be expressed as the coefficient $`c_{A\sigma }`$ of a single Jordan chain mode (corresponding to uniform rotation), and the third is not diagonalized by the mode basis $`𝜻_{A\sigma }`$ and contains cross terms between different modes.

#### a Rotating-frame energy

Consider first the rotating-frame energy. From Eqs. (A7), (A10), (A68) and (K19) we can write this as

$$E_{\mathrm{phys},\mathrm{rot}}^{(2)}[𝜻,𝜻]=\frac{i}{2}𝜻,𝐌𝐓𝜻.$$
(K21)
Inserting the mode expansion (A26) and using the definition (A22) of a right Jordan chain and the definition (A72) of the matrix $`_{A\sigma ,B\tau }`$ gives

$$E_{\mathrm{phys},\mathrm{rot}}^{(2)}[𝜻,𝜻]=\underset{A\sigma ,B\tau }{}c_{A\sigma }^{}c_{B\tau }\left[\frac{\omega _B}{2}_{A\sigma ,B\tau }+\frac{i}{2}_{A\sigma ,B(\tau 1)}\right].$$
(K22)
Since the energy is real and $``$ is Hermitian we can rewrite this as

$$E_{\mathrm{phys},\mathrm{rot}}^{(2)}[𝜻,𝜻]=\frac{1}{4}\underset{A\sigma ,B\tau }{}c_{A\sigma }^{}c_{B\tau }\left[(\omega _A^{}+\omega _B)_{A\sigma ,B\tau }+i_{A\sigma ,B(\tau 1)}i_{A(\sigma 1),B\tau }\right].$$
(K23)
Since the matrix $``$ is block diagonal (see Sec. A 5 above), this rotating-frame energy can be written as a sum of three contributions, from (i) real-frequency modes, (ii) complex-frequency (unstable) modes, and (ii) non-trivial Jordan chain modes. The term corresponding to real frequencies is

$$\frac{1}{2}\underset{A}{}\epsilon _A|c_A|^2,$$
(K24)
where we have assumed that one is using a basis which diagonalizes the real-frequencies portion of the matrix $`_{AB}`$ (cf. Sec. II C 3 above), and $`\epsilon _A`$ is the mode rotating-frame-energy at unit amplitude, given from Eqs. (45) and (50) by

$$\epsilon _A=\omega _A_{AA}=\omega _Ab_A.$$
(K25)
Using the notations (57) we can rewrite the energy (K24) as

$$\underset{\alpha }{}|c_\alpha |^2\omega _\alpha b_\alpha =\underset{\alpha }{}|c_\alpha |^2\epsilon _\alpha .$$
(K26)
To write down the portion of the expression (K23) corresponding to complex frequencies, we need to introduce some more notation. The modes come in pairs $`(𝝃,\omega )`$ and $`(g_{}𝝃^{},\omega ^{})`$; see Secs. A 4 and A 5 above. We write $`A=(\widehat{a}\delta k)`$, where $`\widehat{a}`$ and $`\delta `$ label the frequency $`\omega _{\widehat{a}\delta k}=\omega _{\widehat{a}\delta }`$, and $`k`$ labels the different (degenerate) eigenvectors associated with that frequency. The index $`\delta `$ can take on the values $`1`$ and $`1`$ such that

$`\omega _{\widehat{a}1}`$ $`=`$ $`\omega _{\widehat{a}}`$ (K27)
$`\omega _{\widehat{a}(1)}`$ $`=`$ $`\omega _{\widehat{a}}^{}`$ (K28)
$`𝝃_{\widehat{a}1k}`$ $`=`$ $`𝝃_{\widehat{a}k}`$ (K29)
$`𝝃_{\widehat{a}(1)k}`$ $`=`$ $`(g_{}𝝃_{\widehat{a}k})^{}.`$ (K30)
It is possible to choose the mode basis such that the matrix $``$ takes the form

$`_{\widehat{a}\delta k,\widehat{b}ϵl}`$ $`=`$ $`\delta _{\widehat{a}\widehat{b}}\delta _{\delta ,ϵ}\delta _{kl}\left[\delta _{\delta ,1}𝒟_{\widehat{a}k}+\delta _{\delta ,1}𝒟_{\widehat{a}k}^{}\right];`$ (K31)
cf. Eq. (A92) above. Inserting this into the expression (K23) gives for the contribution to $`E_{\mathrm{phys},\mathrm{rot}}^{(2)}`$ from complex-frequency modes

$$\underset{\widehat{a}}{}\mathrm{Re}\left[\underset{k}{}\omega _{\widehat{a}1}c_{\widehat{a}1k}^{}c_{\widehat{a}(1)k}𝒟_{\widehat{a}k}\right].$$
(K32)
We see that if only a single mode is excited, the corresponding rotating-frame energy vanishes, in agreement with the result of Friedman and Schutz . However, if a mode $`(𝝃,\omega )`$ and also its conjugate pair $`(g_{}𝝃^{},\omega ^{})`$ are both excited, then there is a nonzero contribution to the rotating-frame energy.

Finally, there is a contribution to $`E_{\mathrm{phys},\mathrm{rot}}^{(2)}`$ from non-trivial Jordan chain modes. For simplicity, we restrict attention here to the case of zero-frequency Jordan chains of length one. For this case the expression (K23) reduces to

$$\frac{i}{2}\underset{AB}{}_{A1,B0}c_{A1}^{}c_{B1},$$
(K33)
where we have used the identities (A120)–(A122). As discussed in Sec. A 5 a above, we can always choose the mode basis to diagonalize the matrix

$`\beta _{AB}`$ $`=`$ $`i_{A1,B0}=i𝜻_{A1},𝐌𝜻_{B0}`$ (K34)
$`=`$ $`𝝃_{A0},𝝃_{B0}𝝃_{A1},𝐂𝝃_{B1},`$ (K35)
and the expression (K33) therefore becomes

$$\frac{1}{2}\underset{A}{}\beta _A|c_{A1}|^2,$$
(K36)
where $`\beta _A=\beta _{AA}`$ are the eigenvalues of the matrix $`𝜷`$. Note that there is no rotating-frame energy associated with the $`c_{A0}`$ coefficients.

To summarize, the total rotating-frame energy is given by the sum of the expressions (K26), (K32), and (K36).

#### b Linear piece of angular momentum

The linear piece of the $`z`$-component of the angular momentum (proportional to the linear piece of the physical, inertial-frame energy) can be written as the coefficient of a single Jordan-chain mode. That Jordan-chain mode corresponds to the degree of freedom of uniform spin-up $`\mathrm{\Omega }\mathrm{\Omega }+\mathrm{\Delta }\mathrm{\Omega }`$ of the star, and is a zero-frequency Jordan chain of length one. As explained in Sec. D above, the space of such zero-frequency, length one Jordan-chain modes contains all the differential-rotation modes and is infinitely degenerate. In Sec. A 5 above, we explained how to compute the basis $`𝝌_{A\sigma }`$ of left Jordan chains of this space, starting from a specified basis $`𝜻_{A\sigma }`$ of right Jordan chains. Here, however, it turns out to be more convenient to specify the basis by first choosing the basis $`𝝌_{A\sigma }`$ of left Jordan chains, satisfying Eq. (A23), and by defining $`𝜻_{A\sigma }`$ to be the dual basis defined by Eq. (A24).

The uniform-rotation mode, specified in terms of its left Jordan chain $`𝝌_\sigma `$, $`\sigma =1,2`$, is as follows. Let $`𝝃_r`$ be the uniform-rotation mode function (K9), and define $`𝝉_1`$ to be any solution to the equation

$$𝐂𝝉_1=𝐁𝝃_r.$$
(K37)
Define

$`𝝌_0=\left[\begin{array}{c}𝐁𝝃_r/2\\ 𝝃_r\end{array}\right]`$ (K40)
and

$`𝝌_1=\left[\begin{array}{c}𝝃_r𝐁𝝉_1/2\\ 𝝉_1\end{array}\right].`$ (K43)
Then one can verify using the relation $`𝐂𝝃_r=0`$ that the left Jordan chain relations $`𝐓^{}𝝌_0=0`$ and $`𝐓^{}𝝌_1=𝝌_0`$ are satisfied. The basis $`𝝌_{A\sigma }`$ of left Jordan chain modes can be completed by adding to the above chain any linearly independent set of left Jordan chains describing the differential rotation modes. From Eqs. (A27), (K8) and (K40), the coefficient $`c_\sigma `$ with $`\sigma =1`$ corresponding to the uniform-rotation mode is

$`c_1`$ $`=`$ $`𝝌_0,𝜻=𝝃_r,𝝅+𝐁𝝃/2`$ (K44)
$`=`$ $`E_{\mathrm{phys}}^{(1)}[𝜻]=\mathrm{\Omega }J_{\mathrm{phys}}^{(1)}[𝜻].`$ (K45)
#### c Quadratic piece of angular momentum

From Eq. (K17), the quadratic piece of the angular momentum can be written as

$$\mathrm{\Omega }J_{\mathrm{phys}}^{(2)}[𝜻,𝜻]=𝜻,𝐊𝜻=\underset{A\sigma ,B\tau }{}c_{A\sigma }^{}c_{B\tau }𝒦_{A\sigma ,B\tau },$$
(K46)
where

$`𝐊={\displaystyle \frac{1}{4}}\left[\begin{array}{cc}0& 𝐁\\ 𝐁& 0\end{array}\right]`$ (K49)
and $`𝒦_{A\sigma ,B\tau }=𝜻_{A\sigma },𝐊𝜻_{B\tau }`$. For non-Jordan chain, real-frequency modes we can write the coefficients $`𝒦_{AB}`$, using the orthogonality relation (35) and Eq. (A44), as

$`𝒦_{AB}`$ $`=`$ $`{\displaystyle \frac{1}{4}}(\omega _A+\omega _B)𝝃_A,i𝐁𝝃_B+{\displaystyle \frac{1}{4}}𝐁𝝃_A,𝐁𝝃_B`$ (K50)
$`=`$ $`{\displaystyle \frac{1}{2}}\omega _Ab_A\delta _{AB}{\displaystyle \frac{1}{4}}(\omega _A+\omega _B)^2𝝃_A,𝝃_B+{\displaystyle \frac{1}{4}}𝐁𝝃_A,𝐁𝝃_B.`$ (K51)
For regular modes in a rotating star, in the limit $`\mathrm{\Omega }0`$, the second term in Eq. (K50) scales as $`O(\mathrm{\Omega }^2)`$ and is therefore small compared to the first term which is linear in $`\mathrm{\Omega }`$. Therefore in general the off-diagonal elements of $`𝒦_{AB}`$ are non vanishing; see Eq. (99) above.

### 3 Energy deposited by an external force

Turn now to the question of how to compute the energy deposited in a rotating star by an externally applied acceleration $`𝐚_{\mathrm{ext}}(𝐱,t)`$, to quadratic order in $`𝐚_{\mathrm{ext}}`$. The rotating-frame energy can be easily computed, just as one computes the energy deposited in a non-rotating star. Assuming that there are no unstable modes, one evolves the mode amplitudes using Eq. (46) or Eqs. (A131)–(A132), and one computes the energy from Eqs. (K26) and (K36).

The angular momentum contributions to the total inertial-frame energy (K20) are more difficult to compute. As shown above, the quadratic piece of the angular momentum contains cross terms between different modes. Also, in evaluating the linear term $`E_{\mathrm{phys}}^{(1)}[𝜻]`$ given by Eq. (K45), one needs to solve the equation of motion for the coefficient $`c_1(t)`$ of the uniform-rotation mode to quadratic order in $`𝐚_{\mathrm{ext}}`$, i.e., one needs to solve for the back reaction of all the other stellar modes on the uniform-rotation mode. The linear equation of motion for the coefficient $`c_1(t)`$ is given by combining Eqs. (A13), (A30), (K40), which gives the equation for angular momentum conservation to linear order:

$$\dot{c}_1=𝝃_r,𝐚_{\mathrm{ext}}.$$
(K52)
To quadratic order in $`𝐚_{\mathrm{ext}}`$ (or equivalently in $`𝝃`$), one can obtain the equation of motion for $`c_1(t)`$ by using the following form of the second-order equation of motion

$`\ddot{𝝃}+𝐁\dot{𝝃}+𝐂𝝃=𝐚^{(2)}[𝝃,𝝃]+𝐚_{\mathrm{ext}}(𝐱+𝝃,t)+O(𝝃^3),`$ (K53)
which generalizes Eqs. (7) and (112) above. Replacing the external acceleration on the right hand side of Eq. (K52) with the right hand side of Eq. (K53), and using the identity

$$𝝃_r,𝐚^{(2)}[𝝃,𝝃]=\frac{1}{2}𝐁𝝃,𝐂𝝃$$
(K54)
which is derived below, one obtains the second-order equation of motion for $`c_1`$ in the form

$$\dot{c}_1=𝝃_r,𝐚_{\mathrm{ext}}+𝝃_r,(𝝃\mathbf{})𝐚_{\mathrm{ext}}+\frac{1}{2}𝐁𝝃,𝐂𝝃+O(𝜻^3).$$
(K55)
The third term on the right hand side here describes the excitation of the uniform-rotation mode by the other stellar modes.

It is possible to get around these difficulties in the computation of the angular momentum deposited in the star, by computing the total torque acting on the star, and by not attempting to compute how the angular momentum deposited is distributed between the various stellar modes. The total torque is given by§\*§\*§\*If one takes a time derivative of the expression (K15) for the total angular momentum and uses the equation of motion (K53), one obtains the result (K57) with the two additional terms $`𝝃_r,𝐚^{(2)}[𝝃,𝝃]𝐁𝝃,𝐂𝝃/2.`$ This allows one to derive the identity (K54), which can also be verified directly using the expressions (8), (10), (119) and (120).

$`𝛀\dot{𝐉}`$ $`=`$ $`𝛀{\displaystyle d^3x\rho _0(𝐱)\left[𝐱+𝝃(𝐱)\right]\times 𝐚_{\mathrm{ext}}(𝐱+𝝃)}`$ (K56)
$`=`$ $`𝝃_r,𝐚_{\mathrm{ext}}+𝐁𝝃/2,𝐚_{\mathrm{ext}}+𝝃_r,(𝝃\mathbf{})𝐚_{\mathrm{ext}}+O(𝐚_{\mathrm{ext}}^3,𝝃^3).`$ (K57)
Suppose now that the external acceleration can be expanded as

$$𝐚_{\mathrm{ext}}=\epsilon 𝐚_{\mathrm{ext}}^{(1)}+\epsilon ^2𝐚_{\mathrm{ext}}^{(2)}+O(\epsilon ^3),$$
(K58)
and that the response of the stellar modes is similarly expanded as

$$c_{A\sigma }(t)=\epsilon c_{A\sigma }^{(1)}(t)+\epsilon ^2c_{A\sigma }^{(1)}(t)+O(\epsilon ^3).$$
(K59)
The second-order angular momentum deposited can now be computed as follows. First, solve for the linear response $`c_{A\sigma }^{(1)}(t)`$ of the star. Second, compute the angular momentum deposited due to the first term on the right hand side of Eq. (K57), which is the standard expression for torque. Third, add to this the angular momentum deposited by the remaining two terms in the total torque (K57), which can be written as

$$\epsilon ^2\underset{A\sigma }{}c_{A\sigma }^{(1)}(t)_{A\sigma }(t)+O(\epsilon ^3),$$
(K60)
with

$$_{A\sigma }(t)=d^3x\rho _0(𝐱)𝝃_{A\sigma }(𝐱)\mathbf{}\left[𝛀𝐱\times 𝐚_{\mathrm{ext}}^{(1)}(𝐱)\right].$$
(K61)
This procedure yields a expression for the total angular momentum deposited which is a sum over modes $`(A\sigma )`$. Note, however, that the term $`A\sigma `$ in the expression (K60) is not the angular momentum deposited into the mode $`𝜻_{A\sigma }`$; as noted above the total angular momentum contains cross terms between different modes.