text:
# Saturation of the R-mode Instability

## 1  Introduction

What sets the observed spin rates of neutron stars?

Theoretically, we expect neutron stars can rotate up to $`10^3\mathrm{Hz}`$ without breaking apart (Cook et al. 1994b, a; Fryer & Heger 2000; Heger et al. 2000). However, for the rapidly accreting, weakly magnetic LMXB’s, oscillations seen during type I X-ray bursts (Strohmayer et al. 1996), as well as quasi-periodic oscillations (van der Klis 1998), seem to indicate spin frequencies narrowly clustered near $`300\mathrm{Hz}`$. If LMXB’s are the progenitors of millisecond pulsars, and the timescale over which they should be spun up by accretion is only $`10^7\mathrm{yr}`$ for high accretion rates, why aren’t more stars spun up near $`10^3\mathrm{Hz}`$ over their $`>10^9`$ year lifetime?

Wagoner (1984) proposed that for weakly magnetic neutron stars, the spin up torque due to accretion is balanced by spin down torque from gravitational radiation reaction. There are currently two distinct models to explain the non-axisymmetric deformation of the star producing the radiation. The first mechanism involves mass quadrupole deformations of the neutron star crust (Bildsten 1998; Ushomirsky et al. 2000) while the second involves mass-current quadrupole emission from the r-mode instability (Bildsten 1998; Andersson et al. 1999b, 2000), which will be examined in detail in this paper.

Many young neutron stars associated with supernova remnants also seem to be spinning slowly, in spite of the theoretical expectation (Fryer & Heger 2000; Heger et al. 2000) that typical $`825M_{}`$ progenitors lead to neutron stars rotating with periods of order 111 This result depends sensitively on the poorly understood coupling between the core and envelope of the progenitor. Angular momentum transport mechanisms due to, for instance, weak magnetic fields may decrease the rotation rate of the core prior to collapse. $`1\mathrm{msec}`$. Kaspi & Helfand (2002) cite the following examples for the inferred initial spin period $`P_{\mathrm{init}}`$ and age $`T`$ of some of the fastest rotators: the Crab pulsar with $`P_{\mathrm{init}}=19\mathrm{msec}`$ and $`T=948\mathrm{yr}`$; PSR J0537-6910 in host remnant N157B with $`P_{\mathrm{init}}10\mathrm{msec}`$ and $`T=5000\mathrm{yr}`$; PSR B1951+32 in CTB 80 has $`P_{\mathrm{init}}39\mathrm{msec}`$ and $`T=10^5\mathrm{yr}`$. The Crab is by far the most certain estimate for $`P_{\mathrm{init}}`$, with a known age from the historical supernova and measured braking index. However, Kaspi & Helfand (2002) also note several slow rotators, such as PSR J1811-1925 in G11.2-0.3 with $`P_{\mathrm{init}}=62\mathrm{msec}`$ and age $`T=2000\mathrm{yr}`$.

The apparent discrepancy between the theoretically expected fast rotation rates and the observed slow rotation rates could be reconciled if some mechanism could slow down fast rotators, effectively preventing them from reaching millisecond spin rates. The r-mode instability is a possible mechanism.

This instability was discovered by Andersson (1998) and Friedman & Morsink (1998) showed that all rotating, inviscid stars are unstable because of this general relativistic effect. The instability arises when certain stellar oscillation modes, called Rossby waves (or r-modes), are driven unstable by the emission of gravitational waves. As a result, the rotational energy of the star is converted into both mode energy and gravitational waves, causing the star to spin down. Detailed calculations (Lindblom, Owen & Morsink 1998; Andersson et al. 1999a; Kokkotas & Stergioulas 1999; Lindblom et al. 1999; Bildsten & Ushomirsky 2000; Levin & Ushomirsky 2001; Lockitch & Friedman 1999) show that viscous dissipation by large scale shear, boundary layer shear at the crust-core interface, and modified URCA bulk viscosity are likely insufficient to counter this driving in rapidly rotating neutron stars. However, Lindblom & Owen (2002a) point out an interesting mechanism for bulk viscosity arising from hyperon interactions which may overcome the driving. Mendell (2001) has investigated the effects of magnetic fields on the boundary layer, finding that large fields can significantly increase the damping rate. Lastly, the work of Levin & Ushomirsky (2001) shows that damping from the crust-core boundary layer leads to a double-valued instability curve, which may explain why LMXB spin frequencies are lower than those of the millisecond pulsars.

The instability may be important in two respects. First, r-modes in any neutron star rotating faster than some critical rate will become unstable, causing the star to rapidly spin down. Hence, r-modes may set a maximum rotation rate for neutron stars. Second, the enormous amount of energy radiated in gravitational waves may be detectable by LIGO.

In section 2 we review how nonlinear saturation occurs in the limits of weak and strong driving. We derive formal expressions for the saturation amplitude, which depend on the microphysical details of the nonlinear interaction and damping rates. Section 2.1 contains a review of nonlinear coupling of just three oscillation modes, with emphasis on amplitude saturation by the parametric instability. Section 2.2 reviews amplitude saturation by a continuum of modes in which a well defined inertial range exists. In section 3, we discuss the modes present in rapidly rotating neutron stars, arguing that the buoyancy and elastic restoring forces are weak compared to the Coriolis force. We compute WKB inertial eigenmodes in section 3.2. The nonlinear coupling coefficients are computed in section 4, and damping rates in section 5. The saturation amplitude for the discrete case is discussed in section 6, and the continuum case in section 7. Neutron star spin evolution due to the r-mode instability is discussed in section 9. Our results are compared to those of previous investigators in section 8. We discuss the detectability of the gravitational wave signal in section 10, and give a summary and conclusions in section 11. Two appendices give detailed calculations of the turbulent cascade for stellar inertial modes, and the nonlinear force coefficients.

## 2  Saturation by Nonlinear Mode Coupling

We start by reviewing the equations of motion for the mode amplitudes, and then specialize to the weak and strong driving limits.

We will work in a reference frame co-rotating with the star. Expansion of the fluid displacements, relative to the co-rotating frame, in terms of the linear eigenmodes

$`\left[\begin{array}{c}𝝃(t)\\ \dot{𝝃}(t)\end{array}\right]={\displaystyle \underset{\alpha }{}}q_\alpha (t)\left[\begin{array}{c}𝝃_\alpha \\ i\omega _\alpha 𝝃_\alpha \end{array}\right]`$ (5)
leads to the following system of coupled oscillator equations for the dimensionless complex amplitudes $`q_\alpha (t)`$ (Schenk et al. 2002):

$`\dot{q}_\alpha +i\omega _\alpha q_\alpha `$ $`=`$ $`\pm \gamma _\alpha q_\alpha +{\displaystyle \frac{i}{2}}\omega _\alpha {\displaystyle \underset{\beta \gamma }{}}\kappa _{\alpha \beta \gamma }^{}q_\beta ^{}q_\gamma ^{}.`$ (6)
The left hand side of eq.6 represents an unforced oscillator of rotating frame frequency $`\omega _\alpha `$, while the terms on the right hand side are the driving ($`+`$) or damping ($``$) term and the nonlinear term, which is quadratic in $`q`$. In our notation, $`\kappa _{\alpha \beta \gamma }`$ is roughly the ratio of interaction energy to mode energy at unit amplitude. The rotating frame mode energy is $`E=2E_{\mathrm{unit}}|q|^2`$, where $`E_{\mathrm{unit}}`$ is a (arbitrary) unit of energy which we find convenient to set to $`E_{\mathrm{unit}}=0.5Mr_{}^2\mathrm{\Omega }^2`$ 222 The nonlinear interaction energy also scales as the rotation energy of the star. Here $`M`$ and $`r_{}`$ are the stellar mass and radius, and $`𝛀=\mathrm{\Omega }𝒆_z`$ is the angular velocity. The sum over modes $`_\beta `$ involves a sum over the mode $`(\omega _\beta ,𝝃_\beta )`$, with amplitude $`q_\beta `$, as well as its complex conjugate $`(\omega _\beta ,𝝃_\beta ^{})`$, with amplitude $`q_\beta ^{}`$ (see Schenk et al. 2002 for a detailed derivation, but note that the type of index denoted there by $`A`$ is denoted here by $`\alpha `$.)

### 2.1  the discrete limit

In the regime where the driving rate of the unstable mode is smaller than the damping rates of low order modes, the instability can be saturated by a transferal of energy to a small number of damped modes. We will begin by discussing the coupling between the “parent” r-mode, and two damped “daughter” modes. Although an idealization, this basic problem is soluble, and indicates which modes couple most strongly to the r-mode. We review the dynamics of such 3-mode networks, including the parametric instability, the equilibrium solution, and the linear and nonlinear stability of the equilibrium solution (Wersinger et al. 1980; Wu & Goldreich 2001; Dziembowski & Krolikowska 1985; Dimant 2000; Abarbanel et al. 1993).

In terms of the real amplitude and phase variables, defined by $`q=A\mathrm{exp}(i\phi )`$, the equations for a system of three modes are

$`\dot{A}_1`$ $`=`$ $`\gamma _1A_1\omega _1\kappa A_2A_3\mathrm{sin}\phi `$
$`\dot{A}_2`$ $`=`$ $`\gamma _2A_2\omega _2\kappa A_3A_1\mathrm{sin}\phi `$
$`\dot{A}_3`$ $`=`$ $`\gamma _3A_3\omega _3\kappa A_1A_2\mathrm{sin}\phi `$
$`\dot{\phi }`$ $`=`$ $`\delta \omega \kappa \mathrm{cos}\phi A_1A_2A_3\left[{\displaystyle \frac{\omega _1}{A_1^2}}+{\displaystyle \frac{\omega _2}{A_2^2}}+{\displaystyle \frac{\omega _3}{A_3^2}}\right]`$ (7)
where the index 1 refers to the parent and $`2`$ and $`3`$ refer to the two daughter modes. We have defined the frequency detuning $`\delta \omega =\omega _1+\omega _2+\omega _3`$, coupling coefficients $`\kappa _{123}=\kappa \mathrm{exp}(i\delta )`$, and the relative phase $`\phi =\delta +\phi _1+\phi _2+\phi _3`$. The qualitative features of the time evolution, such as equilibrium and stability, depend only on the three dimensionless parameters $`(\gamma _2+\gamma _3)/\gamma _1`$, $`\gamma _2/\gamma _3`$, and $`\delta \omega /\gamma _1`$.

The parametric instability (Landau & Lifshitz 1969; Dziembowski & Krolikowska 1985; Kumar & Goodman 1996; Wu & Goldreich 2001) is a mechanism by which the daughter mode amplitudes will grow exponentially if the parent mode amplitude exceeds a certain threshold 333 A simple example of parametric instability is a pendulum in which the length of the string is being varied periodically. See Landau & Lifshitz (1969). . The result is that energy can be quickly taken out of the parent mode and transferred to the daughter modes, providing a means to limit the parent mode’s amplitude. Furthermore, the growth rate of the daughters is larger than the growth rate of the parent so that the daughters can “catch up” to the parent even if they start from a lower amplitude.

The parametric instability can be derived in the approximation where the parent mode’s amplitude is much larger than the daughter modes’ amplitudes so that the influence of the daughters on the parent can be neglected. Performing a linear stability analysis of eq.7 (Landau & Lifshitz 1969), one finds that the daughters grow exponentially like $`\mathrm{exp}(t/\tau )`$ when the parent mode amplitude $`A_1`$ exceeds the critical value given by

$`A_1^2`$ $`=`$ $`{\displaystyle \frac{\overline{\gamma }_2\overline{\gamma }_3}{\kappa ^2\omega _2\omega _3}}\left[1+\left({\displaystyle \frac{\delta \omega }{\overline{\gamma }_2+\overline{\gamma }_3}}\right)^2\right],`$ (8)
where $`\overline{\gamma }_{2,3}=\gamma _{2,3}+1/\tau `$. In particular, the threshold at which the instability first starts to operate is given by eq. (8) at $`\tau =\mathrm{}`$,

$`A_1^2`$ $`=`$ $`{\displaystyle \frac{1}{\kappa ^2Q_2Q_3}}\left[1+\left({\displaystyle \frac{\delta \omega }{\gamma _2+\gamma _3}}\right)^2\right]`$ (9)
where $`Q_2=\omega _2/\gamma _2`$ and $`Q_3=\omega _3/\gamma _3`$ are the quality factors of the daughter modes. In the limit $`\gamma _2,\gamma _30`$ of negligible damping of the daughter modes, it is useful to consider in addition the threshold above which the daughter mode’s growth rate will exceed that of the parent mode. This is given by eq. (8) at $`\tau =\gamma _1^1`$, $`\gamma _2=\gamma _3=0`$:

$`A_1^2`$ $`=`$ $`{\displaystyle \frac{1}{\kappa ^2\omega _2\omega _3}}[\gamma _1^2+\delta \omega ^2/4].`$ (10)
We give an example showing the parametric instability in Fig. 1.

Once the parametric instability occurs, the daughter modes start to grow rapidly. We now discuss the conditions under which the subsequent evolution leads to a saturation of the parent mode in the three mode system.

Setting the time derivatives in eq.7 to zero, one finds the equilibrium solution for the parent

$`A_1^2`$ $`=`$ $`{\displaystyle \frac{1}{\kappa ^2Q_2Q_3}}\left[1+\left({\displaystyle \frac{\delta \omega }{\gamma _2+\gamma _3\gamma _1}}\right)^2\right],`$ (11)
and daughter mode energies $`A_{2,3}^2/A_1^2=Q_{2,3}/Q_1`$, where $`Q=\omega /\gamma `$ is the quality factor of the mode. Naively, one expects that energy transfer from the parent to the daughters occurs only if the daughter modes have a lower energy than the parent, implying a lower daughter mode quality factor. This expectation is verified by a stability analysis (Wu & Goldreich 2001; Dimant 2000) which shows that the equilibrium solution is stable only when $`\gamma _2+\gamma _3\mathrm{¿}\mathrm{}\gamma _1`$ is satisfied.

More precisely, there are three different regimes in the three dimensional space of parameters $`(\gamma _2+\gamma _3)/\gamma _1`$, $`\gamma _2/\gamma _3`$ and $`\delta \omega /\gamma _1`$. First, the equilibrium solution is linearly stable to small perturbations if two conditions are met (Wersinger et al. 1980): (i) the ratio of damping to driving is sufficiently large $`\gamma _2+\gamma _3\mathrm{¿}\mathrm{}\gamma _1`$, and (ii) the detuning is sufficiently large, $`\delta \omega \mathrm{¿}\mathrm{}(\gamma _2+\gamma _3)/2`$. Second, in the regime where $`\gamma _2+\gamma _3\mathrm{¿}\mathrm{}\gamma _1`$ but where the detuning is small, $`\delta \omega \mathrm{¡}\mathrm{}(\gamma _2+\gamma _3)/2`$, the the amplitudes and phase undergo limit cycles characterized by bounded, quasiperiodic orbits, as shown by Fig. 2 of Wu & Goldreich (2001). Those limit cycle solutions have time averaged parent mode amplitudes comparable to the equilibrium amplitude (11), so the equilibrium amplitude still characterizes the motion. Third, if the daughter mode damping is insufficient, $`\gamma _2+\gamma _3\mathrm{¡}\mathrm{}\gamma _1`$, all three amplitudes rise without bound and the solution is nonlinearly unstable (Dimant 2000). For our purposes, any solution which is nonlinearly stable can saturate the growth of the r-mode, so that the effective stability criterion is

$`\gamma _2+\gamma _3\mathrm{¿}\mathrm{}\gamma _1.`$ (12)
In the regime where the equilibrium solution is stable, it acts like an attractor, and the system tends to evolve into this equilibrium after the daughter mode amplitudes become comparable to that of the parent mode. The example shown in Fig. 1 exhibits this behavior, even though the system is started well away from equilibrium. Note that the equilibrium parent mode amplitude (11) is always approximately equal to the threshold amplitude (9), in the regime (12) where the energy transfer is stable.

The parametric instability can provide a means for saturating the r-mode amplitude. Suppose that a daughter pair exists which is parametrically unstable for a certain value $`A_1`$ of the parent mode amplitude, and that no other daughter pairs are unstable at that amplitude. Then, if the resonance is sharp, it is plausible that only the parent and two daughter modes are relevant, and if the condition (12) is satisfied so that the transfer of energy is stable, then driving of the r-mode by gravitational radiation reaction can be balanced by nonlinear energy transfer to the pair of daughter modes. Thus, the daughter mode pair for which the instability threshold (9) is lowest sets the saturation amplitude for the r-mode, if the stability constraint (12) is satisfied for that daughter mode pair. Daughter pairs with higher thresholds will not be excited because the parent’s amplitude cannot rise much above the lowest threshold (see Fig. 1).

The task of finding the saturation amplitude in the weak driving regime involves searching through all possible daughter mode pairs to minimize the parametric threshold (9). This amounts to maximizing $`\kappa `$ while minimizing the mismatch $`\delta \omega ^2+(\gamma _2+\gamma _3)^2`$ subject to the stability constraint. Once this “best” daughter mode pair has been found, the saturation amplitude is

$`A_1^2`$ $``$ $`{\displaystyle \frac{1}{\kappa ^2Q_2Q_3}}|_{\mathrm{best}\mathrm{pair}},`$ (13)
assuming that the strong resonance condition $`\delta \omega \mathrm{¡}\mathrm{}\gamma _2+\gamma _3`$ is satisfied. We can state the following rule of thumb: for coupling coefficients of order unity, the r-mode will saturate to an amplitude less than unity if the best daughter pair are high Q oscillators. Quality factors of low order global modes in neutron stars can easily be $`10^6`$ or larger.

Finding the saturation amplitude in the weak driving regime has now been reduced to the following physics problem. First determine the oscillation modes present in the star. Calculate their damping and driving rates, as well as the nonlinear coupling coefficients between daughter pairs and the r-mode. Once the magnitude and scalings of these quantities are known, reliable estimates of the parametric threshold can be made (see Sec. 6 below).

Finally, note that nonlinear coupling terms such as $`\dot{q}_2\kappa _{211}q_1^2`$ which couple the parent mode twice with a daughter mode have been ignored in eq.7. Since these terms scale as $`A_1^2`$, instead of $`A_1`$ as for the parametrically driven modes, they are smaller in the weakly nonlinear regime. In addition, the coupling coefficients drop off rapidly for this type of coupling as the wavenumber of mode $`2`$ is increased (see appendix B). Hence, only daughter modes with comparable wavenumber to the parent couple well. However, the resonance condition cannot be finely tuned for comparable wavenumber modes, since there are so few of them. As opposed to the couplings $`\kappa _{211}q_1^2`$, parametric type couplings have the double advantage of allowing coupling of the parent mode with daughter modes of arbitrarily large wavenumber, and the resonance condition becomes satisfied to a higher degree of precision for large daughter mode wavenumber.

### 2.2  the continuum limit

In the above “discrete” scenario, the saturation amplitude of the driven mode scales as $`AQ_d^1`$, where $`Q_d`$ is the quality factor of a damped mode. In the “continuum” picture that we now discuss, the saturation amplitude is independent of the linear damping rates, since the energy is transferred by nonlinear interactions. In this cascade picture, both the shape and normalization of the wave energy spectrum are set only by the detailed nonlinear interaction between waves, and the energy input to the system.

How does the cascade arise? Imagine starting with a system in the weak driving limit and adiabatically increasing the driving. When $`\gamma _{gr}`$ becomes greater than the damping rate $`\gamma _2+\gamma _3`$ of the daughter pair with the lowest threshold, the equilibrium solution for that pair is no longer stable and the energy of all three modes will begin to grow. When the energy has grown to the point that additional parametric thresholds are crossed, and if the energy transfer to these pairs is stable, the driven mode will again be saturated. As the driving is increased, this process will continue until many daughter modes are excited with large amplitude. Since linear damping is smaller than driving over a certain range of daughter mode lengthscales, an *inertial range* has formed where nonlinear forces are dominant. We now proceed to give a heuristic derivation of instability saturation in this continuum case, leaving the detailed derivation appendix A.

Since many modes are excited, we treat the quantum numbers for each mode as a continuum. Introducing the “occupation number” 444 To get a quantity with the units of action, multiply $`N_\alpha `$ by $`E_{\mathrm{unit}}`$. (quasi-particle number) for mode $`\alpha `$

$`N_\alpha =|q_\alpha |^2/|\omega _\alpha |0,`$ (14)
the mode energy becomes

$`E_\alpha `$ $`=`$ $`2E_{\mathrm{unit}}|q_\alpha |^2=2E_{\mathrm{unit}}|\omega _\alpha |N_\alpha .`$ (15)
Eq.6 describes both the fast variation of each individual mode, as well as the slow variation due to nonlinear interactions between modes. We may average over the fast oscillations using the random phase approximation (Zakharov et al. 1992; Kumar & Goldreich 1989; Wu 1998) if the phase randomization time set by the wave dispersion is shorter than the nonlinear interaction timescale. Since the dispersion time is comparable to the mode period for inertial waves, this is equivalent to the weak nonlinearity condition. The resultant kinetic equation for the wave amplitudes is (Zakharov et al. 1992; Kumar & Goldreich 1989)

$`\dot{N}_\alpha `$ $`=`$ $`I_\alpha +\mathrm{\Gamma }_\alpha N_\alpha `$ (16)
where $`I_\alpha `$ represents the rate of change of $`N_\alpha `$ due to nonlinear interactions, and $`\mathrm{\Gamma }_\alpha `$ is the rate of driving ($`>0`$) or damping ($`<0`$). The interaction term has the form

$`I_\alpha `$ $`=`$ $`s_\alpha \pi {\displaystyle \underset{\beta \gamma }{}}|\kappa _{\alpha \beta \gamma }|^2|\omega _\alpha \omega _\beta \omega _\gamma |\delta (\delta \omega _{\alpha \beta \gamma })`$ (17)
$`\times `$ $`\left(s_\alpha N_\beta N_\gamma +s_\beta N_\gamma N_\alpha +s_\gamma N_\alpha N_\beta \right)`$
where $`s_\alpha `$ is the sign of the frequency of mode $`\alpha `$.

To proceed further, we must introduce a few properties of the oscillation modes to be derived in section 3. Let $`n`$, $`k`$, and $`m`$ denote the perpendicular (to $`𝛀`$), parallel, and azimuthal number of nodes, respectively. Since inertial mode oscillation frequencies are proportional to the rotation frequency, we write the mode frequency as $`\omega =2\mathrm{\Omega }\mu `$, where $`\mu ^21`$ is the dimensionless frequency.

Approximate stationary solutions of eq.16 are found in two steps (see Zakharov et al. (1992) for detailed derivations). First, one ignores the driving and damping, so that the energy flux is conserved. In this case, the energy flux $``$ is defined by

$`\omega {\displaystyle 𝑑mI}`$ $``$ $`\mathbf{}_k`$ (18)
where $`\mathbf{}_k`$ is the gradient in momentum space and we have integrated over the $`m`$ quantum number. In appendix A, we show that stellar inertial waves support a flux of energy to large $`n`$ and small $`\mu `$. A schematic drawing of this cascade is given in fig.2. The occupation number 555The scaling of this expression for the occupation number can be simply derived from dimensional analysis together with the fact that 3-mode interactions dominate over 4-mode and higher order interactions. See, e.g., Sec. 3.3.1 of Zakharov et al. (1992).  for each mode is

$`N`$ $`=`$ $`\mathrm{\Omega }^1N_0n^4|\mu |^{1/2}n^{7/2}|k|^{1/2}`$ (19)
where the normalization constant is related to the fluxes $`^n`$ and $`^\mu `$ in the $`n`$ and $`\mu `$ directions, respectively, by

$`^n`$ $`=`$ $`8\alpha _nN_0^2n^1|\mu |^1\mathrm{\Omega }E_{\mathrm{unit}}`$
$`^\mu `$ $`=`$ $`8\alpha _\mu N_0^2n^2{\displaystyle \frac{\mu }{|\mu |}}\mathrm{\Omega }E_{\mathrm{unit}}.`$ (20)
The constants $`\alpha _n`$ and $`\alpha _\mu `$ are order unity and positive. A surface of constant energy in momentum space has $`\mu n^8`$, showing that the energy cascades to small frequencies quite rapidly with wavenumber, because of the strong dependence of the coupling coefficients on frequency (see appendix B).

The final step is to match the inertial range solution to the driving range. In other words, we need to find steady state solutions to the equation $`I+\gamma _{gr}N=0`$, where $`\gamma _{gr}`$ is the driving rate by gravitational radiation. We approximate $`N`$ in the driving region by extending the inertial range solution. The power input by the instability is given by 666 Even if a relatively narrow region $`\mathrm{\Delta }n^3`$ in phase space is being driven, Zakharov et al. (1992) find that one should use the whole volume $`n_d^3`$ instead of $`\mathrm{\Delta }n^3`$ since a peak develops in the driving region. Since the r-mode has relatively small wavenumber, the width of the driving region occupied by the r-mode may be considered relatively wide ($`\mathrm{\Delta }n/n1/3`$, $`\delta k/k1`$, $`\delta m/m1/2`$).

$`\text{input power}={\displaystyle 𝑑n𝑑k𝑑m\gamma _{gr}E}`$ (21)
$``$ $`{\displaystyle \frac{2}{3}}n_r^3\gamma _{gr}E_r={\displaystyle \frac{8}{3}}\gamma _{gr}E_{\mathrm{unit}}N_0n_r^1\mu _r^{1/2}`$
where the $`r`$ subscript denotes the driven r-mode and we have approximated the (dimensionless) volume in phase space being driven as $`n_r^3`$. The energy escaping the driving region is given by integrating up the flux through each boundary. Roughly, this is given by

output power $`=`$ $`{\displaystyle _{\mu _r(11/n_r)}^{\mu _r}}𝑑\mu n^n`$ (22)
$`+`$ $`{\displaystyle _{n_r}^{n_r+1}}𝑑nn|^\mu |`$
$``$ $`8n_r^1N_0^2(\alpha _n+\alpha _\mu )\mathrm{\Omega }E_{\mathrm{unit}}.`$
Equating the input to output power we can solve for the normalization constant

$`N_0`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{\alpha _n+\alpha _\mu }}{\displaystyle \frac{\gamma _{gr}}{\mathrm{\Omega }}}\mu _r^{1/2}.`$ (23)
Given the normalization, we find

$`{\displaystyle \frac{E_r}{E_{\mathrm{unit}}}}`$ $`=`$ $`{\displaystyle \frac{E(rmode)}{0.5Mr_{}^2\mathrm{\Omega }^2}}`$ (24)
$`=`$ $`{\displaystyle \frac{2n_r^4\mu _r}{3(\alpha _n+\alpha _\mu )}}{\displaystyle \frac{\gamma _{gr}}{\mathrm{\Omega }}}\alpha _e{\displaystyle \frac{\gamma _{gr}}{\mathrm{\Omega }}}`$
where we parametrize our inexact treatment of the matching condition with the parameter $`\alpha _e=2n_r^4\mu _r/3(\alpha _n+\alpha _\mu )`$. If we use the quantum numbers of the r-mode, $`n_r=3`$ and $`\mu _r=1/3`$, we find $`\alpha _e4\times 10^4`$. However, we choose to be very cautious about this factor since we are extrapolating a WKB treatment into the regime of low order modes 777 The detuning may become non-negligible for low order modes. In addition, the resonance width from $`\gamma _{\mathrm{gr}}`$ can become important for coupling directly to the unstable r-mode. . A more conservative estimate would be to set $`n_r=\mu _r=1`$, giving $`\alpha _e0.1`$. We will use the more conservative result for numerical work in the rest of this paper, but recall that it may overestimate the saturation amplitude by up to three orders of magnitude. Using the r-mode driving rate from eq. 62, the final result for the saturation energy is

$`{\displaystyle \frac{E(rmode)}{0.5Mr_{}^2\mathrm{\Omega }^2}}`$ $`=`$ $`10^6\left({\displaystyle \frac{\alpha _e}{0.1}}\right)\nu _{\mathrm{khz}}^5`$ (25)
where $`\nu _{\mathrm{khz}}`$ is the spin frequency in units of $`1000\mathrm{Hz}`$.

Why is the saturation amplitude so small? The factor $`\gamma _{gr}/\mathrm{\Omega }`$ is inevitable 888 If the energy transfer is local in frequency space, this scaling will also hold for interaction with other wave families, such as inertial-gravity modes. since the only quantity with the units of frequency in the nonlinear interaction rate is $`\mathrm{\Omega }`$. The numerical factor $`\alpha _e`$ depends on considerations such as the effective volume and area of the driving region, and the power in the driving region relative to the largest scale (energy bearing) waves.

Eq.25 is one of the central results of this paper. It applies when nonlinear energy transfer is faster than linear damping. If nonlinear energy transfer becomes slower than linear damping, the discrete limit of section 2.1 is recovered. Note that the saturation amplitude decreases very rapidly with stellar spin frequency.

## 3 Oscillation modes in rapidly rotating neutron stars

In this section we discuss the oscillation modes present in rapidly rotating neutron stars. We argue that at the rapid rotation rates of interest for the r-mode instability, the buoyancy and elastic restoring forces can be ignored in comparison with the Coriolis force. The resulting modes which are restored by the Coriolis force are called inertial modes, of which the r-modes are a subset.

### 3.1  motivation for inertial waves

Within a minute after their birth in a supernova, neutron star cores have become transparent to neutrinos and cooled down sufficiently to form a degenerate gas of neutrons, with a small admixture of electrons and protons determined by beta equilibrium. As shown clearly by Reisenegger & Goldreich (1992), the varying electron fraction $`y_e6.0\times 10^3(\rho /2.8\times 10^{14}\mathrm{g}\mathrm{cm}^3)`$ in the star causes a stable stratification and resulting buoyancy force: Since the neutron pressure $`p_n[\rho /(1+y_e)]^{5/3}`$, displacing a fluid element upward on timescales slower than the sound crossing time and faster than the timescale of the beta reactions results in the fluid element being heavier than its surroundings, since it came from a region of larger $`y_e`$. An oscillatory motion results, with maximum frequency of order the Brunt-Vaisala frequency (Reisenegger & Goldreich 1992) $`N_{bv}(0.5y_eg/H_p)^{1/2}500\mathrm{sec}^12\pi \times 100\mathrm{Hz}`$, where $`g`$ is the local gravity and $`H_p`$ is the local pressure scale height.

The buoyancy force on a fluid element is just $`F_{\mathrm{buoy}}=g\delta \rho \rho N_{bv}^2\xi ^r`$, where $`\xi ^r`$ is the radial component of the Lagrangian fluid displacement. The Coriolis force is given by $`F_{\mathrm{cor}}2\rho \mathrm{\Omega }\omega \xi `$, so that the ratio of these two forces for $`\omega \mathrm{\Omega }`$ is roughly

$`{\displaystyle \frac{F_{\mathrm{buoy}}}{F_{\mathrm{cor}}}}`$ $``$ $`\left({\displaystyle \frac{N_{bv}}{\mathrm{\Omega }}}\right)^2\left({\displaystyle \frac{100\mathrm{Hz}}{\nu _{\mathrm{spin}}}}\right)^2.`$ (26)
In a detailed study of the solutions to the fluid perturbation equations for rotating stars including buoyancy, Yoshida & Lee (2000) showed that in the limit of $`F_{\mathrm{buoy}}/F_{\mathrm{cor}}1`$ the solutions are approximated very well by the r- and g-modes. This limit of large buoyancy force was examined by Morsink (2002) who showed that the nonlinear couplings between r-modes are too small to cause saturation to occur. In the limit of $`F_{\mathrm{buoy}}/F_{\mathrm{cor}}1`$ Yoshida & Lee (2000) have shown that the solutions of the perturbation equations are well approximated by the inertial modes. As long as we restrict our calculations to stars spinning at a frequency greater than 100 Hz, the inertial modes with zero buoyancy are very good approximation. As we are interested in the possibility of mode saturation at spin frequencies at least as large as 300 Hz, the inertial modes are the most relevant modes and it is possible for us to ignore the buoyancy force as a first approximation. This enables us to find simple solutions for the modes if we further approximate the shape of the star as spherical, a valid approximation for rotation rates well below the breakup rate. However, we expect that the qualitative results found here will hold even in the case when buoyancy is included. The reason is that the approximations made still provide a dense spectrum of modes that may be arbitrarily resonant with the r-mode in the continuum limit. Although the numerical value of the coupling coefficients and damping rates may change because we don’t have exactly the correct shape of the eigenfunctions, we are confident that the essential qualitative features present in our simple example will carry over.

Levin & Ushomirsky (2001) have shown that the elastic restoring force in the neutron star crust becomes small compared to the Coriolis force above a rotation rate of $`50\mathrm{Hz}`$. The net result is that core modes can penetrate into the crust, with only a small discontinuity at the crust-core boundary because of the impedance mismatch. We will ignore crustal elasticity for the remainder of this paper.

We have not included superfluid effects in our calculations. The principal new effect would be dissipation due to mutual friction (the modes themselves are not expected to be changed very much; see, e.g., Lindblom & Mendell (2000)). However, we note that our estimate of the saturation amplitude does not depend on the dissipation rate if an inertial cascade forms (although the outer scale of the inertial range might be affected). Thus, if saturation involves a cascade of energy to numerous inertial modes, we still expect our estimates to hold. Our estimates would change if dissipation via mutual friction is strong enough that only a few modes are excited parametrically. We postpone a thorough examination of this case, which would depend on uncertain mutual friction coefficients, for another paper. However, either way, the saturation amplitude will still be very small.

In the next subsection, we discuss inertial mode eigenfunctions in weakly stratified stars.

### 3.2  stellar inertial modes

We solve the Euler and continuity equations for adiabatic perturbations of a rotating star. The background star is taken to be spherically symmetric with uniform rotation rate $`𝛀=\mathrm{\Omega }𝒆_z`$ and negligible stable stratification. Perturbation modes of the form $`\mathrm{exp}(im\varphi i\omega t)`$ are found using the Cowling approximation.

The Euler, continuity, and state equations are

$`\rho \ddot{𝝃}+2\rho 𝛀\times \dot{𝝃}`$ $`=`$ $`\mathbf{}\delta p+𝒈\delta \rho `$ (27)
$`\delta \rho +\mathbf{}\left(\rho 𝝃\right)`$ $`=`$ $`0`$ (28)
$`\delta p`$ $`=`$ $`c^2\delta \rho ,`$ (29)
where we have ignored terms involving the Brunt-Vaisala frequency

$`N_{bv}`$ $`=`$ $`\left(g{\displaystyle \frac{d\mathrm{ln}(p^{1/\mathrm{\Gamma }_1}/\rho )}{dr}}\right)^{1/2},`$ (30)
a valid assumption for $`\omega N_{bv}`$ and $`\mathrm{\Omega }N_{bv}`$. Here we have defined the adiabatic index $`\mathrm{\Gamma }_1`$. In this limit, the adiabatic sound speed $`c`$ and density scale height $`H`$ are related by $`c(gH)^{1/2}`$. Substituting the assumed dependence on $`\varphi `$ and $`t`$, and eliminating $`\delta \rho `$, we find that eqns. 2729 become

$`𝝃+iq𝒆_z\times 𝝃`$ $`=`$ $`\mathbf{}\psi `$ (31)
$`\mathbf{}𝝃+{\displaystyle \frac{\omega ^2}{c^2}}\psi `$ $`=`$ $`{\displaystyle \frac{\xi ^r}{H}}.`$ (32)
Here we have replaced the Eulerian pressure perturbation by the quantity $`\psi `$ defined by $`\delta p=\rho \omega ^2\psi ,`$ and defined the dimensionless inverse frequency $`q=2\mathrm{\Omega }/\omega `$. We will also heavily use the dimensionless frequency $`\mu =1/q=\omega /2\mathrm{\Omega }`$. We drop the term $`\omega ^2\psi /c^2\mathrm{\Omega }^2/(GM/r_{}^3)`$ since we are working to leading order in $`\mathrm{\Omega }`$; this is consistent with our assumption that the background star is spherical and suffices to compute the mode functions to leading order in $`\mathrm{\Omega }`$.

The Euler equation 31 can be solved 999 The determinant of this transformation is singular only if $`\omega ^2=4\mathrm{\Omega }^2`$. for $`𝝃`$ in terms of $`\psi `$:

$`𝝃`$ $`=`$ $`(1q^2)^1\left[\mathbf{}\psi q^2𝒆_z(𝒆_z\mathbf{}\psi )iq𝒆_z\times \mathbf{}\psi \right].`$ (33)
Substituting eq. 33 into the continuity equation 32 gives the wave equation

$`\mathbf{}^2\psi q^2{\displaystyle \frac{^2\psi }{z^2}}`$ $`=`$ $`H^1\left({\displaystyle \frac{\psi }{r}}q^2\mathrm{cos}\theta {\displaystyle \frac{\psi }{z}}{\displaystyle \frac{mq}{r}}\psi \right).`$ (34)
The boundary condition near the surface is that the Lagrangian change in the pressure vanish, $`\mathrm{\Delta }p=\delta p+𝝃\mathbf{}p=0`$, so that $`\delta p\rho g\xi ^r`$. This is just the statement that the surface layer is hydrostatic, a consequence of the vanishing sound crossing time across a scale height for small depth.

Equation 34 does not appear to be solvable by separation of variables101010By separation of variables, we mean that (1) the differential equation is separable, and (2) the boundary conditions are applied on a surface where one of the coordinates is constant.. This motivates us to examine approximate solutions valid for short wavelengths. Our solution generalizes the exact solution of Bryan (1889) for the constant density star; in fact our solution is just Bryan’s solution divided by $`\sqrt{\rho }`$.

Defining

$`\psi (𝒙)`$ $`=`$ $`\psi _0\left({\displaystyle \frac{\rho _0}{\rho }}\right)^{1/2}f(𝒙)`$ (35)
where $`\psi _0`$ is a normalization constant and $`\rho _0`$ is the central density, the differential equation for $`f`$ is

$`\mathbf{}^2fq^2{\displaystyle \frac{^2f}{z^2}}+𝒦^2f`$ $`=`$ $`0`$ (36)
with

$`𝒦^2`$ $`=`$ $`(1q^2\mathrm{cos}^2\theta )\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{dH^1}{dr}}{\displaystyle \frac{1}{4H^2}}\right)`$ (37)
$`+`$ $`{\displaystyle \frac{mq}{rH}}+{\displaystyle \frac{1}{2rH}}\left(2q^2\mathrm{sin}^2\theta \right).`$
The first two terms in eq. 36 are just the usual differential equation for inertial modes, as derived by Bryan. The compression term $`\xi ^r/H`$ in the continuity equation 32 is imaginary in the WKB limit, and leads to the WKB envelope $`\rho ^{1/2}`$. The definition in eq.35 accounts for this envelope, so that the correction terms in eq.37 are now real. The $`𝒦^2f`$ term is most important near the surface, where it scales as $`(|𝒌|H)^2`$ relative to the other terms in eq.36 ($`|𝒌|`$ is the local WKB wavenumber). This term describes a slow variation of the wavenumber with position, and is negligible for short wavelength modes. Henceforth, we set $`𝒦0`$ in the interior of the star.

The short wavelength approximation breaks down when the vertical wavenumber $`k_r`$ 111111 One must use care evaluating $`k_r`$ for inertial modes near the surface of the star, since it can vary strongly with the angle $`\theta `$. Qualitatively, this strong variation occurs because one is imposing a spherical boundary condition on waves with inherent cylindrical symmetry. is comparable to the scale height, at about $`2k_rH1`$, as can be verified directly from eq. 36. When the eigenfunction is constant over a scale height one may picture that the wave attempts to lift an entire scale height of material, which causes reflection. We extend our interior solution to the surface by replacing $`\rho `$ in eq. 35 with a “cutoff” value $`\rho _{cut}`$ which becomes constant about one wavelength from the surface. 121212 If the density profile near the surface is a power law with depth, one can separate variables in the bi-spheroidal coordinates introduced below. These more rigorous solutions close to the surface agree with the cutoff behavior described here for the density. Although one could, in principle, match the interior WKB solution to the surface solution, the cutoff for the density gives an adequate approximation for the problem at hand.

Bryan (1889) found a solution to eq.36 for $`𝒦=0`$ in terms of an ingenious bi-spheroidal coordinate system that depends on the frequency $`\mu `$. Lindblom & Ipser (1999) have given a careful discussion of this coordinate system, paying particular attention to the behavior of the coordinates at the surface of the star. Define the bi-spheroidal coordinates $`|\mu |x_11`$ and $`|\mu |x_2|\mu |`$ by

$`x`$ $`=`$ $`r_{}\left[{\displaystyle \frac{(1x_1^2)(1x_2^2)}{1\mu ^2}}\right]^{1/2}\mathrm{cos}\varphi `$
$`y`$ $`=`$ $`r_{}\left[{\displaystyle \frac{(1x_1^2)(1x_2^2)}{1\mu ^2}}\right]^{1/2}\mathrm{sin}\varphi `$
$`z`$ $`=`$ $`r_{}{\displaystyle \frac{x_1x_2}{|\mu |}}.`$ (38)
In fig. 3 we plot the surfaces of constant bi-spheroidal coordinate in the $`xz`$ plane. In either the $`\mu 1`$ or $`\mu 1`$ limits, the surfaces of constant $`(x_1,x_2)`$ are nearly in the $`z`$ and $`R`$ directions over most of the star, as one would expect for a local plane wave propagating in the $`z`$ or $`R`$ direction. (Here $`R`$ is the cylindrical radius.) However, the coordinate lines near the point $`rr_{}`$ and $`\mathrm{cos}\theta =\pm |\mu |`$ on the surface of the star vary quite rapidly with respect to $`r`$ and $`\theta `$. The result is that the WKB wavenumber becomes quite large near these singular points. (We give a detailed mathematical discussion in appendix B.) As one can see from fig. 3, the coordinate lines come closer to the surface near the singular points, implying the upper turning point is much closer to the surface near the equator (for small $`\mu `$) than the poles. As a result, the wave amplitudes will be much larger near the equator, as we will now show.

Our approximate solution for the interior of the star is to ignore terms of order $`(|𝒌|H)^2`$, so that the differential equation becomes

$`\mathbf{}^2fq^2{\displaystyle \frac{^2f}{z^2}}`$ $`=`$ $`0.`$ (39)
Changing to bi-spheroidal coordinates in eq. 39 gives separable differential equations (see Bryan (1889) and Lindblom & Ipser (1999) for details) . Define the solution $`f(x_1,x_2)=f_1(x_1)f_2(x_2)`$ where both $`f_1`$ and $`f_2`$ satisfy

$`{\displaystyle \frac{}{x}}\left[(1x^2){\displaystyle \frac{f}{x}}\right]+\left(\kappa ^2{\displaystyle \frac{m^2}{1x^2}}\right)f`$ $`=`$ $`0`$ (40)
for separation constant $`\kappa ^2`$. This equation has the Legendre function solutions $`\kappa ^2=n(n+1)`$, $`f_1(x_1)=P_{nm}(x_1)`$ and $`f_2(x_2)=P_{nm}(x_2)`$.

The resulting solution for $`\psi `$ is then

$`\psi (𝒙,t)`$ $`=`$ $`\psi _0\left({\displaystyle \frac{\rho _0}{\rho }}\right)^{1/2}P_{nm}(x_1)P_{nm}(x_2)e^{im\varphi i\omega t}.`$ (41)
Note the important fact that this solution is valid for an arbitrary density profile $`\rho `$, so long as one is safely in the short wavelength limit. This is true even when $`\rho `$ is not a separable function of $`x_1`$ and $`x_2`$, as is generally the case in the interior since $`\rho (r)=\rho \left(\sqrt{1(x_1^2\mu ^2)(\mu ^2x_2^2)/\mu ^2(1\mu ^2)}\right)`$.

The r-modes do not have short wavelengths and hence cannot be described by the above WKB approximation. However, in the leading order approximation of a spherical background star with no buoyancy, the r-mode solutions are given by (Bryan 1889)

$`\psi (𝒙,t)`$ $`=`$ $`\psi _0P_{|m|+1,m}(x_1)P_{|m|+1,m}(x_2)\mathrm{exp}(im\varphi i\omega t)`$ (42)
$``$ $`zR^{|m|}\mathrm{exp}(im\varphi i\omega t)`$
and have frequencies $`\mu =\mathrm{sign}(m)/(|m|+1)`$.

We now derive the WKB limit for the solution in eq. 41. Writing $`x_{1,2}=\mathrm{cos}\theta _{1,2}`$ and substituting $`f\mathrm{sin}^{|m|}\theta \mathrm{exp}(i^\theta 𝑑\theta k)`$ in eq. 40, we find the following standing wave solutions:

$`f(\theta )`$ $``$ $`{\displaystyle \frac{\pi ^1}{\mathrm{sin}^{1/2}\theta }}\mathrm{cos}(p\theta +\alpha )`$ (43)
where the wavenumber is given by

$`p`$ $`=`$ $`\left(n(n+1)|m|[|m|+1]\right)^{1/2}n`$ (44)
for $`n|m|`$ (the WKB limit) and

$`\alpha `$ $`=`$ $`p\pi /2\text{ for even parity}`$
$`\alpha `$ $`=`$ $`(p+1)\pi /2\text{ for odd parity modes}.`$ (45)
We have chosen to normalize the Legendre polynomials to unity over $`4\pi `$ steradians. Note that the nodes are spaced evenly in $`\theta _{1,2}`$. This WKB approximation to the Legendre equation fails within about one wavelength of $`\theta =0,\pi `$. The $`\mathrm{sin}\theta _{1,2}`$ factor causes an increase in amplitude toward $`\mathrm{sin}\theta _{1,2}=0`$. The collected result is then

$`\psi (𝒙)`$ $`=`$ $`{\displaystyle \frac{\psi _0}{\pi ^2}}{\displaystyle \frac{\mathrm{cos}(p\theta _1+\alpha )\mathrm{cos}(p\theta _2+\alpha )}{\left(\frac{\rho }{\rho _0}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\right)^{1/2}}}e^{im\varphi }.`$ (46)
The factor in the denominator $`\rho \mathrm{sin}\theta _1\mathrm{sin}\theta _2\rho R`$ is just the mass element, and enforces roughly equal kinetic energy in between each pair of nodes.

An approximate dispersion relation is easily derived using the eigenfunctions of eq.46. The boundary condition is that the compression term $`\xi ^r/H`$ in eq. 32 must remain finite as $`H0`$, implying $`\xi ^r0`$ in the low frequency approximation. At either surface patch $`x_1=|\mu |`$ or $`|x_2|=|\mu |`$, this condition implies

$`(1x^2){\displaystyle \frac{dP_{nm}}{dx}}+mP_{nm}`$ $`=`$ $`0`$ (47)
at $`x=\pm |\mu |`$. Eq.47 is equivalent to the one given by Bryan (1889), as noted by Lindblom & Ipser (1999). The dependence of the frequency and wavenumber on the background stellar model, as discussed by Lockitch & Friedman (1999), is small in the WKB limit. Substituting the WKB expressions gives

$`\mathrm{sin}\theta \mathrm{tan}(p\theta +\alpha )`$ $`=`$ $`{\displaystyle \frac{m}{p}}.`$ (48)
In the limit $`p|m|`$, the solutions are found by inspection to be $`p\theta +\alpha =k\pi `$, for the mode index $`k`$. Including the finite $`m`$ term to first order gives the solution

$`\mu _{nkm}`$ $`=`$ $`\omega _{nkm}/2\mathrm{\Omega }`$ (49)
$`=`$ $`\mathrm{sin}(k\pi /p+m/p^2\pi \delta /2p)`$
$``$ $`k\pi /n`$ (50)
where $`\delta =0`$ for even parity modes and $`\delta =1`$ for odd parity modes. The $`m`$ term can be dropped except for the very low frequency, even parity mode with frequency $`\mu m/p^2`$. All other modes have $`\mu p^1`$. We find the approximate formula in eq.49 to agree quite well with the exact solutions of eq.47 even for $`p`$ as small as 5. Lastly, we note that Lockitch & Friedman (1999) have checked the eigenmodes found using bi-spheroidal coordinates with those from a numerical code in spherical coordinates, finding agreement.

We choose to normalize the eigenfunctions so that at unit amplitude ($`A=1`$) all modes have the same energy, which we call $`2E_{\mathrm{unit}}`$. We can analytically compute the mode energy in the WKB limit where the eigenfunctions are rapidly oscillating ($`nm`$), with the result

$`E`$ $`=`$ $`\omega ^2{\displaystyle d^3x\rho |𝝃|^2}`$ (51)
$``$ $`4\pi ^2{\displaystyle \frac{p^2\mu ^4}{(1\mu ^2)^{3/2}}}\rho _0\mathrm{\Omega }^2r_{}\psi _0^2.`$
This formula agrees well with numerical integrations. Our normalization convention is that at unit amplitude all modes have the same energy $`2E_{\mathrm{unit}}`$. We then find the value of the normalization constant

$`\psi _0^2`$ $`=`$ $`{\displaystyle \frac{E_{\mathrm{unit}}}{2\pi ^2\frac{p^2\mu ^4}{(1\mu ^2)^{3/2}}\rho _0\mathrm{\Omega }^2r_{}}}p^2\mu ^4.`$ (52)
Modes with rapid spatial variation ($`p1`$) or larger frequency $`\mu `$ have smaller normalization in order for the energy to be the same. As $`\mu 1`$, the wave amplitude goes to zero since inertial modes do not exist outside this range.

Before moving on to discuss the nonlinear force coefficients, we discuss the normalization integral in eq.51. One can easily find the mode energy to leading order in $`\mu `$ by setting $`\mu =0`$ in eq.51. In this limit, the bi-spheroidal coordinates become $`x_1(1R^2)^{1/2}`$ and $`x_20`$. In this limit, the integrand is constant in $`z`$, and varies as $`(1R^2)^{1/2}`$ with $`R`$, which is large near the surface. The kinetic energy then converges as $`(1R^2)^{1/2}`$ from the surface.

## 4  coupling coefficients

The lowest order nonlinear interaction couples three inertial waves, implying quadratic nonlinear terms as in eq.6. The expressions for the nonlinear force coefficients can be derived either from an action principle (Newcomb 1962; Kumar & Goldreich 1989; Kumar & Goodman 1996) or directly from the equation of motion (Schenk et al. 2002). Note that Schenk et.al. have stressed that the form of the coupling coefficient is the same for rotating systems as for nonrotating systems; only the explicit expressions for the eigenfunctions and background stellar model need be modified. Since we are using daughter modes with wavelengths much smaller than a stellar radius, we keep only the largest 131313 Inertial waves in an infinite homogeneous, incompressible medium have a nonlinear coupling as given here. Including terms arising from compressibility or variation of the background stellar quantities then gives terms which are small in the limit $`|𝒌|H1`$. We ignore these terms here for simplicity, although the r-mode is formally a large lengthscale mode.  term in the coupling coefficient in an expansion of $`(|𝒌|H)^1`$. For modes $`𝝃_1`$, $`𝝃_2`$, $`𝝃_3`$, the dimensionless coupling coefficient 141414 see Schenk et.al. for a derivation of eq.6 and the explicit form of the dimensionless coupling coefficient $`\kappa `$ is

$`\kappa _{123}`$ $``$ $`{\displaystyle \frac{1}{2E_{\mathrm{unit}}}}{\displaystyle }d^3x(\xi _1^i\xi _2^j\delta p_{3;ij}`$ (53)
$`+`$ $`\xi _2^i\xi _3^j\delta p_{1;ij}+\xi _3^i\xi _1^j\delta p_{2;ij}).`$
Since $`\delta p\mathrm{\Omega }^2`$, we find that the interaction energy, $`\kappa E_{\mathrm{unit}}`$, scales as the rotational energy of the star. A natural unit of energy is then $`E_{\mathrm{unit}}=0.5Mr_{}^2\mathrm{\Omega }^2`$. In these units, $`A^2`$ is the mode energy in units of $`2E_{\mathrm{unit}}=Mr_{}^2\mathrm{\Omega }^2`$.

In section 4.1 we discuss conservation rules for the nonlinear coupling coefficients. Effectively, these rules pick out the largest possible coupling coefficients. The scalings for $`\kappa `$ are discussed in section 4.2. We confirm a result found in previous studies (Wu & Goldreich 2001) that for waves which satisfy the conservation rules, the coupling coefficients do not become smaller as the daughter mode wavenumber is increased; even though each individual eigenfunction is highly oscillatory, the product is relatively constant. Numerical results are presented in section 4.3 and a detailed analytic calculation is given in appendix B.

### 4.1  energy and momentum conservation

Consider a parent mode with quantum numbers $`(n_1,k_1,m_1)`$ and frequency $`\mu _1`$. We are free to choose daughter mode quantum numbers $`(n_2,k_2,m_2)`$ and $`(n_3,k_3,m_3)`$ in order to find the largest coupling coefficient (see e.g. Wu & Goldreich (2001)). The integrand is highly oscillatory unless the phases of the waves match at each point in the star. If we expand the standing wave solution in eq.46 in terms of travelling waves, a non-oscillatory integrand implies momentum conservation for the three travelling waves. In addition to conservation of the $`m`$ quantum number, due to axisymmetry of the background star, we also have momentum conservation along the $`\theta _1`$ and $`\theta _2`$ directions. For small $`\mu `$, the total number of nodes along $`\theta _1`$ and $`\theta _2`$ simplifies to $`N_1pn`$ and $`N_2|\mu |p/\pi k`$. The approximate conservation laws which lead to large $`\kappa `$ can then be written

$`m_1+m_2+m_3=0\text{ angular momentum}`$
$`|n_2n_3|\mathrm{¡}\mathrm{}n_1\text{ momentum along }\theta _1`$
$`||k_2||k_3||\mathrm{¡}\mathrm{}|k_1|\text{ momentum along }\theta _2\text{ }.`$ (54)
For small frequency, the $`\theta _1`$ and $`\theta _2`$ directions lie nearly along the $`R`$ and $`z`$ directions, so that the second and third momentum conservation rules correspond to conservation of momentum along $`R`$ and $`z`$. In the limit that the daughter modes have much smaller wavelengths than the parent mode, which will turn out to be the important limit, we find the simple result $`n_2n_3`$ and $`|k_2||k_3|`$; momentum conservation implies that the daughter modes have momenta of equal magnitude and oppositely directed.

So far, we have used momentum conservation to determine three of the six daughter mode quantum numbers. In order for energy to be efficiently transferred between modes, the interaction must be as nearly resonant as possible, meaning that the detuning is small:

$`\delta \omega =\omega _1+\omega _2+\omega _30.`$ (55)
There are two simple limits of interest. For short wavelength daughter modes with $`n_2n_3`$ and $`k_2k_3`$, one has $`\mu _2\mu _3\mu _1/2`$; the parent mode interacts with nearly identical daughter modes of half the frequency of the parent. The second solution is where $`n_1n_2`$, and $`n_3n_1,n_2`$. In this case we find $`n_2n_1+n_3`$, $`k_2k_1`$, and $`k_3k_1(n_3/n_1)^2`$. The frequencies are then $`\mu _2\mu _1`$ and $`\mu _3\mu _1(n_3/n_1)\mu _1`$.

### 4.2  analytic estimates

Here we give a back of the envelope estimate for the coupling coefficient, leaving the more detailed calculation for appendix B. We will only consider the important limit of short wavelength, nearly identical daughter modes with $`\mu _2\mu _3\mu _1/2`$. We shall ignore factors of order unity for the present, concentrating only on the scalings. As $`\kappa `$ is dimensionless, we set $`r_{}=1`$ in this section for simplicity.

Incompressibility of the waves implies

$`\kappa {\displaystyle \frac{1}{E_{\mathrm{unit}}}}{\displaystyle d^3x\rho \omega _1^2k_{1z}^2\psi _1\xi _{2z}^2}.`$ (56)
For the polytrope of index $`1`$ we find $`\rho r_{}r\widehat{z}`$ near the surface, where $`\widehat{z}`$ is the distance from the surface. Since the WKB envelope of the waves rises steeply toward the surface, and the factor of $`\rho `$ cancels the $`\rho ^1`$ from $`𝝃_{2z}^2`$, we find that the dominant contribution comes above the turning point for the parent mode, where

$`\psi _1`$ $``$ $`{\displaystyle \frac{1}{p_1\mu _1^2\widehat{z}_1^{1/2}}}.`$ (57)
Here $`\widehat{z}_1`$ is the turning point depth of the parent mode. The daughter mode eigenfunction is strongly peaked in the $`\theta `$ direction due to the wavenumber

$`k_{2z}{\displaystyle \frac{p_2\mu _2x_2}{\left[(\mathrm{cos}^2\theta \mu _2^2)^2+8\mu _2^2\widehat{z}\right]^{1/2}}}`$ (58)
where $`\theta `$ is the polar angle in spherical coordinates. The displacement for the daughter mode is then

$`\xi _{2z}`$ $``$ $`\widehat{z}^{1/2}\left[(\mathrm{cos}^2\theta \mu _2^2)^2+8\mu _2^2\widehat{z}\right]^{1/2}.`$ (59)
For $`\mathrm{cos}\theta |\mu _2|`$, $`k_{1z}p_1`$ since it is well away from the singularity for mode 1 at $`\mathrm{cos}\theta =|\mu _1|`$. Plugging into eq.56 we find

$`\kappa p_1\widehat{z}_1^{1/2}{\displaystyle _{\widehat{z}_2}^{\widehat{z}_1}}𝑑\widehat{z}{\displaystyle _1^1}{\displaystyle \frac{d(\mathrm{cos}\theta )}{(\mathrm{cos}^2\theta \mu _2^2)^2+8\mu _2^2\widehat{z}}}.`$ (60)
The integrand has a width $`d(\mathrm{cos}\theta )\widehat{z}^{1/2}`$ and a height $`(\mu _2^2\widehat{z})^1`$, giving an area $`(\mu _2^2\widehat{z}^{1/2})^1`$. Using $`|\mu _2||\mu _1|`$ and $`\widehat{z}_1\widehat{z}_2`$, the final result is then

$`\kappa `$ $``$ $`{\displaystyle \frac{p_1}{\mu _1^2\widehat{z}_1^{1/2}}}{\displaystyle _{\widehat{z}_2}^{\widehat{z}_1}}{\displaystyle \frac{d\widehat{z}}{\widehat{z}^{1/2}}}p_1\mu _1^2.`$ (61)
The detailed calculation in appendix B confirms that the coefficient is about unity.

We now comment on the scalings for the maximum coupling coefficient in eq.61. The maximum coupling coefficient is found to be independent of the daughter mode quantum numbers. The reason, elucidated by Wu & Goldreich (2001), is that one is integrating over the daughter mode kinetic energy $`\rho \mu _2^2\xi _2^2`$. This quantity is normalized to $`2E_{\mathrm{unit}}`$ when integrated over the whole star, and is roughly $`E_{\mathrm{unit}}\times \widehat{z}_1^{1/2}`$ when integrated over $`0\widehat{z}\widehat{z}_1`$. Next, the factor $`p_1n_1`$ implies shorter wavelength parent modes interact more strongly. This factor would appear for coupling of local waves in a box. However, the factor $`\mu _1^2`$ would not appear for local waves in a box; it arises from the large peak in the integrand near the surface.

One might wonder whether or not the approximate conservation laws for colliding WKB waves will hold since one is integrating over a small region of the star. The dominant contribution to the integrand comes from a region of size $`\widehat{z}\widehat{z}_1|\mu _1|/n_1`$, and the angular size is $`d(\mathrm{cos}\theta )\widehat{z}^{1/2}(|\mu _1|/n_1)^{1/2}`$. The daughter modes have wavelengths $`n_2`$ or $`\mu _2n_2`$, depending on direction, so there are still sufficient oscillations in the important region of the star for large enough $`n_2`$.

### 4.3  numerical calculation

We compute the integral in eq.53 numerically as follows. Choose a point in the star at which to evaluate the integrand. Evaluate $`\psi `$ and $`\delta p`$ on the vertices of a Cartesian cube about this point. The derivatives in eq.33 and 53 can then be taken along Cartesian basis vectors 151515 We evaluate vector quantities along Cartesian basis vectors to avoid “curvature terms” (Wu & Goldreich 2001) arising from differentiating curvilinear basis vectors. Wu and Goldreich found these terms are quite large, and cancel out in the end, so that significant cancellation error can occur. We avoid such cancellation error by using Cartesian basis vectors. , and then appropriate sums over indices taken. The resulting scalar integrand is independent of the coordinate $`\varphi `$ since $`m_1+m_2+m_3=0`$, so that only a two-dimensional integral over $`r`$ and $`\theta `$ remains. We perform this integration with second order accuracy, and increase the number of grid points until the integral converges.

In fig.4, we show the numerical integrations for the coupling coefficient as a function of $`n_1`$, but fixed $`k_1`$ and $`m_1`$. We also fix $`(k_2,m_2)`$ and $`(k_3,m_3)`$ but allow $`n_2`$ and $`n_3`$ to vary. For a given $`n_1`$ we see there is a variation in $`\kappa `$ due to the degree of momentum conservation. However, the upper envelope set by the maximum coupling coefficient agrees to within $`10\%`$ of our analytic formula.

## 5  damping and driving rates

We review the driving rate by gravitational radiation, and derive simple analytic estimates for the damping rates of inertial modes.

### 5.1  driving rate

Gravitational radiation reaction is a driving force if the phase velocity in the azimuthal direction is positive in the inertial frame and negative in the rotating frame; otherwise it damps the mode (Friedman & Morsink 1998). The driving rate falls off extremely rapidly with wavenumber, so that only the very lowest modes have an appreciable driving rate compared to damping. Lockitch & Friedman (1999) have computed these driving rates for the inertial modes of a polytrope of index 1, and identified several low order driven modes. However, since the most unstable mode by far is the $`(n,k,m)=(3,1,2)`$ r-mode, we can ignore all the others to a good approximation.

The driving rate of the $`(n,k,m)=(3,1,2)`$ r-mode for a polytrope of index 1 with $`M=1.4M_{}`$ and $`r_{}=12km`$ is (Lockitch & Friedman 1999)

$`\gamma _{gr}`$ $`=`$ $`0.05\mathrm{sec}^1\nu _{\mathrm{khz}}^6.`$ (62)
### 5.2  bulk viscosity damping

We now compute the damping rate of inertial modes by bulk viscosity damping due to the modified URCA processes. We take the coefficient of bulk viscosity from Sawyer (1989) and Cutler et al. (1990).

The damping rate is

$`\dot{E}_{bulk}`$ $`=`$ $`{\displaystyle d^3x\zeta \omega ^2|\mathbf{}𝝃|^2}.`$ (63)
The Lagrangian compression is

$`\mathbf{}𝝃`$ $`=`$ $`{\displaystyle \frac{g\xi ^r}{c^2}}{\displaystyle \frac{\delta p}{\mathrm{\Gamma }_1p}}{\displaystyle \frac{\xi ^r}{H}}`$ (64)
where the second equality is for low frequency modes. The bulk viscosity coefficient is

$`\zeta `$ $`=`$ $`\zeta _{fid}\omega ^2T_9^6\rho ^2`$ (65)
where

$`\zeta _{fid}`$ $`=`$ $`(6\times 10^{25}\mathrm{g}\mathrm{cm}^1\mathrm{sec}^3)(10^{15}\mathrm{g}\mathrm{cm}^3)^2`$ (66)
$`=`$ $`6\times 10^5\mathrm{g}^1\mathrm{cm}^5\mathrm{sec}^3.`$
Plugging in gives

$`\dot{E}_{bulk}`$ $`=`$ $`\zeta _{fid}T_9^6{\displaystyle d^3x\frac{g^2\rho ^4|\xi _r|^2}{\mathrm{\Gamma }_1^2p^2}}.`$ (67)
We will evaluate this integral for a polytrope of index 1. In this case

$`{\displaystyle \frac{\mathrm{\Gamma }_1p}{\rho ^2}}`$ $`=`$ $`{\displaystyle \frac{GM}{r_{}\rho _0}}`$ (68)
is a constant so that

$`\dot{E}_{bulk}`$ $`=`$ $`\zeta _{fid}T_9^6\left({\displaystyle \frac{\rho _0r_{}}{GM}}\right)^2{\displaystyle d^3xg^2|\xi _r|^2}.`$ (69)
In the WKB limit, this integral is logarithmically divergent at $`r=r_{}`$ and $`\mathrm{cos}\theta =\pm |\mu |`$. This divergence implies that equal contributions to the integrand come per decade of distance from the surface. Since the true eigenfunctions flatten off one wavelength from the surface, we cut off the integrals at this distance. Plugging everything into the integral and approximating slowly varying quantities by their surface values gives the amplitude damping rate

$`\gamma _{bulk}`$ $`=`$ $`{\displaystyle \frac{\dot{E}_{bulk}}{2E_{\mathrm{unit}}}}={\displaystyle \frac{1}{8}}\zeta _{fid}T_9^6{\displaystyle \frac{\rho _0}{r_{}^2\mathrm{\Omega }^2}}{\displaystyle \frac{\mathrm{ln}\mathrm{\Lambda }}{\mu ^2}},`$ (70)
where $`\mathrm{\Lambda }=2|\mu |(1\mu ^2)^{1/2}p`$ is roughly the number of nodes along the rotation axis. Evaluating this expression for a fiducial neutron star with polytrope index $`n=1`$, mass $`M=1.4M_{}`$ and radius $`r_{}=12\mathrm{k}\mathrm{m}`$ we find the numerical value

$`\gamma _{bulk}`$ $`=`$ $`1.7\times 10^{10}\mathrm{sec}^1T_9^6\nu _{\mathrm{khz}}^2{\displaystyle \frac{\mathrm{ln}\mathrm{\Lambda }}{\mu ^2}}.`$ (71)
Note the extremely important fact that this damping rate is very weakly dependent on the wavelength of the mode! The usual picture of a cascade to small scales does not make sense for damping by bulk viscosity. Instead one must carry the energy to small frequency.

For the $`(n,k,m)=(3,1,2)`$ r-mode, the previously calculated value is (Lindblom et al. 1999)

$`\gamma _{bulk}(rmode)`$ $`=`$ $`2.8\times 10^{12}\mathrm{sec}^1T_9^6\nu _{\mathrm{khz}}^2.`$ (72)
The r-mode has a different scaling with $`\mathrm{\Omega }`$ and normalization since the compression is smaller: $`\mathbf{}𝝃\mathrm{\Omega }^2`$ instead of $`\mathbf{}𝝃\xi ^r/H`$.

### 5.3  shear viscosity

The shear viscosity for nuclear matter has been calculated by Flowers & Itoh (1979), with an analytic fit by Cutler et al. (1990) of the form

$`\nu _s`$ $`=`$ $`\nu _{s,fid}(\rho /\rho _0)^{5/4}T_9^2`$ (73)
where $`\nu _{s,fid}=2000\mathrm{cm}^2\mathrm{sec}^1(\rho _0/10^{15}\mathrm{g}\mathrm{cm}^3)^{5/4}`$.

The shear viscosity damping is then

$`\dot{E}_{shear}`$ $`=`$ $`T_9^2{\displaystyle d^3x\rho \nu _s\omega ^2\left(|\xi _{(i,j)}|^2\frac{1}{3}|\mathbf{}𝝃|^2\right)}`$ (74)
$``$ $`T_9^2\nu _{s,fid}\rho _0\omega ^2{\displaystyle d^3x(\rho /\rho _0)^{9/4}k^2|𝝃|^2}`$
where we have kept terms of leading order in $`(|𝒌|H)^1`$, and subscripted brackets denote a symmetrized derivative. Plugging everything in, and approximating the density by a power law with depth appropriate for a polytrope of index 1, we find the damping rate

$`\gamma _{shear}`$ $``$ $`0.6\pi {\displaystyle \frac{p^2}{1\mu ^2}}{\displaystyle \frac{\nu _{s,fid}}{r_{}^2}}T_9^2.`$ (75)
For our fiducial star this becomes

$`\gamma _{shear}`$ $`=`$ $`3.8\times 10^9\mathrm{sec}^1T_9^2{\displaystyle \frac{p^2}{1\mu ^2}}.`$ (76)
The previously computed r-mode shear damping rate is (Lockitch & Friedman 1999)

$`\gamma _{shear}(rmode)`$ $`=`$ $`4\times 10^9\mathrm{sec}^1T_9^2,`$ (77)
which is about a factor of two different from our formula.

As first noted by Bildsten & Ushomirsky (2000), the r-mode is damped much more efficiently by shear in the crust-core boundary layer than by shear over the bulk of the stellar interior. Levin & Ushomirsky (2001) later corrected this damping rate to account for crust with a finite shear modulus. The key parameter is the fractional velocity jump over the boundary layer, called $`\eta `$. Levin and Ushomirsky found the rate of damping to be

$`\gamma _{vbl}(rmode)`$ $`=`$ $`1.5\times 10^3\mathrm{sec}^1\eta ^2\nu _{\mathrm{khz}}^{1/2}T_9^1`$ (78)
with a realistic estimate for the fractional velocity jump of $`\eta 0.1`$. Inclusion of the finite shear modulus of the crust gives much better agreement of the r-mode instability curve with the observations of LMXB’s.

We have neglected damping of the daughter modes by shear in the boundary layer.

## 6  r-mode saturation by discrete modes: the small driving limit

A fundamental plot for the r-mode instability is given in fig.5. The r-mode is unstable for spin frequencies above the thick dashed lines, where $`\gamma _{gr}=\gamma _{shear}`$, $`\gamma _{vbl}`$, or $`\gamma _{bulk}`$. The solid lines show where driving of the r-mode equals damping of daughter modes, indicating marginal stability of the energy transfer. For bulk viscosity, the ratio of driving to damping is

$`\zeta _b`$ $`=`$ $`{\displaystyle \frac{\gamma _{gr}}{\gamma _{bulk}}}=8\times 10^6\nu _{\mathrm{khz}}^8T_9^6`$ (79)
while for shear viscosity

$`\zeta _s`$ $`=`$ $`{\displaystyle \frac{\gamma _{gr}}{\gamma _{shear}}}=1.3\times 10^7\nu _{\mathrm{khz}}^6T_9^2n^2.`$ (80)
In these estimates we have used $`\mu =1/6`$, appropriate for daughter modes with the largest coupling coefficients, and $`n`$ denotes the wavenumber of the daughter mode. Only in the region from the $`\zeta _s=1`$ and $`\zeta _b=1`$ lines to the r-mode instability curve can we possibly have stable energy transfer for the three mode system.

### 6.1  young neutron stars

Nascent, rapidly rotating neutron stars cool into the region of instability (Owen et al. 1998) at fixed spin frequency. For daughter modes mainly damped by bulk viscosity, there is a narrow region near the instability curve in which energy transfer for a single triplet of modes would be stable. However, in this region, the damping is relatively independent of the daughter mode wavenumber. The quality factor of a daughter mode is roughly

$`Q_d`$ $``$ $`{\displaystyle \frac{\omega }{\gamma _{bulk}}}3.4\times 10^5\nu _{\mathrm{khz}}^3T_{10}^6`$ (81)
where we have used the daughter modes with the largest coupling coefficients so that $`\mu _{2,3}=1/6`$. We have also set $`\mathrm{ln}\mathrm{\Lambda }1`$. We found the coupling coefficients are roughly $`\kappa 27`$ for a parent mode with $`n_1=3`$ and $`\mu _1=1/3`$ so that the saturation amplitude for a three mode system is given by

$`A_1^2(bulk)`$ $`=`$ $`{\displaystyle \frac{E(rmode)}{0.5Mr_{}^2\mathrm{\Omega }^2}}{\displaystyle \frac{1}{\kappa ^2Q_d^2}}`$ (82)
$``$ $`10^{14}\nu _{\mathrm{khz}}^6T_{10}^{12}.`$
This formula would imply that nascent neutron stars cooling into the instability curve after a supernova will saturate at a very small fraction of the rotational energy of the star.

However, this formula is not applicable for the following reason. Since the damping rate of the daughter modes is relatively independent of the daughter mode wavenumber, all daughter modes have essentially the same parametric threshold (roughly eq.82) until $`n`$ becomes large enough that shear viscosity becomes comparable to bulk viscosity (see fig.5). We estimate this point to be at $`n_{b=s}10^4T_{10}^4\nu _{\mathrm{khz}}^1`$. For daughter modes with frequency $`\mu 1/6`$, there are roughly $`n_{b=s}^210^8T_{10}^8\nu _{\mathrm{khz}}^2`$ daughter modes parametrically excited to large amplitude by the parent r-mode, so that the discrete limit is not applicable. Thus the r-mode instability in young neutron stars is in the continuum limit discussed in section 7.

### 6.2  LMXB’s

For the LMXB case, neutron stars with temperature $`T3\times 10^8K`$ spin up until they hit the instability curve (Bildsten 1998; Levin 1999). When the instability curve near $`T_9=0.3`$ is set by boundary layer shear viscosity with $`\eta =0.1`$ (Levin & Ushomirsky 2001), we see that if the star stays close to the instability curve, one must go to daughter modes with $`10100`$ nodes before the energy transfer can become stable. As a result, $`10^{24}`$ daughter modes will be parametrically excited to large amplitude, and the continuum limit is more appropriate.

If, however, boundary layer shear viscosity is not operating for some reason, then the discrete mode approximation will be valid near the instability curve. The quality factor of the daughter modes in this case is

$`Q_d`$ $``$ $`5.5\times 10^9\nu _{\mathrm{khz}}T_8^2n^2`$ (83)
giving a saturation amplitude for LMXB’s near the instability curve to be

$`A_1^2(shear)`$ $`=`$ $`{\displaystyle \frac{E(rmode)}{0.5Mr_{}^2\mathrm{\Omega }^2}}{\displaystyle \frac{1}{\kappa ^2Q_d^2}}`$ (84)
$``$ $`10^{21}(0.33/\nu _{\mathrm{khz}})^2T_8^4n^4,`$
which is quite small.

The conclusion we draw in this section is that, for the likely scenario in which either bulk viscosity or boundary layer shear viscosity sets the r-mode instability curve, many modes will be parametrically excited to large amplitude, and the continuum limit discussed in the next section is a better approximation.

## 7  r-mode saturation in the continuum limit

The saturated r-mode energy was found in section 2.2 to be

$`{\displaystyle \frac{E(rmode)}{0.5Mr_{}^2\mathrm{\Omega }^2}}`$ $`=`$ $`10^6\left({\displaystyle \frac{\alpha _e}{0.1}}\right)\nu _{\mathrm{khz}}^5`$ (85)
This solution is valid when there is a clear separation between the inner and outer scales of the turbulence (see Fig. 2). The outer scale is given by the r-mode itself, while the inner scale is where $`\gamma _{\mathrm{nl}}`$, the characteristic rate for amplitude change by nonlinear interactions, is equal to the dissipation rate given by $`\gamma _{\mathrm{shear}}`$ (bulk viscosity is irrelevant for the inner scale; see below). In other words, the inner scale is where the Reynolds number for that scale becomes order unity. We can estimate the nonlinear interaction rate using eqs.16 and 17, with the scalings $`n=n_\alpha n_\beta n_\gamma `$, $`\mu =\mu _\alpha \mu _\beta \mu _\gamma `$, etc. We find

$`\gamma _{\mathrm{nl}}`$ $``$ $`{\displaystyle \frac{I}{N}}{\displaystyle \frac{1}{N}}{\displaystyle \frac{n^5}{\mu }}(\mathrm{\Omega }N)^2\gamma _{\mathrm{gr}}{\displaystyle \frac{n}{n_r}}\left({\displaystyle \frac{\mu _r}{\mu }}\right)^{3/2}`$ (86)
where $`n_r`$ and $`\mu _r`$ set the scale for the driving region at which $`\gamma _{\mathrm{gr}}=\gamma _{\mathrm{nl}}`$.

The expression for shear viscosity from eq.76 can be written $`\gamma _{\mathrm{shear}}=\gamma _sn^2`$ in the $`\mu ^21`$ limit. Equating $`\gamma _{\mathrm{nl}}`$ to $`\gamma _{\mathrm{shear}}`$, we find the inner scale is

$`{\displaystyle \frac{n}{n_r}}\left({\displaystyle \frac{\mu }{\mu _r}}\right)^{3/2}`$ $``$ $`n_r^2{\displaystyle \frac{\gamma _{\mathrm{gr}}}{\gamma _s}}`$ (87)
$``$ $`170n_r^2\left({\displaystyle \frac{\nu _{\mathrm{spin}}}{330\mathrm{Hz}}}\right)^6T_8^2.`$
In the region of the $`(\mathrm{\Omega },T)`$ plane where the right hand side of Eq. (87) is large, many modes are excited with a clear separation between inner and outer scales of turbulence (see Fig. 2). This region includes the entire instability window when the r-mode is damped by a viscous boundary layer. When there is no viscous boundary layer, there is a small region close to the instability curve where (87) is small and where the discrete limit applies instead of the continuum limit.

For bulk viscosity things work differently. The cascade solution will be valid in the region of phase space where $`\gamma _{\mathrm{nl}}\mathrm{¿}\mathrm{}\gamma _{\mathrm{bulk}}`$, or

$`{\displaystyle \frac{\mu }{\mu _r}}\left({\displaystyle \frac{n}{n_r}}\right)^2`$ $`\mathrm{¿}\mathrm{}`$ $`\mu _r^4\left({\displaystyle \frac{\gamma _b}{\gamma _{\mathrm{gr}}}}\right)^2`$ (88)
$`=`$ $`10^5\mu _r^4T_{10}^{12}\nu _{\mathrm{khz}}^{16},`$
where we have written $`\gamma _{\mathrm{bulk}}=\gamma _b\mu ^2`$. In the region of the $`(\mathrm{\Omega },T)`$ plane where the right hand side of Eq. (88) is large compared to unity, bulk viscosity is dominant at the outer scale and no cascade solution exists. When the right hand side of Eq. (88) is small compared to unity, then a cascade can form, but the bulk viscosity is irrelevant for setting the inner scale of the cascade. We note that the boundary (88) of the bulk viscosity dominated regime approximately coincides with the curve $`\zeta _b=1`$ in Fig. 5. A newly born neutron star will very rapidly move from the instability curve to the $`\zeta _b=1`$ curve at which point the cascade can form.

Finally we note 161616 We thank P. Goldreich for bringing this to our attention. that the weak turbulence approximation which underlies the derivation of Eq. (16) eventually breaks down as one goes to small scales. The weak turbulence approximation requires that the nonlinear energy transfer timescale $`1/\gamma _{\mathrm{nl}}20\mathrm{s}\nu _{\mathrm{khz}}^6n^1\mu ^{3/2}`$ be much longer than the mode period $`5\times 10^4\mathrm{s}\mu ^1\nu _{\mathrm{khz}}^1`$, which breaks down in the regime

$$n\mu ^{5/2}\mathrm{¿}\mathrm{}4\times 10^4\nu _{\mathrm{khz}}^5.$$
(89)
Thus our turbulent cascade solution of the equations of motion (16) will likely be replaced by some form of strong turbulence at sufficiently small scales. However, this should not affect our prediction of the r-mode’s saturation amplitude, as the regime (89) in phase space where the approximation breaks down is well separated from the driving regime $`\mu n1`$.

## 8 comparison with previous work

There have been three distinct alternative nonlinear mechanisms proposed to saturate the growth of the r-mode: (1) For large amplitude pulsations, the Fermi energies of the electron- proton-neutron gas become significantly shifted, and the kinetic energy can be rapidly converted to both heat and neutrinos by nonlinear bulk viscosity (Reisenegger 2001). (2) The amplitude grows so large ($`EE_{\mathrm{rotation}}GM^2/r_{}`$) that strong shocks occur, rapidly thermalizing the kinetic energy (Lindblom et al. 2001, 2002); (3) In neutron stars with a crust, a turbulent boundary layer forms at the crust-core interface (Wu et al. 2001). We now discuss each in a bit more detail.

Since the Fermi energy of the electron and neutron have a different dependence on density, the Fermi surfaces are shifted out of beta equilibrium when matter is compressed. The scaling of the resulting neutrino emission rate depends on the ratio of chemical potential imbalance, $`\mathrm{\Delta }\overline{\mu }_\nu =\mu _n\mu _p\mu _e`$, to temperature, which is roughly (Reisenegger 1995)

$`{\displaystyle \frac{\mathrm{\Delta }\overline{\mu }_\nu }{T}}`$ $``$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{E_{F,e}}{T}}{\displaystyle \frac{\mathrm{\Delta }\rho }{\rho }}`$ (90)
where $`E_{F,e}`$ is the Fermi energy of the electron. When this ratio is large, $`5/8`$ of the resulting dissipation heats the star and $`3/8`$ goes into neutrinos. The rate of such dissipation scales as $`\left(\frac{\mathrm{\Delta }\overline{\mu }_\nu /T}{7.33}\right)^8(\mathrm{\Delta }\rho /\rho )^8`$ times the neutrino emissivity of uncompressed matter. The mode damping rate is extremely sensitive to the compression, and can saturate the growth of the r-mode for sufficiently large amplitude. Reisenegger (2001) has done a detailed calculation, finding that the saturation energy is comparable to the stellar rotation energy. This interesting idea gives a larger (less constraining) saturation amplitude compared to the value found in this paper.

Next, Lindblom et al. (2001, 2002) have performed state-of-the-art 3D Newtonian hydrodynamics simulations including a prescription for the radiation reaction force. The only damping mechanism included in the code is numerical viscosity, and of order $`128^3`$ points were used. They were able to follow the linear growth of the r-mode, all the way into the nonlinear regime where $`EE_{rot}`$. Shocks then formed near the surface of the star, rapidly thermalizing the kinetic energy of the mode.

Since the growth of the instability is so slow compared to the dynamical time in the star, they found it necessary to artificially increase the radiation reaction force by a factor of $`4500`$. A natural question is how the mode would saturate if the correct, physical value of the driving force was used. The following physical example is useful to consider. Imagine water waves being driven by wind moving at $`1\mathrm{cm}\mathrm{sec}^1`$, a whisper of a breeze, as compared to $`4500\mathrm{cm}\mathrm{sec}^1`$, a hurricane. For small amplitude water waves, four wave interactions can transport energy to small scales, saturating the growth of the waves. In a hurricane, the wave growth time is so short that waves grow to large amplitudes and break. Since Lindblom et al. have not addressed how saturation might occur for physical values of either driving or damping of the waves, the relevance of their simulations to saturation of the r-mode instability is not clear.

One comparison which can be made is to use our formula in eq.24 to estimate the saturation amplitude seen in Lindblom et al.’s simulations when $`\gamma _{gr}4500\gamma _{gr}`$. We find $`E_{\mathrm{simulation}}/(0.5Mr_{}^2\mathrm{\Omega }^2)=0.1\left(\alpha _e/0.1\right)(4500\gamma _{gr}/\mathrm{\Omega })5\times 10^3\left(\alpha _e/0.1\right)\nu _{\mathrm{khz}}^5`$. This result can be translated into Lindblom et al.’s notation by using $`A_1^20.5\stackrel{~}{J}\alpha _{\mathrm{simulation}}^2=0.008\alpha _{\mathrm{simulation}}^2`$ with the result $`\alpha _{\mathrm{simulation}}0.7\left(\alpha _e/0.1\right)^{1/2}\nu _{\mathrm{khz}}^{5/2}`$. This comparison shows that if one attempted to extrapolate the saturation amplitude over more than three decades in driving force, that the saturation amplitude by mode coupling would indeed be of order unity, in their notation. This does not, of course, explain the saturation amplitude seen in Lindblom et al.’s simulations, which they explain is due to strong shocks near the stellar surface. However, the results of this paper show that their claims of (1) saturation energy of order the rotation energy, and (2) strong shocks as the saturation mechanism, are not supported. They are an artifact of the unphysically large value for the radiation reaction force.

We comment further on the ability of a numerical simulation to accurately reproduce the cascade of energy to small scales as derived in this paper. In simulations with $`128^3`$ points, only a certain number of modes exist because of the finite resolution. Since the detuning is a rapidly decreasing function of wavenumber, secular energy transfer by nearly resonant interactions becomes more important as the number of grid points increases. For instance, daughter modes with half the frequency of the r-mode have quantum numbers in the ratio $`k/n1/20`$, so that one needs nearly $`20`$ times as many nodes in the cylindrical radius as in the z direction in order to find parametrically excited daughter modes. We estimate that only a few of these might have been accurately modeled by Lindblom et al’s simulations. In time evolutions of the mode amplitude equations with a small number of low order, nonresonant modes (Morsink 2001; Arras 2001), large saturation amplitudes were found as compared to the results in this paper. The reason, as can be clearly seen in eq.9 for the parametric threshold, is that the r-mode cannot easily excite daughter pairs with large detuning. However, going to higher order modes with much smaller detuning can give a saturation amplitude orders of magnitude smaller than for arbitrary, low order modes.

Lindblom et al. specifically commented that three-mode coupling is not the saturation mechanism in section H of their paper. Their claim was based on the lack of power observed in certain modes besides the r-mode during their simulation. However, they focused on interactions which couple the r-mode twice to a third mode. As they themselves comment at the end of section H, they have not included parametric excitation of daughter modes in their constraints. As discussed in section 2.1, couplings of the type discussed by Lindblom et al. (2002) are far less important than parametric couplings, because (1) they are down by a factor of parent mode amplitude, which is small, and (2) only a relatively small region of phase space couples well with the r-mode by non-parametric couplings. Hence, Lindblom et al.’s constraints are not useful since they constrain an unimportant process.

The inability of simulations to include very high order modes presumably also explains the results of the fully relativistic simulations of Font & Stergioulas (2001), in which an r-mode with order unity amplitude was observed not to lose any energy to other modes over several dynamical times. More recent simulations by Gressman et al. (2002) show that for slightly larger initial amplitudes, the r-mode decays rapidly into a differentially rotating configuration without shocks forming. These results are not inconsistent with our analyses, but our results indicate that the r-mode never reaches the regime of rapid nonlinear decay seen by Gressman et al. (2002).

Next we discuss the turbulent boundary layer mechanism of Wu et al. (2001) which operates in neutron stars with a crust. Energy dissipation by turbulent drag scales as $`A_1^3`$, leading to saturation of the mode. The attractiveness of this idea is that the turbulent drag force is well understood in magnitude and scaling both from numerical estimates as well as laboratory experiments. These authors considered the effect of such energy dissipation on the crust and thermal history of the star, and go on to discuss the observable spin frequency of the star after it exits the r-mode instability region. For a realistic fractional velocity jump across the crust-core boundary layer, $`\eta 0.1`$, they found the r-mode saturated at a value $`E/(0.5Mr_{}^2\mathrm{\Omega }^2)0.2\nu _{\mathrm{khz}}^{10}`$, which is larger (less constraining) than the value found here, both in normalization, and in the dependence on $`\nu _{\mathrm{khz}}`$. Furthermore, their mechanism does not operate in completely fluid stars without a crust, which is the case for hot young neutron stars.

Lastly, we mention that this paper is a companion paper to that of Morsink (2002), which discusses the nonlinear coupling among r-modes in a star for which buoyancy forces are dominant over Coriolis forces. Morsink found that, because the r-mode frequency decreases with $`m`$, interactions do not become more resonant as the daughter mode $`m`$ increase. As a result, energy transfer among three r-modes is not likely to produce a saturation value as low as in this paper.

We conclude that nonlinear mode coupling to inertial modes provides the most stringent constraints on the r-mode amplitude at this time.

## 9  spin evolution of neutron stars

The spindown torque exerted on the neutron star by gravitational radiation is roughly

$`\tau _{gr}`$ $``$ $`0.12\left({\displaystyle \frac{\mathrm{\Omega }r_{}}{c}}\right)^7{\displaystyle \frac{GM^2}{r_{}}}A_1^2`$ (91)
$``$ $`2\times 10^{42}\mathrm{erg}\left({\displaystyle \frac{\alpha _e}{0.1}}\right)\nu _{\mathrm{khz}}^{12}.`$
The spindown time associated with this torque is

$`t_{\mathrm{spindown}}`$ $``$ $`{\displaystyle \frac{Mr_{}^2\mathrm{\Omega }}{11\tau _{gr}}}0.01\mathrm{yr}\left({\displaystyle \frac{\alpha _e}{0.1}}\right)^1\nu _{\mathrm{khz}}^{11}`$ (92)
$``$ $`2\times 10^3\mathrm{yr}\left({\displaystyle \frac{\alpha _e}{0.1}}\right)^1\left({\displaystyle \frac{330\mathrm{Hz}}{\nu }}\right)^{11}.`$
Since the spindown rate decreases strongly with spin frequency, most of the time is spent at lower frequencies.

The spindown time becomes $`\mathrm{¿}\mathrm{}10^4\mathrm{yr}`$ at the lowest rotation rates inside the instability curve, while it is of order a few days for stars rotating near breakup. This is of interest for certain gamma-ray burst models, such as the “supranova” model (Vietri & Stella 1999) in which core collapse leads to ejection of the stellar envelope, as well as a rapidly rotating neutron star which is above the maximum mass for a nonrotating star. Angular momentum transport can then slow the neutron star down, leading to collapse to a black hole and generation of a powerful gamma-ray burst. Our results imply that the gamma-ray burst should occur within of order a week after the supernova explosion.

Next we turn our attention to neutron stars in LMXB’s. The ratio of spindown torque, due to radiation reaction, to accretion torque $`\tau _{acc}\dot{M}(GMr_{})^{1/2}`$ is roughly

$`{\displaystyle \frac{\tau _{gr}}{\tau _{acc}}}`$ $``$ $`{\displaystyle \frac{2\times 10^{42}(\alpha _e/0.1)\nu _{\mathrm{khz}}^{12}}{10^{34}(\dot{M}/10^8M_{}\mathrm{yr}^1)}}`$ (93)
$``$ $`{\displaystyle \frac{10^8M_{}\mathrm{yr}^1}{\dot{M}}}\left({\displaystyle \frac{\alpha _e}{0.1}}\right)\left({\displaystyle \frac{\nu _{\mathrm{spin}}}{200\mathrm{Hz}}}\right)^{12}.`$
For accretion rates smaller than the Eddington rate $`\dot{M}10^8M_{}\mathrm{yr}^1`$, and spin frequencies above $`\nu _{\mathrm{spin}}=200\mathrm{Hz}`$, the radiation reaction torque is larger than the accretion torque and can halt the further spinup of the neutron star. If the neutron star viscosity is dominated by normal matter, then the star enters into a limit cycle of spinup by accretion and spindown by the r-mode, as discussed by Levin (1999). \[The alternative equilibrium scenario is discussed in Sec. 10 below.\] Since the r-mode is only likely to be unstable for $`\nu _{\mathrm{spin}}\mathrm{¿}\mathrm{}300\mathrm{Hz}`$, the r-mode can halt spinup inside the region of instability.

The observable spin frequency is determined by where the star exits the region of r-mode instability, if no other process spins the star down further. The exact spin frequency at which the star exits the region of r-mode instability depends on the evolution of both the spin frequency and the stellar temperature (Levin 1999; Owen et al. 1998; Wu et al. 2001). We can estimate this terminal frequency (Wu et al. 2001) by equating the neutrino cooling luminosity, $`L_\nu =7.4\times 10^{39}T_9^8\mathrm{erg}\mathrm{sec}^1`$, with the rate of stellar heating due to the r-mode. If we approximate that all the energy input to the r-mode by radiation reaction is damped away as heat, the rate of heating of the star is just given by $`\dot{E}_{\mathrm{heat}}=2\gamma _{gr}E`$, where $`E`$ is the saturation energy found in eq.85. Equating heating and cooling, we find the equilibrium temperature as a function of spin frequency, given by

$`T_{\mathrm{eq}}`$ $``$ $`10^9K\left({\displaystyle \frac{\alpha _e}{0.1}}\right)^{1/8}\left({\displaystyle \frac{\nu _{\mathrm{spin}}}{330\mathrm{Hz}}}\right)^{13/8}.`$ (94)
The crystallization temperature of the crust is (Wu et al. 2001) $`T_{\mathrm{melt}}(510)\times 10^9K`$ so the heating by the r-mode cannot prevent the crust from forming when $`\nu _{\mathrm{spin}}10^3\mathrm{Hz}`$. If the instability curve is set by boundary layer shear viscosity ($`\gamma _{gr}=\gamma _{vbl}`$), the intersection of the equilibrium spin down curve with the r-mode instability curve is given by the terminal frequency

$`\nu _{\mathrm{terminal}}`$ $``$ $`250\mathrm{Hz}\left({\displaystyle \frac{\alpha _e}{0.1}}\right)^{0.02}\left({\displaystyle \frac{\eta }{0.1}}\right)^{0.28}`$ (95)
with a core temperature of roughly

$`T_{\mathrm{terminal}}`$ $``$ $`6\times 10^8K\left({\displaystyle \frac{\alpha _e}{0.1}}\right)^{0.15}\left({\displaystyle \frac{\nu _{\mathrm{spin}}}{330\mathrm{Hz}}}\right)^{0.46}.`$ (96)
Note that the observable spin frequency is very insensitive to the saturation parameter $`\alpha _e`$, as well as to the fractional velocity jump $`\eta `$. The spin frequency found in eq.95 is comparable to the lower end of the observed LMXB’s, consistent with a limit cycle (Levin 1999) of spin-up by accretion and spin-down by the r-mode. The timescale to exit the instability curve is roughly $`t_{\mathrm{spindown}}2000\mathrm{yr}(330/250)^{11}(\alpha _e/0.1)^14\times 10^4\mathrm{yr}(\alpha _e/0.1)^1`$. This spindown timescale is very sensitive to the position of the instability curve.

For young neutron stars with strong magnetic fields, the spindown torque from magnetic dipole radiation is comparable to that from gravitational radiation. Equating the magnetic dipole spindown timescale (Shapiro & Teukolsky 1983) $`t_{\mathrm{md}}30\mathrm{yr}B_{12}^2\nu _{\mathrm{khz}}^2`$ to $`t_{\mathrm{spindown}}`$, we find that gravitational radiation reaction dominates for frequencies above $`\nu _{\mathrm{spin}}400\mathrm{Hz}(0.1/\alpha _e)^{1/9}B_{12}^{2/9}`$, where $`B_{12}`$ is the surface dipole field in units of $`10^{12}\mathrm{G}`$. Hence for typical pulsars with magnetic fields $`10^{12}\mathrm{G}`$, the spindown torque is dominated by the r-mode only for fairly large spin frequencies.

## 10 Detectability of gravitational waves

We now discuss the prospect of detecting gravitational waves from r-modes, based on the saturation amplitude (25). We consider three different scenarios: (i) newly born neutron stars where an optically observed extra-Galactic supernova provides the sky location for the gravitational wave search; (ii) LMXB’s in the spinup-spindown limit cycle first discussed by Levin (1999); and (iii) LMXB’s in spin and thermal equilibrium.

For newly born neutron stars, Brady & Creighton (2000) (BC) discuss the detection likelihood by LIGO assuming a large saturation amplitude. They parameterize the saturation amplitude in terms of a parameter $`\kappa `$ in their Eq. (7.3), which we denote by $`\kappa _{bc}`$. Our result (25) gives $`\kappa _{bc}1.2\times 10^3\nu _{\mathrm{khz}}^5`$, while BC took $`\kappa _{bc}=1`$. In the first year of spindown, $`\nu _{\mathrm{khz}}`$ decreases from $`1`$ to $`0.66`$ \[cf. Eq. (92) above\], and thus the gravitational wave strain amplitude will be a factor $`\sqrt{\kappa _{bc}}2.8\times 10^2`$ smaller than that considered by BC. The distance to which the source can be seen by enhanced LIGO detectors, for fixed integration time (see below), is correspondingly reduced from BC’s estimate of $`8\mathrm{Mpc}`$ to $`200\mathrm{kpc}`$, almost inside the Galaxy. Since the galactic supernova rate is roughly once per $`50100\mathrm{yrs}`$, the probability that LIGO will detect young neutron stars radiating due to r-modes is small.

We now discuss why we can treat the integration time as fixed. The matched filtering signal to noise ratio $`S/N`$ for gravitational waves, when averaged over source orientations and polarizations, depends only on the energy per unit frequency $`dE/df`$ of the waves (Flanagan & Hughes 1998):

$$\frac{S^2}{N^2}=\frac{2G}{5\pi ^2c^3D^2}𝑑f\frac{1}{f^2S_h(f)}\frac{dE}{df}.$$
(97)
Here $`D`$ is the distance to the source and $`S_h(f)`$ is the detector noise spectrum. For waves of fixed azimuthal quantum number $`m`$, using the replacement $`dE=2\pi fdJ/m`$ yields (Blandford 1984; Lindblom & Owen 2002b)

$$\frac{S^2}{N^2}=\frac{4G}{5\pi mc^3D^2}𝑑f\frac{1}{fS_h(f)}\frac{dJ}{df},$$
(98)
where $`J`$ is the $`z`$-component of angular momentum. As noted by Lindblom & Owen (2002b), the expression (98) is independent of how quickly the star looses angular momentum, and hence of the saturation amplitude. Thus, a priori one would not expect our low saturation amplitude (25) to affect very much the detectability of the signal. The problem however is that it is not possible to integrate long enough to accumulate the total signal-to-noise ratio (98).

Using the stellar model discussed before Eq. (62), the relation $`f=4\nu _{\mathrm{spin}}/3`$ between gravitational wave frequency $`f`$ and spin frequency, the broadband LIGO-II noise curve 171717In the relevant frequency range $`f\mathrm{¿}\mathrm{}500\mathrm{Hz}`$, this noise curve is approximately given by $`fS_h(f)1.7\times 10^{44}(f/1000\mathrm{Hz})^3`$. from Gustafson et al. (1999), and neglecting the spin dependence of the moment of inertial of the star, we can evaluate (98) for a spindown from an initial spin frequency $`\nu _{\mathrm{khz},\mathrm{i}}`$ in kHz to a final spin frequency $`\nu _{\mathrm{khz},\mathrm{f}}`$. The result is

$$\frac{S}{N}=5.4\left(\frac{10\mathrm{M}\mathrm{p}\mathrm{c}}{D}\right)\left(\frac{1}{\nu _{\mathrm{khz},\mathrm{f}}^2}\frac{1}{\nu _{\mathrm{khz},\mathrm{i}}^2}\right)^{1/2}.$$
(99)
The complete spindown from say $`1\mathrm{kHz}`$ to $`250\mathrm{H}\mathrm{z}`$ \[cf. Eq. (95) above\] gives $`S/N21`$ at $`10\mathrm{Mpc}`$. The first year of spindown from $`1000\mathrm{Hz}`$ to $`650\mathrm{Hz}`$ \[cf. Eq. (92) above\] gives instead $`S/N6.2`$, which is not too much smaller.

However, the need to perform a search over spindown parameters in practice limits the integration time to $`10^6`$ seconds. BC analyzed the performance of a “stack slide” search method, involving demodulating the signal for many different choices of spindown parameters, dividing the demodulated data into several chunks or “stacks”, computing the power spectrum of each stack and adding the power spectra. The threshold value $`\rho _{\mathrm{th}}`$ of $`S/N`$ for this method, assuming $`1\%`$ false alarm probability, is approximately given by solving the equation

$$\mathrm{\Gamma }(N_s,N_s+\rho _{\mathrm{th}}^2/4)/\mathrm{\Gamma }(N_s)=0.01/(N_bN_p).$$
(100)
where $`\mathrm{\Gamma }(a,b)`$ is the incomplete gamma function, $`N_s`$ is the number of stacks, $`N_b`$ is the number of frequency bins per stack, and $`N_p`$ is the number of points in the space of spindown parameters. BC’s estimate that r-modes are detectable out to $`8\mathrm{Mpc}`$ was based on assuming a Teraflop of computing power, which implied an optimum detection strategy of $`8`$ stacks of $`10^5`$ s duration each, integrating from $`\nu _{\mathrm{spin}}=200\mathrm{Hz}`$ to $`\nu _{\mathrm{spin}}186\mathrm{Hz}`$.

We can modify the BC analysis for our turbulent cascade scenario as follows. Optimum sensitivity is achieved late in the spindown, so we assume that $`\nu _{\mathrm{spin}}650\mathrm{Hz}`$, corresponding to 1 year after the start of the spindown if the initial spin frequency is 1 kHz. We take the parameter values $`\mu _{\mathrm{max}}=0.3`$, $`f_{\mathrm{max}}=(4/3)650\mathrm{H}\mathrm{z}`$ (instead of $`200\mathrm{Hz}`$ as in BC), and a spindown timescale $`\tau _{\mathrm{min}}=1\mathrm{yr}`$, which from Eq. (92) is appropriate after 1 year of spindown. Maximizing over the number of stacks and stack durations as in BC gives that the optimum detection strategy for a Teraflop of computing power is to use $`10`$ stacks of duration $`3\times 10^4`$ s each. The corresponding number of parameter space points is $`N_p6\times 10^8`$, from Eq. (2.20) of BC, which gives from Eq. (100) a threshold value of $`\rho _{\mathrm{th}}15`$. Combining this with Eq. (99), and noting that at $`\nu _{\mathrm{spin}}650\mathrm{H}\mathrm{z}`$ an integration time of $`3\times 10^5`$ s corresponds to $`\mathrm{\Delta }\nu _{\mathrm{spin}}0.5(\alpha _e/0.1)\mathrm{Hz}`$ gives that the source would be detectable to $`200(\alpha _e/0.1)^{1/2}\mathrm{kpc}`$, consistent with our earlier estimate181818 The main reason for the loss of sensitivity compared to BC is the reduction in $`\mathrm{\Delta }\nu _{\mathrm{spin}}`$ from $`14\mathrm{Hz}`$ to $`0.5\mathrm{Hz}`$; the star is spinning down more slowly..

This conclusion, however, is based on the assumption of using the stack-slide search method. It is conceivable that an alternative signal processing strategy (and increased computational power) might enable one to integrate for longer periods and achieve a sensitivity closer to the original BC estimate.

We mention in passing another possible difficulty in searching for the signal from r-modes when a turbulent cascade is present. This difficulty is that the phase of the r-mode will wander randomly in time due to the interaction with the turbulent cascade, on some timescale $`t_\mathrm{c}`$. The peak in the Fourier transform of the demodulated data stream will correspondingly be smeared out over a frequency interval of width $`1/t_\mathrm{c}`$, which will be over several frequency bins if the stack size is larger than $`t_\mathrm{c}`$.

The phase coherence timescale for a typical mode in the cascade will be of order

$$t_c\frac{1}{\gamma _{\mathrm{nl}}},$$
(101)
or smaller, where $`\gamma _{\mathrm{nl}}`$ is the the nonlinear energy transfer rate (86) (Zakharov et al. 1992). For the r-mode this is only $`200\mathrm{s}`$ at $`\nu _{\mathrm{spin}}700\mathrm{H}\mathrm{z}`$. However, one might expect the coherence time of the r-mode to be somewhat longer than the estimate (101), since the r-mode is being pumped coherently and is loosing energy by interacting simultaneously with a large number of different modes. Unless the phase coherence time for the r-mode is $`10^210^3`$ times larger than the estimate (101), the sensitivity of the search will be reduced. Again, it may be possible to modify the data analysis procedure to compensate for the phase wandering 191919The method suggested by BC to compensate for phase wandering requires a stack size shorter than $`t_c`$ and a computational power that scales as $`3^{T/t_c}`$, where $`T`$ is the total integration time. In practice this limits $`T`$ to $`\mathrm{¡}\mathrm{}(2030)t_c`$..

Next, we consider the detectability of r-modes in LMXB’s in the spin up/spin down limit cycle. The signal from the spin down phase is essentially the same as for newborn neutron stars, except that they will typically be seen at a low frequency where most of the spindown time is spent. At lower spin frequencies the search over spin-down parameters becomes significantly easier, since the spindown timescale is longer. The formula (4.3) in BC for computational power, with $`f_{\mathrm{max}}=\nu _{\mathrm{spin}}`$ and $`\tau _{\mathrm{min}}`$ given by Eq. (92), shows that for $`\nu _{\mathrm{spin}}\mathrm{¡}\mathrm{}400\mathrm{Hz}`$ integration times as long as $`10^7`$ s can be achieved with 1 Teraflop of computing power. Combining Eqs. (92) and (99) gives that the signal to noise ratio for a $`10^7`$ s integration starting at $`\nu _{\mathrm{spin}}`$ is

$$\frac{S}{N}9\left(\frac{\nu _{\mathrm{spin}}}{330\mathrm{Hz}}\right)^{9/2}\left(\frac{100\mathrm{kpc}}{D}\right)\left(\frac{\alpha _e}{0.1}\right)^{1/2}.$$
(102)
In the regime $`\nu _{\mathrm{spin}}400\mathrm{H}\mathrm{z}`$, the signal-to-noise threshold from the BC method is $`\rho _{\mathrm{th}}10`$, within a factor of $`2`$, giving that the signal should be visible to a distance

$$D90\mathrm{kpc}\left(\frac{\nu _{\mathrm{spin}}}{330\mathrm{Hz}}\right)^{9/2}\left(\frac{\alpha _e}{0.1}\right)^{1/2}$$
(103)
for $`\nu _{\mathrm{spin}}400\mathrm{Hz}`$.

However, as noted by Levin (1999), the chance of observing a particular source emitting gravitational waves is proportional to the relative length of time spent in the spindown phase of the limit cycle. Using our saturation amplitude we find a duty cycle $`10^3(0.1/\alpha _e)`$, implying that one would need of order $`10^3(\alpha _e/0.1)`$ LMXB’s within the distance (103) in order to overcome the small duty cycle. Nevertheless, as argued by Heyl (2002), there may be enough Galactic LMXB’s that some will be seen in the spin-down phase by enhanced LIGO, especially if $`\alpha _e`$ is smaller than $`0.1`$.

We note that for LMXB’s, the phase wandering of the r-mode due to the turbulent cascade is less of a problem, since the nonlinear energy transfer timescale (86) increases rapidly as $`\nu _{\mathrm{spin}}`$ decreases. There is in addition a phase wandering due to fluctuations in the accretion torque, but this occurs over much longer timescales and can be dealt with in the manner suggested by BC.

The third possibility we consider is when the viscosity of an accreting star is independent of temperature or is an increasing function of temperature. In such a case the star can achieve an equilibrium state where the accretion spin-up torque is stably balanced by the radiation reaction torque due to the r-mode, and r-mode heating is balanced by neutrino cooling (Levin 1999). Such equilibria have been found for stars with hyperon cores (Wagoner 2002) and for strange stars (Andersson, Jones & Kokkotas 2001). \[However, the central densities of neutron stars are sufficiently uncertain that hyperon cores may or may not exist.\] In these scenarios, the equilibrium r-mode amplitude is not set by the turbulent cascade considered here, but instead by the equilibrium conditions. The gravitational wave signal is weaker than the limit cycle case considered above, and its strength can be inferred from the X-ray flux (Bildsten 1998).

## 11  conclusions

In this paper, we have accomplished several objectives, which can be divided into stellar oscillation theory, and phenomenology of neutron star spin evolution.

We have, for the first time, presented a WKB theory of global stellar inertial modes (section 3 and appendix B), including both the rapidly varying phase and amplitude which rises quickly toward the surface. Both the eigenmodes and pulsation frequencies take on a very simple form, which was never clearly elucidated in previous calculations \[e.g., Lindblom & Ipser (1999); Lockitch & Friedman (1999).\] Appendix B gives a detailed mathematical treatment of inertial waves. We have estimated when the affects of buoyancy become important. The damping rates by bulk and shear viscosity appropriate for neutron stars have been derived in section 5, and reduced to simple, accurate formulae giving the scalings with stellar parameters and mode quantum numbers.

Next, we have given a complete review of saturation of an overstable mode in the two different limits of strong and weak driving force. The literature for the weak driving limit, familiar from studies of main sequence or white dwarf pulsators (Dziembowski & Krolikowska 1985; Wu & Goldreich 2001), is reviewed in detail, and shown not to apply to most physical situations in which the r-mode instability operates. The strong driving limit, in which a turbulent cascade forms, has never been applied to stellar oscillations to our knowledge. Therefore, the weak turbulence methods in this paper may find application for amplitude saturation in stars with a driving force strong enough to parametrically excite many modes, but weak enough that shocks do not form. In the strong driving limit, we find that the r-mode saturates at an energy $`E/E_{rotation}\gamma _{gr}/\mathrm{\Omega }`$.

The consequences of these calculations for neutron star spin evolution are as follows. First, the time scale for a rapidly rotating ($`10^3\mathrm{Hz}`$) young neutron star to spin down to a frequency $`\nu _{\mathrm{spin}}`$ by gravitational radiation spindown torque is $`t_{\mathrm{spindown}}2\times 10^3\mathrm{yr}(\alpha _e/0.1)^1(\nu _{\mathrm{spin}}/330\mathrm{Hz})^{11}`$. Hence, the r-mode can be responsible for an initial, rapid spindown, but magnetic dipole radiation will dominate for spin frequencies below $`\nu _{\mathrm{spin}}400\mathrm{Hz}(0.1/\alpha _e)^{1/9}B_{12}^{2/9}`$. Second, in spite of the small saturation amplitude, we find that the gravitational radiation spin-down torque is still sufficiently large to halt the spin-up by accretion for neutron stars in LMXB’s. Hence, our calculation confirms the validity of this assumption by previous investigators. If the viscosity is dominated by normal matter, the star will enter into the limit cycle in spin and temperature discussed by Levin (1999). Finally, we estimate that newly born neutron stars will be visible to $`200(\alpha _e/0.1)^{1/2}\mathrm{kpc}`$ with enhanced LIGO interferometers, and LMXB’s in the spin down phase of the Levin limit cycle out to $`90(\nu _{\mathrm{spin}}/330\mathrm{Hz})^{9/2}(\alpha _e/0.1)^{1/2}\mathrm{kpc}`$ for $`\nu _{\mathrm{spin}}\mathrm{¡}\mathrm{}400\mathrm{Hz}`$.

It is a pleasure to acknowledge many useful conversations on stellar oscillations with Yanqin Wu. P.A. would also like to thank Chris Matzner, Chris Thompson, and Maxim Lyutikov. Part of this work was completed when two of us (P.A. and S.M) were at the Institute for Theoretical Physics at the Spin and Magnetism in Young Neutron Stars workshop. We thank Lars Bildsten for his gracious hospitality and for a number of useful conversations, and John Friedman, Curt Cutler and Peter Goldreich for useful comments on the manuscript. P.A. is supported by an NSERC fellowship. This work was supported in part by NSF grants PHY-9900672 and PHY-0084729 at Cornell University. E.E.F. was supported by NSF grants PHY-9722189 and PHY-0140209 and by the Alfred P. Sloan Foundation. S.M. received support from NSERC.

## Appendix A cascade solution

In this appendix we find scale-free solutions to the kinetic equation 17. The material in this appendix can be found in texts on weak turbulence theory, e.g. Zakharov et al. (1992), section 3.3. For convenience we use dimensionless units with $`E_{\mathrm{unit}}=1`$ and $`\mathrm{\Omega }=1`$.

We approximate the mode frequency as $`\omega =2\mu 2\pi k/n=s2\pi |k|/n`$, where $`s=\pm 1`$ is the sign of the frequency. For notational simplicity, we will use both $`k`$ and $`\mu `$ as positive numbers for the following derivation, using $`s`$ to take into account the sign. The resulting expression will then be written in a form valid for either positive or negative $`\mu `$. The sum over modes is then given by

$`{\displaystyle \underset{\alpha }{}}`$ $``$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=n/\pi }{\overset{n/\pi }{}}}{\displaystyle \underset{m=n}{\overset{n}{}}}\pi ^1{\displaystyle \underset{s=\pm 1}{}}{\displaystyle _0^{\mathrm{}}}𝑑nn{\displaystyle _0^1}𝑑\mu {\displaystyle _n^n}𝑑m.`$ (A1)
Here we used $`\mu =\pi k/n`$ instead of $`k`$ since $`k`$ has an implicit scaling with $`n`$. The coupling coefficients are assumed to be negligible if momentum conservation is not satisfied, so that

$`|\kappa _{\alpha \beta \gamma }|^2`$ $``$ $`|\overline{\kappa }_{\alpha \beta \gamma }|^2\delta (s_{\alpha 1}n_\alpha +s_{\beta 1}n_\beta +s_{\gamma 1}n_\gamma )\delta (s_{\alpha 2}k_\alpha +s_{\beta 2}k_\beta +s_{\gamma 2}k_\gamma )\delta (m_\alpha +m_\beta +m_\gamma ).`$ (A2)
The signs $`s_{\alpha 2}`$ etc. allow for waves moving in either direction (see appendix B). We use the fact that the coupling coefficient is approximately independent of both $`m`$ and the sign of $`\mu `$, and is separately scale invariant in $`n`$ and $`\mu `$ to say

$`\overline{\kappa }_{\alpha \beta \gamma }`$ $`=`$ $`\overline{\kappa }(n_\alpha ,\mu _\alpha ,n_\beta ,\mu _\beta ,n_\gamma ,\mu _\gamma )`$
$`\overline{\kappa }(an_\alpha ,b\mu _\alpha ,an_\beta ,b\mu _\beta ,an_\gamma ,b\mu _\gamma )`$ $`=`$ $`a^ub^v\overline{\kappa }(n_\alpha ,\mu _\alpha ,n_\beta ,\mu _\beta ,n_\gamma ,\mu _\gamma ).`$ (A3)
In appendix B we show that $`u=1`$ and $`v=2`$. Since the mode frequencies and coupling coefficients are approximately independent of the $`m`$ quantum numbers, the $`m`$ dependence can be integrated over giving the function

$`A(n_\alpha ,n_\beta ,n_\gamma )`$ $`=`$ $`{\displaystyle _{n_\alpha }^{n_\alpha }}𝑑m_\alpha {\displaystyle _{n_\beta }^{n_\beta }}𝑑m_\beta {\displaystyle _{n_\gamma }^{n_\gamma }}𝑑m_\gamma \delta (m_\alpha +m_\beta +m_\gamma ).`$ (A4)
This function is symmetric, and is $`A4n_\alpha n_\beta `$ in either the limit $`n_\alpha n_\beta n_\gamma `$ or $`n_\beta n_\alpha n_\gamma `$.

Plugging these definitions into eq.17, and summing over the frequency signs we find

$`{\displaystyle 𝑑m_\alpha I_\alpha }=\pi (2/\pi )^2{\displaystyle 𝑑n_\beta n_\beta 𝑑n_\gamma n_\gamma 𝑑\mu _\beta 𝑑\mu _\gamma A|\overline{\kappa }_{\alpha \beta \gamma }|^2\mu _\alpha \mu _\beta \mu _\gamma \delta (\delta n)\delta (\delta k)N_\alpha N_\beta N_\gamma }`$ (A5)
$`\times `$ $`[({\displaystyle \frac{1}{N_\alpha }}{\displaystyle \frac{1}{N_\beta }}{\displaystyle \frac{1}{N_\gamma }})\delta (\mu _\alpha \mu _\beta \mu _\gamma )({\displaystyle \frac{1}{N_\beta }}{\displaystyle \frac{1}{N_\gamma }}{\displaystyle \frac{1}{N_\alpha }})\delta (\mu _\beta \mu _\gamma \mu _\alpha )`$
$``$ $`({\displaystyle \frac{1}{N_\gamma }}{\displaystyle \frac{1}{N_\alpha }}{\displaystyle \frac{1}{N_\beta }})\delta (\mu _\gamma \mu _\alpha \mu _\beta )].`$
We attempt to find scale-free inertial range solutions of the form

$`N`$ $`=`$ $`N_0n^p\mu ^q.`$ (A6)
After inserting this power law dependence on $`n`$ and $`\mu `$ into the eq.A5, one can make a change of coordinates in the second and third sets of terms in the last parenthesis, in order to make them have the same form as the first term, up to a scaling factor. For instance, in the second set of terms let $`n_\beta =n_\alpha ^2/n_\beta ^{}`$, $`\mu _\beta =\mu _\alpha ^2/\mu _\beta ^{}`$, $`n_\gamma =n_\gamma ^{}n_\alpha /n_\beta ^{}`$, and $`\mu _\gamma =\mu _\gamma ^{}\mu _\alpha /\mu _\beta ^{}`$. After simplifying and collecting terms we find

$`{\displaystyle 𝑑m_\alpha I_\alpha }=\pi N_0^2(2/\pi )^2{\displaystyle 𝑑n_\beta n_\beta 𝑑n_\gamma n_\gamma 𝑑\mu _\beta 𝑑\mu _\gamma A|\overline{\kappa }_{\alpha \beta \gamma }|^2\mu _\alpha \mu _\beta \mu _\gamma \delta (\delta n)\delta (\delta k)\delta (\mu _\alpha \mu _\beta \mu _\gamma )}`$ (A7)
$`\times `$ $`(n_\alpha n_\beta n_\gamma )^p(\mu _\alpha \mu _\beta \mu _\gamma )^q\left(n_\alpha ^p\mu _\alpha ^qn_\beta ^p\mu _\beta ^qn_\gamma ^p\mu _\gamma ^q\right)`$
$`\times `$ $`\left[1\left({\displaystyle \frac{n_\alpha }{n_\beta }}\right)^{2(p_0p)}\left({\displaystyle \frac{\mu _\alpha }{\mu _\beta }}\right)^{2(q_0q)1}\left({\displaystyle \frac{n_\alpha }{n_\gamma }}\right)^{2(p_0p)}\left({\displaystyle \frac{\mu _\alpha }{\mu _\gamma }}\right)^{2(q_0q)1}\right]`$
where $`p_0=u+3`$ and $`q_0=v+5/2`$. We can make the scaling of eq.A7 with $`n_\alpha `$ and $`\mu _\alpha `$ explicit by defining $`x_\beta =n_\beta /n_\alpha `$, $`x_\gamma =n_\gamma /n_\alpha `$, $`y_\beta =\mu _\beta /\mu _\alpha `$, $`y_\gamma =\mu _\gamma /\mu _\alpha `$, $`\kappa _{\alpha \beta \gamma }=n_\alpha ^u\mu _\alpha ^vf_{\alpha \beta \gamma }`$, and $`\overline{A}=A/n_\alpha ^2`$. We find

$`{\displaystyle 𝑑m_\alpha I_\alpha }`$ $`=`$ $`4N_0^2n_\alpha ^{2(p_0p)2}\mu _\alpha ^{2(q_0q)2}(p,q)`$ (A8)
where $`(p,q)`$ is the dimensionless integral

$`(p,q)`$ $`=`$ $`{\displaystyle 𝑑x_\beta x_\beta 𝑑x_\gamma x_\gamma 𝑑y_\beta 𝑑y_\gamma \overline{A}f_{\alpha \beta \gamma }^2y_\beta y_\gamma \delta (\delta x)\delta (\delta xy)\delta (\delta y)(x_\beta x_\gamma )^p(y_\beta y_\gamma )^q\left(1x_\beta ^py_\beta ^qx_\gamma ^py_\gamma ^q\right)}`$ (A9)
$`\times `$ $`\left(1x_\beta ^{2(pp_0)}y_\beta ^{2(qq_0)+1}x_\gamma ^{2(pp_0)}y_\gamma ^{2(qq_0)+1}\right).`$
Stationary solutions to eq. A7 are now easily found by using the delta functions to force either the first or second parenthesis to zero (Zakharov et al. 1992). For instance, using the frequency delta function to set the first parenthesis to zero gives $`N\mu ^1`$, so that all modes have the same energy, i.e., thermodynamic equilibrium. Using the momentum delta functions to set the first parenthesis to zero gives thermodynamic equilibria with respect to a moving reference frame. Setting the second parenthesis to zero using the momentum delta functions gives solutions supporting a constant momentum flux. We are interested in solutions which support an energy flux. The frequency delta function can be used to set the second parenthesis to zero if we let $`p=p_0=u+3=4`$ and $`q=q_0=v+5/2=1/2`$. We now proceed to show that this solution corresponds to a flux of energy to small frequency and large wavenumber.

In the inertial range, driving and damping are negligible, and the conserved energy flux has components $`^n`$ and $`^\mu `$ which satisfy 202020 The assumption of local energy transfer is later shown to be valid by verifying that the flux integrals converge.

$`\omega _\alpha {\displaystyle 𝑑m_\alpha I_\alpha }+\mathbf{}_k_\alpha `$ $`=`$ $`\omega _\alpha {\displaystyle 𝑑m_\alpha I_\alpha }+{\displaystyle \frac{1}{n}}{\displaystyle \frac{}{n}}\left(n_\alpha ^n\right)+{\displaystyle \frac{_\alpha ^\mu }{\mu }}=0.`$ (A10)
Care must be taken in evaluating eq.A8. For finite values of $`n_\alpha `$ and $`\mu _\alpha `$, taking the limit $`(p,q)(p_0,q_0)`$ gives $`𝑑m_\alpha I_\alpha =0`$. However, in the vicinity of $`n_\alpha =0`$ or $`\mu _\alpha =0`$ eq.A8 takes on the indeterminate form $`0/0`$, since $`0`$ in the numerator and either $`n_\alpha `$ or $`\mu _\alpha `$ goes to zero in the denominator. Following Zakharov et.al., we evaluate this expression using the delta function representation $`\text{lim}_{ϵ0}ϵ|x|^{ϵ1}=2\delta (x)`$ to find

$`\mu _\alpha n_\alpha {\displaystyle 𝑑m_\alpha I_\alpha }`$ $`=`$ $`4N_0^2\left({\displaystyle \frac{1}{\mu _\alpha }}{\displaystyle \frac{}{p}}(p_0,q_0)\delta (n_\alpha )+{\displaystyle \frac{1}{n_\alpha }}{\displaystyle \frac{}{q}}(p_0,q_0)\delta (\mu _\alpha )\right).`$ (A11)
Hence, the flux can only have a source for $`n=0`$ or $`\mu =0`$. Plugging eq.A11 into eq.A10, using $`\omega =2\mu `$, and integrating over $`n`$ or $`\mu `$, respectively, gives the two components of the flux

$`^n`$ $`=`$ $`8N_0^2n^1\mu ^1{\displaystyle \frac{}{p}}(p_0,q_0)`$
$`^\mu `$ $`=`$ $`8N_0^2n^2{\displaystyle \frac{}{q}}(p_0,q_0).`$ (A12)
We now estimate the dimensionless flux integrals, verifying that they converge and give have the correct sign. We do this by breaking the integration up into two regimes: large daughter mode wavenumber ($`x_\beta 1`$) and small daughter mode wavenumber ($`x_\beta 1`$). Although our expansions for the integrand are technically only valid in the respective limits, we extend them all the way to $`x_\beta 1`$.

When the derivatives with respect to $`p`$ and $`q`$ are taken in eq.A9, only the last parenthesis need be differentiated since it gives zero for $`(p,q)=(p_0,q_0)`$. The coupling coefficients for the two limits are given in appendix B. The resulting expressions are

$`{\displaystyle \frac{}{p}}(p_0,q_0)`$ $``$ $`2^{3/2}{\displaystyle _1^{\mathrm{}}}𝑑x_\beta x_\beta ^2\mathrm{ln}x_\beta +{\displaystyle \frac{128}{\pi ^3}}{\displaystyle _0^1}𝑑x_\beta x_\beta ^{1/2}\left(\mathrm{ln}x_\beta +1\right)`$ (A13)
$`=`$ $`2^{3/2}+{\displaystyle \frac{256}{9\pi ^3}}3.7\alpha _n`$
and

$`{\displaystyle \frac{}{q}}(p_0,q_0)`$ $`=`$ $`2^{3/2}\mathrm{ln}2{\displaystyle _1^{\mathrm{}}}𝑑x_\beta x_\beta ^2+{\displaystyle \frac{128}{\pi ^3}}{\displaystyle _0^1}𝑑x_\beta x_\beta ^{1/2}\mathrm{ln}x_\beta `$ (A14)
$`=`$ $`2^{3/2}\mathrm{ln}2{\displaystyle \frac{512}{9\pi ^3}}=3.8\alpha _\mu .`$
We note that the contribution to $`^n`$($`^\mu `$) from both large and small $`x_\beta `$ are positive (negative).

The energy flux is toward larger $`n`$ and smaller $`\mu `$. We now restore the sign of $`\mu `$ in order to have an expression valid for either sign. The final result for the energy flux is then

$`^n`$ $`=`$ $`8\alpha _nN_0^2n^1|\mu |^1`$
$`^\mu `$ $`=`$ $`8\alpha _\mu N_0^2n^2{\displaystyle \frac{\mu }{|\mu |}}.`$ (A15)
We also note the final answer for the occupation number

$`N`$ $`=`$ $`N_0n^4\mu ^{1/2}`$ (A16)
where $`N_0`$ is related to the energy flux by eq.A15.

## Appendix B  analytic estimate of the maximum coupling coefficient

Here we give a detailed analytic calculation of the maximum coupling coefficient. We make the following approximations: (1) $`nm`$, the WKB limit; (2) $`k_\varphi k_R`$, which follows from (1); (3) short wavelength daughter modes with $`\mu _\beta \mu _\gamma \mu _\alpha /2`$; (4) an $`n=1`$ polytropic background star. For convenience, we use units with $`r_{}=M=1`$, and we use the normalization condition $`E_{\mathrm{unit}}=0.5Mr_{}^2\mathrm{\Omega }^2`$. In these units, the density profile near the surface is $`\rho =\rho _0\widehat{z}`$, where the central density is $`\rho _0=\pi /4`$.

We decompose $`\psi `$ in eq.46 as a sum of plane waves as

$`\psi (𝒙)`$ $`=`$ $`{\displaystyle \frac{\psi _0}{(2\pi )^2}}\left({\displaystyle \frac{\rho }{\rho _0}}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\right)^{1/2}{\displaystyle \underset{s_1,s_2=\pm 1}{}}\mathrm{exp}(i\chi _{s_1s_2})`$ (B1)
where the WKB phase is

$`\chi _{s_1s_2}`$ $`=`$ $`p(s_1\theta _1+s_2\theta _2)+m\varphi +\alpha (s_1+s_2).`$ (B2)
In the small $`\mu `$ limit, we can write $`\theta _2\pi /2|\mu |ϵ`$ where $`1ϵ1`$. Hence the effective wavenumber in the $`\theta _2`$ direction is $`p|\mu |=\pi |k|`$. From eq.33 we find the displacement vector

$`𝝃(𝒙)`$ $`=`$ $`{\displaystyle \frac{\psi _0}{(2\pi )^2}}\left({\displaystyle \frac{\rho }{\rho _0}}\mathrm{sin}\theta _1\mathrm{sin}\theta _2\right)^{1/2}{\displaystyle \underset{s_1,s_2=\pm 1}{}}𝜺_{s_1s_2}\mathrm{exp}(i\chi _{s_1s_2})`$ (B3)
where

$`\epsilon _R`$ $`=`$ $`{\displaystyle \frac{i}{1q^2}}\left(k_R+iqk_\varphi \right)i\mu ^2k_R`$
$`\epsilon _\varphi `$ $`=`$ $`{\displaystyle \frac{i}{1q^2}}\left(k_\varphi iqk_R\right)\mu k_R\epsilon _R`$
$`\epsilon _z`$ $`=`$ $`ik_z\epsilon _\varphi .`$ (B4)
The wavenumber is defined by $`𝒌_{s_1s_2}=\chi _{s_1s_2}/𝒙`$ with components found from eq.38 to be

$`k_R`$ $`=`$ $`{\displaystyle \frac{p(1\mu ^2)^{1/2}}{x_1^2x_2^2}}\left(s_1x_1(1x_2^2)^{1/2}s_2x_2(1x_1^2)^{1/2}\right)`$
$`k_\varphi `$ $`=`$ $`{\displaystyle \frac{m}{R}}={\displaystyle \frac{m(1\mu ^2)^{1/2}}{[(1x_1^2)(1x_2^2)]^{1/2}}}`$
$`k_z`$ $`=`$ $`s_1s_2{\displaystyle \frac{|\mu |}{(1\mu ^2)^{1/2}}}k_R.`$ (B5)
The factor $`(x_1^2x_2^2)/|\mu |`$ is the volume element for the bi-spheroidal coordinates, and becomes zero at coordinate singularities. In the limit $`\widehat{z}1`$ and $`\mu 1`$, it has the simple form

$`(x_1^2x_2^2)^2`$ $``$ $`(y^2y_0^2)^2+\mathrm{\Delta }^2`$ (B6)
where $`y=\mathrm{cos}\theta `$, $`y_0(\mu ^22\widehat{z})^{1/2}`$, and $`\mathrm{\Delta }=(8\mu ^2\widehat{z})^{1/2}`$. This formula shows that there is a narrow peak for the wavenumber (for $`\widehat{z}\mu ^2/2`$) near the singular point $`(r,\mathrm{cos}\theta )=(1,\pm |\mu |)`$; there is no pronounced peak for $`\widehat{z}>\mu ^2/2`$ and the integrand decreases strongly.

Plugging the WKB travelling wave forms into eq.53 and keeping only the largest terms in the WKB limit gives

$`\kappa _{\alpha \beta \gamma }`$ $`=`$ $`{\displaystyle \frac{4\mathrm{\Omega }^2\rho _0}{2E_{\mathrm{unit}}}}{\displaystyle \frac{\psi _{0\alpha }\psi _{0\beta }\psi _{0\gamma }}{(2\pi )^6}}{\displaystyle \underset{\stackrel{}{s}_\alpha ,\stackrel{}{s}_\beta ,\stackrel{}{s}_\gamma }{}}{\displaystyle d^3x\frac{\rho }{\rho _0}\left(\mathrm{sin}\theta _{\alpha 1}\mathrm{sin}\theta _{\alpha 2}\mathrm{sin}\theta _{\beta 1}\mathrm{sin}\theta _{\beta 2}\mathrm{sin}\theta _{\gamma 1}\mathrm{sin}\theta _{\gamma 2}\right)^{1/2}}`$ (B7)
$`\times `$ $`\left({\displaystyle \frac{\rho _0^3}{\rho _\alpha \rho _\beta \rho _\gamma }}\right)^{1/2}M_{\alpha \beta \gamma }\mathrm{exp}[i(\chi _\alpha +\chi _\beta +\chi _\gamma )]`$
where the matrix element is defined by

$`M_{\alpha \beta \gamma }`$ $`=`$ $`\mu _\alpha ^2𝒌_\alpha 𝜺_\beta 𝒌_\alpha 𝜺_\gamma +\mu _\beta ^2𝒌_\beta 𝜺_\gamma 𝒌_\beta 𝜺_\alpha +\mu _\gamma ^2𝒌_\gamma 𝜺_\alpha 𝒌_\gamma 𝜺_\beta `$ (B8)
and $`\rho _\alpha `$ is the cutoff version of the WKB envelope.

To evaluate this expression, we first note that the integral is rapidly oscillating unless the conservation rules of eq. 54 are satisfied. These relations are written

$`s_{\alpha 1}p_\alpha +s_{\beta 1}p_\beta +s_{\gamma 1}p_\gamma `$ $`=`$ $`0`$
$`s_{\alpha 2}|k_\alpha |+s_{\beta 2}|k_\beta |+s_{\gamma 2}|k_\gamma |`$ $`=`$ $`0`$
$`m_\alpha +m_\beta +m_\gamma `$ $`=`$ $`0.`$ (B9)
In the limit of large daughter mode wavenumber, we find $`p_\beta p_\gamma `$, $`s_{\gamma 1}=s_{\beta 1}`$, $`|k_\gamma ||k_\beta |`$, $`s_{\gamma 2}=s_{\beta 2}`$, and $`m_\gamma m_\beta `$; in other words, the wavevectors of the daughter modes are equal in magnitude and opposite in direction. In this limit, the remaining spatially constant phase factor is

$`\mathrm{exp}[i(\chi _\alpha +\chi _\beta +\chi _\gamma )]`$ $`=`$ $`i^{\left[\delta _\alpha (s_{\alpha 1}+s_{\alpha 2})+(\delta _\beta \delta _\gamma )(s_{\beta 1}+s_{\beta 2})\right]}`$ (B10)
where $`\delta =0`$ for an even parity mode and $`\delta =1`$ for an odd parity mode. Using the incompressibility condition and momentum conservation leads to the simplification $`𝒌_\beta 𝜺_\gamma =𝒌_\alpha 𝜺_\gamma 𝒌_\alpha 𝜺_\beta `$. Plugging these relations in gives

$`M_{\alpha \beta \gamma }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mu _\alpha ^2𝒌_\alpha 𝜺_\beta \left(2𝒌_\alpha 𝜺_\beta +𝒌_\beta 𝜺_\alpha \right).`$ (B11)
In the limit of small $`\mu `$ and $`k_\varphi k_R`$ we find

$`𝒌_\alpha 𝜺_\beta `$ $``$ $`i{\displaystyle \frac{\mu _\alpha ^2}{4}}k_{\alpha R}k_{\beta R}+ik_{\alpha z}k_{\beta z}`$
$`𝒌_\beta 𝜺_\alpha `$ $``$ $`i\mu _\alpha ^2k_{\alpha R}k_{\beta R}+ik_{\alpha z}k_{\beta z}.`$ (B12)
Performing the spin sums, we find zero for net odd parity, while for even parity we get

$`{\displaystyle \underset{\stackrel{}{s}_\alpha ,\stackrel{}{s}_\beta ,\stackrel{}{s}_\gamma }{}}M_{\alpha \beta \gamma }i^{\left[\delta _\alpha (s_{\alpha 1}+s_{\alpha 2})+(\delta _\beta \delta _\gamma )(s_{\beta 1}+s_{\beta 2})\right]}={\displaystyle \frac{\mu _\alpha ^6p_\alpha ^2p_\beta ^2}{(x_{\alpha 1}^2x_{\alpha 2}^2)^2(x_{\beta 1}^2x_{\beta 2}^2)^2}}`$ (B13)
$`\times `$ $`(\{x_{\alpha 1}^2(1x_{\alpha 2}^2)+x_{\alpha 2}^2(1x_{\alpha 1}^2)\}\{x_{\beta 1}^2(1x_{\beta 2}^2)+x_{\beta 2}^2(1x_{\beta 1}^2)\}`$
$`+`$ $`4x_{\alpha 1}x_{\alpha 2}x_{\beta 1}x_{\beta 2}\left\{(1x_{\alpha 1}^2)(1x_{\alpha 2}^2)(1x_{\beta 1}^2)(1x_{\beta 2}^2)\right\}^{1/2}).`$
So far, we have

$`\kappa _{\alpha \beta \gamma }`$ $`=`$ $`{\displaystyle \frac{4\mathrm{\Omega }^2\rho _0}{2E_{\mathrm{unit}}}}{\displaystyle \frac{\psi _{0\alpha }\psi _{0\beta }^2}{(2\pi )^6}}\mu _\alpha ^6p_\alpha ^2p_\beta ^2𝒥={\displaystyle \frac{p_\alpha 𝒥}{4\pi ^{7/2}}}`$ (B14)
where the dimensionless integral is

$`𝒥`$ $`=`$ $`{\displaystyle d^3x\frac{\rho }{\rho _0}\left(\frac{\rho _0^3}{\rho _\alpha \rho _\beta ^2}\right)^{1/2}\left(\mathrm{sin}\theta _{\alpha 1}\mathrm{sin}\theta _{\alpha 2}\mathrm{sin}^2\theta _{\beta 1}\mathrm{sin}^2\theta _{\beta 2}\right)^{1/2}\frac{1}{(x_{\alpha 1}^2x_{\alpha 2}^2)^2(x_{\beta 1}^2x_{\beta 2}^2)^2}}`$ (B15)
$`\times `$ $`(\{x_{\alpha 1}^2(1x_{\alpha 2}^2)+x_{\alpha 2}^2(1x_{\alpha 1}^2)\}\{x_{\beta 1}^2(1x_{\beta 2}^2)+x_{\beta 2}^2(1x_{\beta 1}^2)\}`$
$`+`$ $`4x_{\alpha 1}x_{\alpha 2}x_{\beta 1}x_{\beta 2}\left\{(1x_{\alpha 1}^2)(1x_{\alpha 2}^2)(1x_{\beta 1}^2)(1x_{\beta 2}^2)\right\}^{1/2}).`$
To evaluate the dimensionless integral $`𝒥`$, we first note that if one attempted to take the $`\mu 0`$ limit, as for the mode energy, one would find an integral

$`𝒥`$ $``$ $`{\displaystyle _0^1}{\displaystyle \frac{dR}{R^{1/2}(1R^2)^2}}`$ (B16)
which implies a linear divergence at $`R=1`$. The divergence implies a strong dependence of the integral on some characteristic lengthscale near the surface.

To proceed with a more detailed calculation, first notice that the integrand is sharply peaked near the surface. Below the daughter mode turning point $`\widehat{z}_\beta `$, a factor of density cancels out. The parent mode WKB amplitude is largest for above the turning point $`\widehat{z}\widehat{z}_\alpha `$, where the cutoff density for the parent mode $`\rho _\alpha /\rho _0\widehat{z}_\alpha `$. Near the daughter mode singular point $`\mathrm{cos}\theta =y=|\mu _\beta |`$, the bi-spheroidal coordinates become $`x_{\alpha 1}|\mu _\alpha |`$, $`x_{\alpha 2}|\mu _\beta |`$, $`x_{\beta 1}x_{\beta 2}|\mu _\beta |`$. A similar expression can be found near the parent mode singular point. The integral $`𝒥`$ then becomes

$`𝒥`$ $`=`$ $`16\pi \widehat{z}_\alpha ^{1/2}{\displaystyle _{\widehat{z}_\beta }^{\widehat{z}_\alpha }}𝑑\widehat{z}{\displaystyle _1^1}{\displaystyle \frac{dy}{(y^2y_{\alpha 0}^2)^2+\mathrm{\Delta }_\alpha ^2}}+4\pi \widehat{z}_\alpha ^{1/2}{\displaystyle _{\widehat{z}_\beta }^{\widehat{z}_\alpha }}𝑑\widehat{z}{\displaystyle _1^1}{\displaystyle \frac{dy}{(y^2y_{\beta 0}^2)^2+\mathrm{\Delta }_\beta ^2}}`$ (B17)
$`=`$ $`16\pi \widehat{z}_\alpha ^{1/2}{\displaystyle _{\widehat{z}_\beta }^{\widehat{z}_\alpha }}𝑑\widehat{z}{\displaystyle \frac{\pi }{\mathrm{\Delta }_\alpha y_{\alpha 0}}}+4\pi \widehat{z}_\alpha ^{1/2}{\displaystyle _{\widehat{z}_\beta }^{\widehat{z}_\alpha }}𝑑\widehat{z}{\displaystyle \frac{\pi }{\mathrm{\Delta }_\beta y_{\beta 0}}}`$
$`=`$ $`{\displaystyle \frac{4\sqrt{2}\pi ^2}{\widehat{z}_\alpha ^{1/2}\mu _\alpha ^2}}{\displaystyle _{\widehat{z}_\beta }^{\widehat{z}_\alpha }}{\displaystyle \frac{\widehat{z}}{\widehat{z}^{1/2}}}+{\displaystyle \frac{4\sqrt{2}\pi ^2}{\widehat{z}_\alpha ^{1/2}\mu _\alpha ^2}}{\displaystyle _{\widehat{z}_\beta }^{\widehat{z}_\alpha }}{\displaystyle \frac{\widehat{z}}{\widehat{z}^{1/2}}}{\displaystyle \frac{16\sqrt{2}\pi ^2}{\mu _\alpha ^2}}.`$
The final result for the coupling coefficient for nearly identical, short lengthscale daughter modes is then

$`\kappa _{\alpha \beta \gamma }`$ $`=`$ $`={\displaystyle \frac{p_\alpha }{4\pi ^{7/2}}}{\displaystyle \frac{16\sqrt{2}\pi ^2}{\mu _\alpha ^2}}=1.02n_1\mu _1^2.`$ (B18)
In appendix A when finding the cascade solution, one also needs the coupling coefficient in the limit that one of the daughter modes has a small wavenumber. In section 4.1 we found that the approximate solutions in the limit $`n_\gamma n_\alpha n_\gamma `$ are: $`n_\gamma n_\alpha +n_\beta `$, $`k_\gamma k_\alpha `$, $`k_\beta k_\alpha (n_\beta /n_\alpha )^2`$, $`\mu _\gamma \mu _\alpha `$, and $`\mu _\beta \mu _\alpha (n_\beta /n_\alpha )\mu _\alpha `$. Since the method is exactly the same as the detailed calculation already given, we merely quote the answer for the maximum coupling coefficient to be

$`\kappa _{\alpha \beta \gamma }`$ $`=`$ $`\left({\displaystyle \frac{8}{\pi ^3}}\right)^{1/2}{\displaystyle \frac{p_\alpha }{|\mu _\alpha \mu _\beta |}}\left(\delta _O+\delta _S\mu _\beta /2\mu _\alpha \right).`$ (B19)
Here $`\delta _O=1`$ if modes $`\alpha `$ and $`\gamma `$ have the opposite parity, and $`\delta _O=0`$ otherwise. Similarly, $`\delta _S=1`$ if modes $`\alpha `$ and $`\gamma `$ have the same parity, and $`\delta _S=0`$ otherwise.