text:
# The rotational modes of relativistic stars: Numerical results

## I Introduction

That gravitational waves can drive various modes of oscillation in a rotating neutron star unstable was first suggested by Chandrasekhar chandra . The detailed mechanism behind this instability was explained by Friedman and Schutz fs78a ; fs78b , who showed that the instability sets in when an originally retrograde mode becomes prograde (according to an inertial observer) due to the rotation of the star. They also showed that the instability is generic, that, for arbitrarily slow rotation, any perfect-fluid stellar model has unstable modes with sufficiently large values of the azimuthal eigenvalue $`m`$ (the mode depends on $`\phi `$ as $`\mathrm{exp}(im\phi )`$). The existence of this radiation-driven instability is potentially important since it could limit the attainable spin rate of astrophysical neutron stars jlf83 . However, a detailed assessment of its astrophysical relevance is complicated by the fact that viscosity tends to counteract the growth of an unstable mode. One must account not only for the familiar hydrodynamic shear and bulk viscosities il , but also exotic mechanisms like the mutual friction that is relevant in superfluid neutron stars. Once a star has cooled below the superfluid transition temperature, mutual friction appears to suppress the instability of f-modes lm95 , and gravitational waves from f-modes appeared to set a limit on rotation only slightly more stringent that the maximum spin of an equilibrium model.

Recent work has modified this conclusion in two ways. First, numerical studies of the marginally stable “neutral” f-modes of fully relativistic, rapidly rotating polytropes by Stergioulas and Friedman sf  showed that relativistic effects destabilize the modes of a rotating star considerably. Most importantly, one finds that in general relativity the $`m=2`$ f-mode may have a neutral point for attainable rates of rotation. This is in contrast with the Newtonian result that the $`m=2`$ mode is unlikely to become unstable in uniformly rotating stars and it is important since the quadrupole mode is the most efficient emitter of gravitational radiation. Hence one would expect it to lead to the fastest growing instability. For realistic equations of state it has been shown msb98  that in a typical $`1.4M_{}`$ star the $`m=2`$ f-mode has a neutral point near $`\mathrm{\Omega }0.85\mathrm{\Omega }_K`$ where $`\mathrm{\Omega }_K`$ represents the mass-shedding limit.

The second piece of evidence follows from the fact that the inertial r-modes are unstable at any rate of rotation in a perfect fluid star na98 ; fm98 . The great surprise of a few years ago was the discovery that, despite the fact that they radiate mainly through the current multipoles, the unstable r-modes could potentially limit the spin of a rotating neutron star significantly lom98 ; aks99 . Since its original discovery the r-mode instability has been discussed in a number of papers, and we refer the interested reader to recent review articles for detailed discussions of the literature fl1 ; ak ; lindblom ; fredlamb ; fl2 . At the present time key issues concern the nonlinear evolution of an unstable mode sfont ; tohl1 ; tohl2 ; arras ; lin , the possible presence of hyperons in the neutron star core (since hyperons may lead to a very strong bulk viscosity) pbj ; lo02 ; hae  and the role of superfluidity lm00 ; ac01 .

The present paper concerns the effect that general relativity has on the instability of the r-modes and other inertial modes. Intuitively, one would expect this to be a relevant issue since the relativistic framedragging is an order $`\mathrm{\Omega }`$ effect which would affect inertial modes at leading order. In the last couple of years significant progress toward an understanding of the nature of inertial modes in general relativity has been made. The picture that is emerging is largely similar to that in Newtonian gravity. In the Newtonian case the r-modes are a purely axial parity subset of the large class of inertial modes which generally have velocity fields that are described by a hybrid mixture of axial and polar components to leading order in $`\mathrm{\Omega }`$. Non-barotropic Newtonian stars have an infinite set of r-modes for each permissible $`l`$ and $`m`$ (corresponding to the spherical harmonic $`Y_l^m`$ used to describe the velocity field), while barotropic models retain only a vestigial r-mode for each $`l=m`$, see lf ; yl  for detailed discussions. In relativity one can prove that barotropic stars (where the perturbations are described by the same one-parameter equation of state as the background star 111In Paper I such stars were refered to as “isentropic”.) have no purely axial inertial modes laf1 . All inertial modes of such stars are hybrids.

As in Newtonian theory ak , the relativistic non-barotropic r-mode problem is significantly different. Like their Newtonian counterparts, relativistic non-barotropic stars have modes that are purely axial in the spherical limit. These modes are determined by solving a single ordinary differential equation for one of the perturbed metric components koj , cf. eq (41). We have previously shown laf1  that discrete mode-solutions to this equation exist for uniform density stars. For more realistic equations of state, the problem is complicated by the fact that this equation corresponds to a singular eigenvalue problem which also admits a continuous spectrum koj ; bk . The dynamical role of this continuous spectrum, or indeed if it remains present when more physics is included in the model, is not clear at the present time. In order to determine a purely axial mode of a non-barotropic star one must typically identify a discrete mode embedded in the continuous spectrum. This technical difficulty has led to suggestions that the r-modes may not even “exist” for certain relativistic stars yosh ; kr01 . However, from eq. (41), it is clear that the slow-rotation approximation is no longer consistent in regions where $`\alpha \stackrel{~}{\omega }O(\mathrm{\Omega }^2)`$ or smaller. This means that the problem likely requires a “boundary layer” approach regul  where either $`\mathrm{\Omega }^2`$ terms, viscosity or the coupling to polar perturbations are included in the analysis of the region near the singular point. Recent results by Yoshida and Lee yl02  and Ruoff et al. rsk02  support this view.

So far there have been two studies of the growth timescale for the unstable relativistic r-modes rk02 ; yf01 . Both concern the non-barotropic problem. The results suggest that post-Newtonian estimates of the instability growth time are surprisingly good, but also indicate a weakening of the instability once the star becomes very compact. The aim of the present paper is to extend these results by considering the effect of radiation reaction on various inertial modes of barotropic stellar models. We expect the results obtained from a barotropic model to be relevant even though real neutron stars will have internal stratification associated with composition gradients (eg. due to the varying proton fraction). As discussed by Reisenegger and Goldreich rg  this will lead to a Brunt-Väisälä frequency $`N500\text{ s}^1(\rho /\rho _{\mathrm{nuclear}})^{1/2}`$ and the presence of core g-modes with frequencies in the range 100-200 Hz. In a very slowly rotating star the buoyancy force will dominate the Coriolis force and the g-modes will remain largely unchanged, but when $`\mathrm{\Omega }>>N`$ the situation will be reversed and the g-modes will be almost entirely rotationally restored. In other words, one would expect all low frequency modes to be well described by the inertial modes of a barotropic model for rotation periods shorter than 2-3 ms. This means that the results obtained in this paper should be a reasonable approximation for neutron stars in Low-Mass X-ray binaries (expected to have rotation rates in the range 250-500 Hz), and an accurate representation of the inertial modes of a newly born neutron star spinning at or near the break-up limit.

The layout of the paper is as follows: In section II we summarize the equations that describe the oscillations of a slowly rotating relativistic star. Section III provides a description of the numerical method we use to solve the eigenvalue problem, while the obtained results are discussed in Section IV. Section V is devoted to a discussion of our method for estimating the radiation reaction timescale for inertial modes. Our conclusions are given in Section VI.

Throughout the paper we use the following conventions: We refer to ref. lf  as Paper I, ref. laf1  as Paper II , while Paper III is regul . Our numbering convention for equations from these papers is such that e.g. (II,2.5) means Eqn. (2.5) from Paper II. Unless otherwise stated we use geometrized units $`c=G=1`$.

## II Perturbations of Slowly Rotating Stars

### II.1 The Equilibrium Star

We consider a perfect fluid star rotating slowly with uniform angular velocity $`\mathrm{\Omega }`$. By slow rotation we mean the assumption that $`\mathrm{\Omega }`$ is small compared to the Kepler velocity, $`\mathrm{\Omega }_K0.67\sqrt{\pi G\overline{ϵ}}`$, at which the star is unstable to mass shedding at its equator ($`\overline{ϵ}`$ represents the average energy density in the star). In particular, we neglect all quantities of order $`\mathrm{\Omega }^2`$ or higher. In this approximation the star retains its spherical shape, because the centrifugal deformation of its figure is an order $`\mathrm{\Omega }^2`$ effect h67 . The only new order $`\mathrm{\Omega }`$ effect that arises because of general relativity is the rotational framedragging, denoted by $`\omega (r)`$ below.

To first order in $`\mathrm{\Omega }`$ the equilibrium state is described by a stationary, axisymmetric spacetime with metric, $`g_{\alpha \beta }`$, of the form h67 ; cm74

$$ds^2=e^{2\nu (r)}dt^2+e^{2\lambda (r)}dr^2+r^2d\theta ^2+r^2sin^2\theta d\phi ^22\omega (r)r^2sin^2\theta dtd\phi $$
(1)
with perfect fluid matter source

$$T_{\alpha \beta }=(ϵ+p)u_\alpha u_\beta +pg_{\alpha \beta }.$$
(2)
Here, $`ϵ`$ and $`p`$ are, respectively, the total energy density and pressure of the fluid as measured by an observer moving with unit 4-velocity,

$$u^\alpha =e^\nu (t^\alpha +\mathrm{\Omega }\phi ^\alpha );$$
(3)
$`t^\alpha =(_t)^\alpha `$ and $`\phi ^\alpha =(_\phi )^\alpha `$ being, respectively, the timelike and rotational Killing vectors of the spacetime.

The metric and fluid variables are required to satisfy Einstein’s equation, $`G_{\alpha \beta }=8\pi T_{\alpha \beta }`$. This reduces to the well-known Tolman-Oppenheimer-Volkov (TOV) equations for the “spherical quantities” (those of order $`\mathrm{\Omega }^0`$) together with Hartle’s equation h67  for the order $`\mathrm{\Omega }`$ quantity,

$$\overline{\omega }(r)\mathrm{\Omega }\omega (r),$$
(4)
that governs the dragging of inertial frames induced by the rotation of the star. For the exact form of these equations, we refer the reader to Paper II; Eqs. (II,3.4-3.7) and (II,4.3).

To complete our specification of the equilibrium star we must provide an equation of state (EOS) relating the density and pressure. In this paper, we will always require our equilibrium solution to satisfy a one-parameter EOS, $`ϵ=ϵ(p)`$. This is an accurate assumption for equilibrium neutron stars since their temperature is likely to be significantly below the Fermi temperature $`T_F10^{12}`$ K. For simplicity we use the polytropic EOS,

$`p`$ $`=`$ $`K\rho ^{1+\frac{1}{n}}`$
$`ϵ`$ $`=`$ $`\rho +np,`$ (5)
where $`\rho `$ is the rest-mass density, $`n`$ is the polytropic index and $`K`$ is the polytropic constant. We use a set of polytropic indices ($`n=0.0,0.5,1.0,1.5`$) that span the range of compressibilities of proposed realistic neutron star equations of state glend ; akmal .

Although our equilibrium model obeys a one-parameter EOS, we do not necessarily require the perturbed fluid to satisfy the same EOS. For an adiabatic perturbation of an equilibrium star obeying a one-parameter EOS, the perturbed pressure and energy density are customarily related by

$$\frac{\delta p}{\mathrm{\Gamma }_1p}=\frac{\delta ϵ}{(ϵ+p)}+\xi ^\alpha A_\alpha $$
(6)
where $`\mathrm{\Gamma }_1(r)`$ is the adiabatic index, $`\xi ^\alpha `$ is the Lagrangian fluid displacement and where the Schwarzschild discriminant,

$$A_\alpha \frac{1}{(ϵ+p)}_\alpha ϵ\frac{1}{\mathrm{\Gamma }_1p}_\alpha p,$$
(7)
governs convective stability in the star. In general, the adiabatic index $`\mathrm{\Gamma }_1`$ need not be equal to the constant

$$\mathrm{\Gamma }\frac{(ϵ+p)}{p}\frac{dp}{dϵ}=1+\frac{1}{n}$$
(8)
associated with the equilibrium (polytropic) EOS. In terms of this constant we have

$$A_\alpha =\left(\frac{1}{\mathrm{\Gamma }}\frac{1}{\mathrm{\Gamma }_1}\right)\frac{1}{p}_\alpha p.$$
(9)
We will call a model barotropic if and only if the perturbed configuration satisfies the same one-parameter EOS as the unperturbed configuration. In this case $`\mathrm{\Gamma }_1\mathrm{\Gamma }`$ and the Schwarzschild discriminant vanishes identically. Such stars are marginally stable to convection, and since they have no internal stratification they do not admit finite frequency modes restored by buoyancy (g-modes) rg . In this paper we will consider perturbations of both barotropic and nonbarotropic stars, but will focus mainly on the barotropic case in our numerical work. The non-barotropic case has already been discussed in detail in laf1 ; yosh ; kr01 ; regul .

### II.2 The Perturbation Equations

We now consider the rotationally restored (inertial) modes of a slowly rotating relativistic star. The equations governing such perturbations were derived in detail in Paper II. Here we will simply quote the main results needed for the present analysis.

Since the equilibrium spacetime is stationary and axisymmetric, we may decompose our perturbations into oscillation modes proportional to $`e^{i(\sigma t+m\phi )}`$. For convenience, we will always choose $`m0`$, since the complex conjugate of an $`m<0`$ mode with real frequency $`\sigma `$ is an $`m>0`$ mode with frequency $`\sigma `$. Note that $`\sigma `$ is the mode frequency measured by an inertial observer at infinity.

In the Lagrangian perturbation formalism f78 ; fi92 , the basic variables are the metric perturbation $`h_{\alpha \beta }`$ and the Lagrangian displacement $`\xi ^\alpha `$. We begin by expanding these variables in vector and tensor spherical harmonics. The Lagrangian displacement vector can be written

$`\xi ^\alpha {\displaystyle \frac{1}{i\kappa \mathrm{\Omega }}}{\displaystyle \underset{l=m}{\overset{\mathrm{}}{}}}\left\{{\displaystyle \frac{1}{r}}W_l(r)Y_l^mr^\alpha +V_l(r)^\alpha Y_l^miU_l(r)P_\mu ^\alpha ϵ^{\mu \beta \gamma \delta }_\beta Y_l^m_\gamma t_\delta r\right\}e^{i\sigma t},`$ (10)
where we have defined,

$$P_\mu ^\alpha e^{(\nu +\lambda )}\left(\delta _\mu ^\alpha t_\mu ^\alpha t\right)$$
(11)
and introduced the “comoving” frequency,

$$\kappa \mathrm{\Omega }\sigma +m\mathrm{\Omega }.$$
(12)
The exact form of expression (10) has been chosen for convenience. In particular, we have used the gauge freedom inherent in the Lagrangian formalism cq76 ; ss77  to set $`\xi _t0`$. Note also the chosen relative phase between the terms in (10) with polar parity (those with coefficients $`W_l`$ and $`V_l`$) and the terms with axial parity (those with coefficients $`U_l`$).

Working in the Regge-Wheeler gauge rw57 , we express our metric perturbation as

$`h_{\mu \nu }=e^{i\sigma t}{\displaystyle \underset{l=m}{\overset{\mathrm{}}{}}}\left[\begin{array}{cccc}H_{0,l}(r)e^{2\nu }Y_l^m& H_{1,l}(r)Y_l^m& h_{0,l}(r)(\frac{m}{sin\theta })Y_l^m& ih_{0,l}(r)\mathrm{sin}\theta _\theta Y_l^m\\ H_{1,l}(r)Y_l^m& H_{2,l}(r)e^{2\lambda }Y_l^m& h_{1,l}(r)(\frac{m}{sin\theta })Y_l^m& ih_{1,l}(r)\mathrm{sin}\theta _\theta Y_l^m\\ \text{symm}& \text{symm}& r^2K_l(r)Y_l^m& 0\\ \text{symm}& \text{symm}& 0& r^2sin^2\theta K_l(r)Y_l^m\end{array}\right]`$ (17)
Again, note the choice of phase between the polar-parity components (those with coefficients $`H_{0,l}`$, $`H_{1,l}`$, $`H_{2,l}`$ and $`K_l`$) and the axial-parity components (those with coefficients $`h_{0,l}`$ and $`h_{1,l}`$).

In Paper II, we also found it convenient to make use of the Eulerian perturbation formalism, whose basic variables are the metric perturbation $`h_{\alpha \beta }`$ and the perturbed density, $`\delta ϵ`$, pressure, $`\delta p`$, and fluid 4-velocity, $`\delta u^\alpha `$. Since they represent scalar quantities the Eulerian changes in the density and pressure may be written as

$$\delta ϵ=\underset{l=m}{\overset{\mathrm{}}{}}\delta ϵ_l(r)Y_l^me^{i\sigma t}$$
(18)
and

$$\delta p=\underset{l=m}{\overset{\mathrm{}}{}}\delta p_l(r)Y_l^me^{i\sigma t},$$
(19)
respectively, while the Eulerian change in the fluid velocity may be expressed in terms of $`\xi ^\alpha `$ and $`h_{\alpha \beta }`$ as defined above,

$$\delta u^\alpha =q_\beta ^\alpha \mathrm{\pounds }_u\xi ^\beta +\frac{1}{2}u^\alpha u^\beta u^\gamma h_{\beta \gamma },$$
(20)
where $`q^{\alpha \beta }g^{\alpha \beta }+u^\alpha u^\beta `$.

We showed in Paper II that the rotationally restored modes of a slowly rotating star have a fundamentally different character depending on whether the star is barotropic or nonbarotropic. The difference pertains to the character of the modes in the limit as the star’s angular velocity, $`\mathrm{\Omega }`$, goes to zero. This is yet another reason why it is appropriate to consider the problem within the slow rotation approximation. In particular, we proved in Paper II that a relativistic barotrope does not admit distinct classes of r-modes or g-modes (modes whose limit as $`\mathrm{\Omega }0`$ are purely axial or purely polar, respectively). Instead, the generic inertial mode of such a star is a hybrid mixture of axial and polar components to lowest order in $`\mathrm{\Omega }`$. In contrast, non-barotropic stars allow distinct g-modes already in the non-rotating case and purely axial r-modes may exist at lowest order in $`\mathrm{\Omega }`$. To reflect the fundamental difference between the two cases, we organized our slow-rotation expansion by requiring our perturbation variables to obey the following ordering in powers of $`\mathrm{\Omega }`$,

$$U_l,h_{0,l}O(1)$$

$$W_l,V_l,H_{1,l}\{\begin{array}{cc}O(1)\hfill & \text{barotropic stars}\hfill \\ O(\mathrm{\Omega }^2)\hfill & \text{nonbarotropic stars}\hfill \end{array}$$

$$H_{0,l},H_{2,l},K_l,h_{1,l},\delta ϵ_l,\delta p_l,\sigma O(\mathrm{\Omega }).$$
(21)
The perturbation equations may then be grouped in powers of $`\mathrm{\Omega }`$ and solved order by order. To compare our results with the Newtonian r-modes and hybrid modes we need only find the leading order mode solutions; that is, we need only find the mode frequency to order $`\mathrm{\Omega }`$ and the eigenfunctions to $`O(1)`$. As discussed in Paper II, these turn out to be determined by a subset of the perturbation equations up to first order in $`\mathrm{\Omega }`$.

It is relevant to point out that the equations derived by this procedure will be somewhat different in different regions of spacetime (see Table 1). In the “near zone,” the region in which $`\sigma r<<1`$, we will be able to ignore second time derivatives, whereas we cannot do this in the “wave zone,” the region in which $`\sigma r>>1`$. Because the inertial modes that we are primarily interested in are restored by the Coriolis force, their frequencies scale with the angular velocity of the star, $`\sigma \mathrm{\Omega }`$. For slow rotation, this implies that the near zone extends far away from the star into the nonrelativistic region ($`M/r<<1`$) and that the wave zone will be located entirely within the nonrelativistic region ip71 . This will be important in Sect. IV when we calculate the energy radiated in gravitational waves and the timescales on which gravitational radiation reaction drives the unstable modes. For now, we will focus on the equations that are relevant in the near zone, which were derived in Paper II. These will allow us to find the eigenvalues and eigenfunctions of the modes that we are interested in. In Sect. IV.2 we will consider the equations more generally when we derive an expression for the radiated energy. The various spacetime regions are depicted schematically in Table 1.

The ordering (21) for the barotropic case is slightly more general than for the nonbarotropic case because the polar-parity coefficients $`W_l`$, $`V_l`$ and $`H_{1,l}`$ are not assumed to be negligible compared to the axial-parity coefficients $`U_l`$ and $`h_{0,l}`$. We will retain all of these variables in presenting the relevant equations and then specialize to each of the two cases. (Because $`h_{1,l}`$ is an order $`\mathrm{\Omega }`$ variable, we ignore it and drop the “0” subscript on $`h_{0,l}`$, writing it as $`h_l`$. Only in Sect. IV.2 will it be necessary to restore this distinction.)

As presented in Paper II, the relevant $`O(1)`$ equations (those that apply to the perturbed spherical star) are:

Eq. (II,3.20),

$$H_{1,l}+\frac{16\pi (ϵ+p)}{l(l+1)}e^{2\lambda }rW_l=0.$$
(22)
Eq. (II,3.23),

$$V_l=\frac{e^{(\nu +\lambda )}}{l(l+1)(ϵ+p)}\left[(ϵ+p)e^{\nu +\lambda }rW_l\right]^{}.$$
(23)
and Eq. (II,3.22),

$`h_l^{^{\prime \prime }}(\nu ^{}+\lambda ^{})h_l^{}+\left[{\displaystyle \frac{(2l^2l)}{r^2}}e^{2\lambda }{\displaystyle \frac{2}{r}}(\nu ^{}+\lambda ^{}){\displaystyle \frac{2}{r^2}}\right]h_l={\displaystyle \frac{4}{r}}(\nu ^{}+\lambda ^{})U_l`$ (24)
where a prime denotes a derivative with respect to $`r`$. Throughout this work we have used Eq. (22) \[i.e., Eq. (II,3.20)\] to eliminate the metric variable $`H_{1,l}(r)`$ in favor of $`W_l(r)`$.

To close the system of equations we need to retain only two of the equations that arise at $`O(\mathrm{\Omega })`$. As discussed in Paper II, the relevant pair are the two independent components of Eq. (II,4.16), which enforces the conservation of vorticity in constant entropy surfaces \[see Eq. (40) below\]. The $`[\theta \phi ]`$ component of Eq. (II,4.16) leads to Eq. (II,4.53),

$$0=\begin{array}{c}\left[l(l+1)\kappa \mathrm{\Omega }(h_l+U_l)2m\overline{\omega }U_l\right]\hfill \\ \\ +(l+1)Q_l\left[\frac{e^{2\nu }}{r}_r\left(r^2\overline{\omega }e^{2\nu }\right)W_{l1}2(l1)\overline{\omega }V_{l1}\right]\hfill \\ \\ lQ_{l+1}\left[\frac{e^{2\nu }}{r}_r\left(r^2\overline{\omega }e^{2\nu }\right)W_{l+1}+2(l+2)\overline{\omega }V_{l+1}\right],\hfill \end{array}$$
(25)
while the $`[r\theta ]`$ component of Eq. (II,4.16) gives Eq. (II,4.54),

$`0`$ $`=`$ $`{\displaystyle \frac{A_r}{(ϵ+p)}}\left[(l1)Q_l\left(i\delta p_{l1}+{\displaystyle \frac{_rp}{\kappa \mathrm{\Omega }r}}W_{l1}\right)(l+2)Q_{l+1}\left(i\delta p_{l+1}+{\displaystyle \frac{_rp}{\kappa \mathrm{\Omega }r}}W_{l+1}\right)\right]`$ (38)
$`+(l2)Q_{l1}Q_l\left[2_r\left(\overline{\omega }e^{2\nu }U_{l2}\right)+{\displaystyle \frac{(l1)}{r^2}}_r\left(r^2\overline{\omega }e^{2\nu }\right)U_{l2}\right]`$
$`+Q_l[\begin{array}{c}(l1)\kappa \mathrm{\Omega }_r\left(e^{2\nu }V_{l1}\right)2m_r\left(\overline{\omega }e^{2\nu }V_{l1}\right)\hfill \\ \\ +\frac{m(l1)}{r^2}_r\left(r^2\overline{\omega }e^{2\nu }\right)V_{l1}+(l1)\kappa \mathrm{\Omega }e^{2\nu }(\frac{16\pi r(ϵ+p)}{(l1)l}\frac{1}{r})e^{2\lambda }W_{l1}]\hfill \end{array}`$
$`+[\begin{array}{c}m\kappa \mathrm{\Omega }_r\left[e^{2\nu }(h_l+U_l)\right]+2_r\left(\overline{\omega }e^{2\nu }U_l\right)\left((l+1)Q_l^2lQ_{l+1}^2\right)\hfill \\ \\ +\frac{1}{r^2}_r\left(r^2\overline{\omega }e^{2\nu }\right)U_l[m^2+l(l+1)(Q_{l+1}^2+Q_l^21)]]\hfill \end{array}`$
$`Q_{l+1}[\begin{array}{c}(l+2)\kappa \mathrm{\Omega }_r\left(e^{2\nu }V_{l+1}\right)+2m_r\left(\overline{\omega }e^{2\nu }V_{l+1}\right)\hfill \\ \\ +\frac{m(l+2)}{r^2}_r\left(r^2\overline{\omega }e^{2\nu }\right)V_{l+1}+(l+2)\kappa \mathrm{\Omega }e^{2\nu }(\frac{16\pi r(ϵ+p)}{(l+1)(l+2)}\frac{1}{r})e^{2\lambda }W_{l+1}]\hfill \end{array}`$
$`+(l+3)Q_{l+1}Q_{l+2}\left[2_r\left(\overline{\omega }e^{2\nu }U_{l+2}\right)+{\displaystyle \frac{(l+2)}{r^2}}_r\left(r^2\overline{\omega }e^{2\nu }\right)U_{l+2}\right].`$
The constants $`Q_l`$ were defined in Paper II to be

$$Q_l\left[\frac{(l+m)(lm)}{(2l1)(2l+1)}\right]^{1/2}.$$
(39)
Notice that in writing Eq. (38) we have retained (in the first line) the term containing the Schwarzschild discriminant $`A_r`$ from the right hand side of Eq. (II,4.16),

$$i\kappa \mathrm{\Omega }e^\nu \mathrm{\Delta }\omega _{\alpha \beta }=\frac{2}{n}A_r_{[\alpha }r_{\beta ]}\mathrm{\Delta }p.$$
(40)
We have done this in order to emphasize the difference between the barotropic and nonbarotropic cases. To find the rotational modes to leading order in $`\mathrm{\Omega }`$ we need a complete set of equations involving only our $`O(1)`$ variables. If $`A_r`$ is not identically zero, the retained term in (38) brings about a coupling between these variables and the $`O(\mathrm{\Omega })`$ variable $`\delta p_l`$. In a barotropic star, $`A_r`$ is identically zero, so this coupling between $`O(1)`$ and $`O(\mathrm{\Omega })`$ variables vanishes. In this case, the five equations (22)-(38) do, indeed, involve only the $`O(1)`$ variables $`h_l`$, $`U_l`$, $`W_l`$, $`V_l`$ and $`H_{1,l}`$. These equations, therefore, comprise a complete set and fully determine our normal mode eigenvalue problem. In a nonbarotropic star, however, $`A_r0`$, and the coupling between $`O(1)`$ and $`O(\mathrm{\Omega })`$ variables in Eq. (38) does not vanish. In this case, the equations (22)-(38) do not involve only the variables $`h_l`$, $`U_l`$, $`W_l`$, $`V_l`$ and $`H_{1,l}`$. There are then two options: The first is reminiscent of the situation for r-modes of non-barotropic Newtonian stars. The problem would become well posed if we extended the analysis to one order higher in $`\mathrm{\Omega }`$ ak  and thus determined also $`\delta p_l`$ etcetera. The second possibility would be to obtain a well-defined eigenvalue problem by assuming the non-barotropic ordering (21) in which only the axial variables $`U_l`$ and $`h_l`$ are $`O(1)`$. With this choice, one obtains an eigenvalue problem from Eqs. (24) and (25) \- dropping $`W_l`$ and $`V_l`$ from the latter as $`O(\mathrm{\Omega }^2)`$ quantities. These equations can be reexpressed as Eq. (II,4.24): Kojima’s koj  eigenvalue equation governing $`h_l`$,

$$(\alpha \stackrel{~}{\omega })\left\{e^{\nu \lambda }\frac{d}{dr}\left[e^{\nu \lambda }\frac{dh_l}{dr}\right]\left[\frac{l(l+1)}{r^2}\frac{4M}{r^3}+8\pi (ϵ+p)\right]h_l\right\}+16\pi (ϵ+p)\alpha h_l=0,$$
(41)
where $`\stackrel{~}{\omega }=\overline{\omega }/\mathrm{\Omega }`$ and $`\alpha =l(l+1)\kappa /2m`$ and where $`U_l`$ is determined (once the eigenvalue, $`\alpha `$, and $`h_l`$ are known) by Eq. (25), which becomes,

$$U_l=\frac{\alpha }{(\stackrel{~}{\omega }\alpha )}h_l.$$
(42)
To summarize: The eigenvalue problem governing the rotational modes of a barotropic relativistic star to lowest order in $`\mathrm{\Omega }`$ is determined by the set of equations (22)-(38) for the variables $`h_l`$, $`U_l`$, $`W_l`$, $`V_l`$ and $`H_{1,l}`$, which are all assumed to be zeroth order in $`\mathrm{\Omega }`$. By contrast, the eigenvalue problem governing the rotational modes of a non-barotropic star to lowest order in $`\mathrm{\Omega }`$ is determined by Eqs. (41) and (42) for the variables $`h_l`$ and $`U_l`$, which alone are assumed to be zeroth order in $`\mathrm{\Omega }`$.

To complete the specification of our eigenvalue problem we must impose appropriate boundary conditions. We require our variables to be regular everywhere in the spacetime, which implies that they vanish at the origin, $`r=0`$. The fluid variables $`U_l(r)`$ and $`V_l(r)`$ are otherwise unconstrained - apart from the fact that they vanish outside the star. In Paper II, we derived the boundary condition (II,4.64) on the fluid variable $`W_l(r)`$ at the surface of the star, $`r=R`$,

$$W_l(R)=0\text{(all }l\text{)}.$$
(43)
This condition and Eq. (22) imply that the metric variable $`H_{1,l}(r)`$ vanishes at the surface, and in the exterior, of the star. By inspection, it is clear that all but one of the perturbation equations vanish in the vacuum exterior to the star ($`r>R`$). The one surviving equation is that governing the metric variable $`h_l(r)`$, i.e., Eq. (24) for barotropic stars or Eq. (41) for nonbarotropic stars. In the near zone, where these equations apply, they both reduce to the same equation outside the star: Eq. (II,4.66), or,

$$(1\frac{2M}{r})\frac{d^2h_l}{dr^2}\left[\frac{l(l+1)}{r^2}\frac{4M}{r^3}\right]h_l=0,$$
(44)
where $`M`$ is the total gravitational mass of the star. This equation has a regular singular point at $`r=\mathrm{}`$, so it has at least one regular series expansion about this point. The other, linearly independent, solution turns out to be singular at $`r=\mathrm{}`$, growing like $`r^{l+1}`$. As we will see in Sect. IV.2 we may ignore this second solution in the near zone because it is of order $`(\sigma r)^{2l+1}`$ relative to the nonsingular solution. To find the mode eigenvalues and eigenfunctions to lowest order in $`\mathrm{\Omega }`$ we therefore need only the regular solution. As discussed in Paper II the regular solution can be immediately written down as Eq. (II,4.67), or

$$h_l(r)=\underset{s=0}{\overset{\mathrm{}}{}}\widehat{h}_{l,s}\left(\frac{R}{r}\right)^{l+s},$$
(45)
where the coefficients $`\widehat{h}_{l,s}`$ are given by the recursion relation (II,4.68),

$$\widehat{h}_{l,s}=\left(\frac{2M}{R}\right)\frac{(l+s2)(l+s+1)}{s(2l+s+1)}\widehat{h}_{l,s1}$$
(46)
or, equivalently, by,

$$\widehat{h}_{l,s}=\frac{(l+s2)!(l+s+1)!(2l+1)!}{s!(l2)!(l+1)!(2l+s+1)!}\left(\frac{2M}{R}\right)^s\widehat{h}_{l,0}$$
(47)
with $`\widehat{h}_{l,0}`$ an arbitrary normalization constant. This known exterior solution must be matched to the interior solution for $`h_l(r)`$ at the surface of the star. This provides a boundary condition on the interior solution. We require that the solutions be continuous at the surface,

$$\underset{\epsilon 0}{lim}\left[h_l(R\epsilon )h_l(R+\epsilon )\right]=0,$$
(48)
for all $`l`$ (which fixes the normalization constant $`\widehat{h}_{l,0}`$), and that the Wronskian of the interior and exterior solutions vanish at $`r=R`$, i.e. that

$$\underset{\epsilon 0}{lim}\left[h_l(R\epsilon )h_l^{}(R+\epsilon )h_l^{}(R\epsilon )h_l(R+\epsilon )\right]=0,$$
(49)
for all $`l`$.

Finally, we note that since we are working in linearized perturbation theory there is a scale invariance to the equations. If $`(\xi ^\alpha ,h_{\alpha \beta })`$ is a solution to the perturbation equations then $`(S\xi ^\alpha ,Sh_{\alpha \beta })`$ is also a solution, for constant $`S`$. We will sometimes find it convenient to impose the following normalization condition in addition to the boundary and matching conditions just discussed:

$$\{\begin{array}{cc}U_m(r=R)=1\hfill & \text{for axial-led inertial modes and r-modes,}\hfill \\ U_{m+1}(r=R)=1\hfill & \text{for polar-led inertial modes.}\hfill \end{array}$$
(50)
## III Mode Frequencies and Eigenfunctions

From the equations in the previous section it is clear that the task of determining the relativistic analogues of the Newtonian r-modes is different depending on whether the star is barotropic or non-barotropic. Each problem presents it’s own computational challenge. For non-barotropic stars one must deal with the fact that Kojima’s equation (41) represents a singular eigenvalue problem. The barotropic case is conceptually easier because one does not have to worry about singularities. On the other hand, all inertial modes of a barotropic star will have a hybrid nature which complicates the numerical solution of the eigenvalue problem considerably. In this paper we focus our attention on barotropic stars since: i) one would expect the inertial modes of a more complex stellar model (with internal stratification) to be similar to those of a barotropic star for sufficiently rapid rotation, ii) the nonbarotropic case has already been discussed by several authors yosh ; kr01 ; regul .

We have used the spectral method described in Appendix A to solve numerically for a sample of rotationally restored modes of slowly rotating, fully relativistic barotropes. The relevant set of equations for this problem is Eqs. (22)-(38), which comprise a system of ordinary differential equations for the variables $`U_l(r)`$, $`V_l(r)`$, $`W_l(r)`$ $`H_{1,l}(r)`$ and $`h_l(r)`$ (for all allowed $`l`$). Together with the boundary and matching conditions at the surface of the star, these equations form a nonlinear eigenvalue problem for the parameter $`\kappa `$, the dimensionless mode frequency in the rotating frame. \[In practice, because the metric variable $`H_{1,l}`$ is related algebraically to the fluid variable $`W_l`$ by Eq. (22), we eliminate $`H_{1,l}`$ from the other equations and solve the system (23)-(38).\]

For simplicity, we have restricted our study to relativistic polytropes even though it would be straightforward to generalise our calculation to tabulated realistic equations of state. This is an important step that should eventually be taken, but we feel that we should first try to understand the overall effect that general relativity has on the inertial modes of a compact star.

In Newtonian barotropic stars there remained a large set of modes that were purely axial to lowest order in $`\mathrm{\Omega }`$: the r-modes with spherical harmonic indices $`l=m`$. The $`l=m=2`$ r-mode is the one expected to dominate the gravitational wave-driven instability of sufficiently hot and rapidly rotating neutron stars lom98 ; aks99 . We showed in Paper II, however, that relativistic barotropes do not admit such modes. Pure r-modes with $`l=m2`$ are not allowed by the perturbation equations (24)-(38). The corresponding modes are instead axial-led hybrid modes. We have previously solved explicitly for these important hybrid modes to first post-Newtonian order in uniform density stars \[cf. Eqs. (II,5.33)-(II,5.40)\]. We now report on a more general numerical study of these and other hybrid inertial modes in fully relativistic stars.

Apart from the replacement of the Newtonian r-modes with hybrids, the structure of the inertial mode spectrum in relativistic stars appears to be identical to that in Newtonian stars. That is, there is a one-to-one correspondence between the rotationally restored modes of a Newtonian star and those of the corresponding relativistic model. One can make a more formal comparison by constructing sequences of relativistic models of increasing compactness, $`M/R`$, a suitable measure of the importance of relativistic effects. A relativistic star with small $`M/R`$ agrees well in all of its physical properties with the Newtonian model having the same EOS and central density. In each star along the relativistic sequence we search for inertial modes using the method described in Appendix A. We find that for each of the Newtonian modes considered, there exists a family of relativistic modes parametrized by $`M/R`$ that approaches the Newtonian mode as $`M/R0`$.

To test the accuracy of our numerical code we compared its results with our post-Newtonian solution from Paper II in the small $`M/R`$ regime (see Fig. 1). The solution (II,5.33)-(II,5.40) gives explicitly those axial-led inertial modes of a relativistic uniform density star that correspond to the $`l=m`$ Newtonian r-modes. For these modes the dimensionless comoving frequency in the Newtonian star is simply

$$\kappa _\text{N}=\frac{2}{(m+1)},$$
(51)
while the post-Newtonian calculation gives Eq. (II,5.33), or,

$$\kappa _{\text{1pN}}=\frac{2}{(m+1)}\left[1\frac{8(m1)(2m+11)}{5(2m+1)(2m+5)}\left(\frac{M}{R}\right)\right].$$
(52)
Fig. 1 shows that the numerically computed eigenvalues agree well with the post-Newtonian solution; the differences being of order $`(M/R)^2`$ as expected. The numerically computed eigenfunctions also agree well with the post-Newtonian eigenfunctions (II,5.34)-(II,5.40). This agreement with our analytic solution gives us confidence that our code is able to find the relativistic modes. Thus, we may now explore the fully relativistic regime and, in particular, consider modes for which we have not worked out a post-Newtonian solution. Fig. 1 also shows the numerically computed eigenvalues for highly relativistic uniform density models. We have used the Newtonian frequency $`\kappa _\mathrm{N}`$ to normalize the results. This makes it easy to see that all modes clearly approach their Newtonian values as $`M/R0`$. Furthermore, the results show that the mode-frequency tends to decrease as the star becomes more compact. It is anticipated that general relativity will have this effect laf1 . The gravitational redshift will tend to decrease the fluid oscillation frequencies measured by a distant inertial observer. Also, because these modes are rotationally restored they will be affected by the dragging of inertial frames induced by the star’s rotation. Specifically, since the Coriolis force is determined by the fluid’s angular velocity relative to that of the local inertial frame $`\overline{\omega }(r)=\mathrm{\Omega }\omega (r)`$ it decreases, and the modes oscillate less rapidly, as the dragging of inertial frames becomes more pronounced. This is, of course, not too surprising given 52.

Next we consider how the relativistic analogues to the Newtonian r-modes are affected by changes in the stiffness of the equation of state. We can easily do this by determining the modes for a varying polytropic index, $`n`$, at a given compactness. Fig. 2 provides such results for $`M/R=0.15`$ which would be a typical value for realistic neutron stars. From this data we see that the mode frequencies differ from the Newtonian values by $`1015\%`$, depending on the stiffness of the equation of state. In order to confirm that these results are typical, we have carried out calculations for a much larger set of modes (inluding also polar-led inertial modes). Results for all modes whose Newtonian frequencies $`\kappa _\text{N}`$ were listed in Tables 5 and 6 of Paper I, show shifts in the frequency that accord well with the results in Fig. 2.

We now turn to a discussion of the eigenfunctions. Because the Newtonian studies have found the $`l=m`$ r-modes to be the most unstable to gravitational radiation-reaction, we would like to understand the nature of these modes in relativistic stars. As we have already mentioned, these modes become axial-led hybrids in relativistic barotropic models. In terms of the spherical harmonic expansion of the fluid displacement, Eq. (10), the statement that the Newtonian modes are pure r-modes with $`l=m`$ means that only the coefficient $`U_m(r)`$ is nonzero (to lowest order in $`\mathrm{\Omega }`$). It has the form

$$U_m(r)=\left(\frac{r}{R}\right)^{m+1}.$$
(53)
\[We have normalized the mode according to Eq. (50) so that $`U_m(r)=1`$ at the surface of the star, $`r=R`$.\] In a relativistic barotrope, on the other hand, other coefficients in Eq. (10) will be nonzero to lowest order in $`\mathrm{\Omega }`$. In addition we can examine the nonzero functions $`h_l(r)`$ ($`h_{0,l}`$) and $`H_{1,l}(r)`$ in Eq. (17), the spherical harmonic expansion of the perturbed metric.

Figs. 3 and 4 show some of these nonvanishing coefficients for the axial-led inertial mode whose Newtonian counterpart is the $`l=m=2`$ r-mode. The mode is shown for two relativistic models: a uniform density star ($`n=0`$ polytrope) and an $`n=1`$ polytrope; both with compactness $`M/R=0.15`$. At this compactness, the mode of the $`n=0`$ model has eigenvalue $`\kappa =0.5991`$ and the mode of the $`n=1`$ model has eigenvalue $`\kappa =0.5907`$. \[In Newtonian gravity the eigenvalue is $`\kappa _\text{N}=2/(m+1)=2/3`$.\] Fig. 3 shows the fluid displacement functions $`U_l(r)`$, $`W_l(r)`$, and $`V_l(r)`$ for $`l4`$ as well as the Newtonian r-mode function $`U_2(r)=(r/R)^3`$ (dashed curve). Meanwhile, Fig. 4 shows the metric functions $`h_l(r)`$ and $`H_{1,l}(r)`$ for $`l4`$. The amplitudes of all of these functions are determined by the normalization condition, Eq. (50), and reveal that the relativistic corrections to the Newtonian r-mode eigenfunction are only of the order of a few percent even for these highly relativistic models. The fact that $`h_2(r)`$ dominates the perturbed metric reveals directly that this mode couples most strongly to current quadrupole radiation, cf. estimates based on the Newtonian mode results lom98 ; aks99 .

Given the existence of a large number of inertial modes of a rotating star, it is interesting to ask whether other such modes may be driven unstable by gravitational radiation. That this is, indeed, the case was shown in Paper I (see also yl ). A particular axial-led mode which is interesting since one would a priori expect it to couple strongly to current quadrupole radiation is shown in the series of figures 5-7. This mode is the relativistic counterpart of one of the Newtonian $`m=2`$ axial-led hybrids studied in Paper I. (In particular, it is the mode whose Newtonian frequency is $`\kappa _\text{N}=0.4669`$ in the uniform density model and $`\kappa _\text{N}=0.5173`$ in the $`n=1`$ polytrope, cf. Paper I.)

Fig. 5 shows the Newtonian mode in a uniform density model together with its relativistic counterpart at the nearly Newtonian compactness of $`M/R=10^4`$. The Newtonian and weakly relativistic fluid functions $`U_l(r)`$, $`W_l(r)`$, and $`V_l(r)`$ for $`l4`$ are indistinguishable - the relativistic corrections being of order $`10^5`$, as one would expect. This scale for the relativistic corrections is also indicated by the leading metric functions $`h_l(r)`$ and $`H_{1,l}(r)`$ and by the fact that the functions with $`l>4`$ are smaller than those shown in both panels by a factor of order $`10^5`$ or smaller.

As discussed in Paper I, it is possible to find exact, analytic expressions for the Newtonian inertial modes of a uniform density star. The explicit forms of the fluid functions shown in Fig. 5 are given in Table 3 of Paper I (but with a different normalization). The functions with $`l>4`$ are identically zero for the Newtonian uniform density model, so the eigenfunction shown is indeed an exact solution. (For the Newtonian $`n=1`$ model the fluid eigenfunctions are similar but not identical to those shown and have nonvanishing higher $`l`$ terms of order $`0.5\%`$ or smaller.)

The fact that $`U_2(r)`$ is of order unity relative to the other fluid functions suggests that this mode ought to lead to significant current quadrupole radiation. Since it satisfies the CFS instability criterion the mode will be driven unstable by gravitational radiation and one might even expect it to make as important a contribution to the gravitational-wave driven spin-down of a hot, rapidly rotating neutron star as the much-discussed $`l=m=2`$ r-mode. However, a more detailed calculation of the growth timescale of the unstable mode does not support these expectations. Instead, one finds that the current quadrupole radiation from this mode is negligible compared to higher $`l`$ multipole moments for the Newtonian models. In Paper I, it was shown that the current quadrupole associated with this mode actually vanishes identically for the uniform density model and that it is negligibly small for the $`n=1`$ polytrope. The relativistic calculation provides us with a way to understand this result since we can now examine the perturbed metric in the exterior spacetime for a weakly relativistic model. From the results in Fig. 5 we see directly that the metric function $`h_2(r)`$ essentially vanishes outside the star compared to the function $`h_4(r)`$, thus implying that the exterior perturbed metric (i.e. the emerging radiation) is dominated by the $`l=4`$ current multipole.

It is interesting to consider whether this result will still be true for strongly relativistic stars. In Figs. 6 and 7 we show the same mode in two stellar models with compactness $`M/R=0.15`$. For this compactness, the mode of the $`n=0`$ model (Fig. 6) has eigenvalue $`\kappa =0.3879`$ and the mode of the $`n=1`$ model (Fig. 7) has eigenvalue $`\kappa =0.4313`$. That is, the frequencies are 17% smaller than the Newtonian ones. For clarity, we display the coefficients of the axial and polar-parity terms in the spherical harmonic expansions (10) and (17) in separate panels. The upper left panel of each figure shows the axial-parity fluid functions $`U_l(r)`$ for $`l6`$, while the upper right panel shows the polar parity fluid functions $`W_l(r)`$, and $`V_l(r)`$ for $`l5`$. (Their Newtonian counterparts are also shown for comparison.) The lower left and right panels show, respectively, the axial metric functions $`h_l(r)`$ and the polar metric functions $`H_{1,l}(r)`$ for $`l6`$. These figures indicate that the relativistic corrections to this mode are only of the order of $`1\%`$ even for these strongly relativistic models. The results also show that the metric function $`h_2(r)`$ is not completely dominated by $`h_4(r)`$ in the exterior spacetime. For both stars, $`h_2`$ is of order $`0.1\%`$ at the surface of the star (where, as usual, the normalization is fixed relative to $`U_2(R)=1`$). This is comparable to $`h_4`$ for the $`n=1`$ star and only a factor of 10 smaller than $`h_4`$ for the $`n=0`$ star. This is an interesting result since it suggests that the Newtonian prediction that the radiation field is dominated by the $`l=4`$ current multipole is not borne out by the fully relativistic calculation. Indeed, in Sect. IV.3 we will see that the $`l=2`$ current multipole dominates over the $`l=4`$ multipole. We will also show that the growth timescale of the mode is nevertheless much longer than that of the $`l=m=2`$ r-mode.

Finally, we present a set of results for the fastest growing unstable polar-led inertial mode in Figs. 8-10. This mode is the relativistic counterpart of one of the Newtonian $`m=2`$ polar-led hybrids studied in Paper I. (In particular, it is the mode whose Newtonian frequency is $`\kappa _\text{N}=1.232`$ in the uniform density model and $`\kappa _\text{N}=1.100`$ in the $`n=1`$ polytrope.) Fig. 8 shows the Newtonian mode in a uniform density model together with its relativistic counterpart at the nearly Newtonian compactness of $`M/R=10^4`$. Again, the Newtonian and weakly relativistic fluid functions $`U_l(r)`$, $`W_l(r)`$, and $`V_l(r)`$ for $`l3`$ are indistinguishable — the relativistic corrections and leading metric functions $`h_l(r)`$ and $`H_{1,l}(r)`$ being of order $`10^4`$. As with the axial-led mode shown in Fig. 5, an exact expression for the mode eigenfunction in the Newtonian uniform density star is given in Paper I (Table 4) with the fluid functions vanishing identically for $`l>3`$. (Again, for the Newtonian $`n=1`$ model the fluid eigenfunctions are similar but not identical to those shown in Fig. 8 and have nonvanishing higher $`l`$ terms of order $`1\%`$ or smaller.)

In Figs. 9-10 we show the mode in two stellar models with compactness $`M/R=0.15`$. At this compactness, the mode of the $`n=0`$ model (Fig. 9) has eigenvalue $`\kappa =1.028`$ and the mode of the $`n=1`$ model (Fig. 10) has eigenvalue $`\kappa =0.9068`$. This means that the frequencies are again about 17% smaller than the Newtonian results. For clarity, we display separately the coefficients of the axial and polar-parity terms in the spherical harmonic expansions (10) and (17). These figures show that the relativistic corrections to this mode are, again, only of the order of $`1\%`$ even for these strongly relativistic models. We can also see that the metric function $`h_3(r)`$ dominates the perturbed metric exterior to the star, which implies that this mode couples most strongly to current octupole radiation.

## IV Radiation Reaction Timescales

Almost all previous estimates of the strength of the r-mode instability are based on Newtonian and post-Newtonian calculations. One of the main goals of this work is to determine whether general relativity will have a significant effect on the stability of the modes, and, if so, to determine whether it will make them more or less unstable. This is obviously a crucial issue. In particular since it is known that relativistic effects change the instability point for the f-modes considerably sf . This fact underlines why radiation driven instabilities need to be studied in the full framework of general relativity. As far as unstable inertial modes are concerned, there has so far been two studies of the associated growth times. Both concern the r-modes of non-barotropic relativistic stars. Yoshida and Futamase yf01  implemented the near-zone boundary conditions discussed by Lindblom, Mendell and Ipser nearzone , while Ruoff and Kokkotas rk02  solved the complex eigenvalue problem 222Correcting an earlier calculation by Andersson na98  that was inconsistent in that it included some terms of order $`\mathrm{\Omega }^2`$, while neglecting others. posed by the slow-rotation equations when a condition of outgoing radiation at infinity is imposed. These two studies provide useful insights into the growth of the unstable r-modes in general relativity, but the problem is still far from well understood. In particular, there are no estimates of the radiation reaction timescales for general inertial modes of a fully relativistic star in the current literature. This is an unfortunate gap since it seems reasonable to expect that barotropic models (for which all relativistic inertial modes are hybrids) will be relevant for rapidly spinning stars. In addition, it was shown in Paper I that, many of the hybrid inertial modes are unstable due to the emission of gravitational radiation. Their growth timescales were estimated using a post-Newtonian calculation and found to be considerably longer than that associated with the purely axial r-mode. It is clearly relevant to ask whether these results remain accurate when relativistic effects are included in the analysis. In fact, we have already discussed why this will not be the case for the particular axial-led mode illustrated in Figs. 5-7.

In this section, we present the first calculation of the gravitational radiation-reaction timescales for general inertial modes of fully relativistic barotropic stars. Our method for determining these timescales is significantly different from those used in yf01 ; rk02 . It is close to the post-Newtonian approach in spirit, but based on a fully relativistic analysis. When the energy radiated per cycle is small compared to the energy of the mode, $`E`$, the imaginary part of the mode frequency is accurately approximated by the expression,

$$\frac{1}{\tau }=\frac{1}{2E}\frac{dE}{dt}.$$
(54)
The rate of energy change is determined by a competition between dissipative effects such as viscosity (which tend to damp these modes) and gravitational radiation reaction (which drives the unstable modes). To estimate the contribution to $`\tau `$ from gravitational radiation we must calculate both the mode energy and the rate at which gravitational waves carry energy away from the star.

### IV.1 Calculation of the mode energy

We compute the mode energy using the Lagrangian perturbation formalism fi92 ; f78 ; fs75 . The appropriate quantity to calculate is the canonical energy, $`E_c`$. For canonical displacements $`E_c`$ agrees with the physical second order change in the energy associated with the mode to lowest order in perturbation theory fs78a . For a mode with behavior $`e^{i(\sigma t+m\phi )}`$ the canonical energy and angular momentum are related by $`E_c=(\sigma /m)J_c`$. Thus, instead of computing the canonical energy directly we perform the (much simpler) calculation of $`J_c`$ using the following expression fi92 ; f78 ; fs75 ,

$$J_c=\frac{1}{2}\mathrm{}n_\alpha \left\{U^{\alpha \beta \gamma \delta }\mathrm{\pounds }_\phi \xi _\beta ^{}_\gamma \xi _\delta +V^{\gamma \delta \alpha \beta }\mathrm{\pounds }_\phi \xi _\beta ^{}h_{\gamma \delta }\frac{1}{32\pi }ϵ^{\alpha \gamma \mu \nu }ϵ_\nu ^{\beta \delta \rho }\mathrm{\pounds }_\phi h_{\gamma \delta }^{}_\beta h_{\mu \rho }\right\}𝑑V$$
(55)
where $``$ denotes complex conjugation, $`n_\alpha =e^\nu _\alpha t`$ is the past-directed timelike normal to a $`t=\text{constant}`$ hypersurface, $`dV=\sqrt{{}_{}{}^{3}g}d^3x`$ is the volume element on that surface and

$`U^{\alpha \beta \gamma \delta }`$ $`=`$ $`(ϵ+p)u^\alpha u^\gamma q^{\beta \delta }+p(g^{\alpha \beta }g^{\gamma \delta }g^{\alpha \delta }g^{\beta \gamma })\mathrm{\Gamma }pq^{\alpha \beta }q^{\gamma \delta }`$
$`2V^{\alpha \beta \gamma \delta }`$ $`=`$ $`(ϵ+p)(u^\alpha u^\gamma q^{\beta \delta }+u^\beta u^\gamma q^{\alpha \delta }u^\alpha u^\beta q^{\gamma \delta })\mathrm{\Gamma }pq^{\alpha \beta }q^{\gamma \delta }.`$
This expression for $`J_c`$ and the corresponding expression for $`E_c`$ below are the canonical angular momentum and energy for the physical perturbation, for the real (or imaginary) part of the complex perturbation $`[h_{\alpha \beta },\xi ^\alpha ]`$. Care must be taken in comparing their values to those of other papers fi92 ; f78 ; lom98  that define $`J_c`$ (or $`E_c`$) for a complex perturbation to be the sum of its values for the real and imaginary parts.

For our slowly rotating star, we need to compute the mode energy only to leading order in $`\mathrm{\Omega }`$. Taking into account the ordering (21), a laborious but straightforward calculation gives,

$$E_c=\frac{\sigma }{m}J_c=\frac{1}{4}\sigma (ϵ+p)\left\{\kappa \mathrm{\Omega }e^\nu \xi ^\alpha \xi _\alpha ^{}+\frac{1}{2}iu^\alpha \left(\xi ^\beta h_{\alpha \beta }^{}\xi ^\beta h_{\alpha \beta }\right)\right\}𝑑V.$$
(56)
\[Note that it is necessary to use the fact that $`\xi ^\alpha `$ is a canonical displacement to reduce expression (55) to this form.\]

We now substitute for $`\xi ^\alpha `$ and $`h_{\alpha \beta }`$ their spherical harmonic expansions (10) and (17) and perform the angular integration. This leaves us with an expression for the canonical energy involving the variables $`U_l(r)`$, $`V_l(r)`$, $`W_l(r)`$ and $`h_l(r)`$ associated with the hybrid mode eigenfunctions,

$$E_c=\frac{\sigma R\overline{ϵ}}{4\kappa \mathrm{\Omega }}I,$$
(57)
with $`\overline{ϵ}`$ the average energy density and

$$I\underset{l=m}{\overset{\mathrm{}}{}}_0^1e^{\lambda \nu }\left(\frac{ϵ+p}{\overline{ϵ}}\right)\left\{\left[e^{2\lambda }\frac{4r(\nu ^{}+\lambda ^{})}{l(l+1)}\right]W_l^2+l(l+1)V_l^2+l(l+1)U_l(U_l+h_l)\right\}d\left(\frac{r}{R}\right).$$
(58)
Given one of our numerical solutions for a particular eigenmode it is straightforward to perform this radial integral numerically and compute the mode energy. It is also straightforward to show that (57) reduces to Eq. (64) of Paper I in the Newtonian limit, once we account for the difference between the inertial and the rotating frames and for the fact that (I,64) involves the sum of the real and imaginary parts of a complex perturbation, giving it an extra factor of two relative to (57).

### IV.2 Calculation of the radiated power

We now compute the rate at which energy is emitted in gravitational waves, simplifying the analysis by using a gauge-invariant expression (derived in Appendix B) for the asymptotic power radiated by an axial mode. We otherwise closely follow a calculation of Ipser ip71  (our bibliography notes two minor corrections to his equations.)

Recall that the perturbation equations we derived in Sect. II.2 are relevant only in the near zone. Because the mode frequency is proportional to the star’s angular velocity, however, the slow rotation approximation allows one to extend the near zone arbitrarily far beyond the star by making the angular velocity sufficiently small (see Table 1). Hence, to find the more general equations valid in the wave zone as well as in the near zone, we need only concern ourselves with metric perturbations of the exterior vacuum spacetime. Furthermore, the equations relevant here are those governing the perturbations of a spherical background. The leading rotational terms are of order $`\omega r`$, and outside the star the metric function $`\omega (r)`$ has the form,

$$\omega (r)=\frac{2J}{r^3},$$
(59)
with $`J`$ the angular momentum of the star h67 . Such terms can be neglected in the slow rotation approximation because

$$\omega r\left(\frac{M}{r}\right)\left(\frac{R}{r}\right)\mathrm{\Omega }R<<1.$$
(60)
Hence, the relevant components of the perturbed Einstein equation ($`\delta G_r^\phi =0`$ and $`\delta G_\theta ^\phi =0`$), are

$`0`$ $`=`$ $`i\sigma \left[h_{0,l}^{}{\displaystyle \frac{2}{r}}h_{0,l}\right]+\left[\sigma ^2{\displaystyle \frac{(l^2+l2)}{r^2}}\left(1{\displaystyle \frac{2M}{r}}\right)\right]h_{1,l}`$ (61)
$`0`$ $`=`$ $`\left(1{\displaystyle \frac{2M}{r}}\right)^2h_{1,l}^{}+\left(1{\displaystyle \frac{2M}{r}}\right){\displaystyle \frac{2M}{r}}h_{1,l}i\sigma h_{0,l},`$ (62)
where we have restored the “0” subscript on the metric function $`h_{0,l}(h_l)`$ to distinguish it from $`h_{1,l}`$.

These equations may be combined to give wave equations for either $`h_{0,l}`$ or $`h_{1,l}`$. It is not difficult to show that the resulting equation for $`h_{0,l}`$ reduces in the near zone to Eq. (44), which we used to impose boundary conditions on our interior solution. The wave equation for $`h_{1,l}`$ is,

$`0`$ $`=`$ $`\left(1{\displaystyle \frac{2M}{r}}\right)^2h_{1,l}^{\prime \prime }{\displaystyle \frac{2}{r}}\left(1{\displaystyle \frac{2M}{r}}\right)\left(1{\displaystyle \frac{5M}{r}}\right)h_{1,l}^{}`$ (63)
$`+\left[\sigma ^2{\displaystyle \frac{(l^2+l2)}{r^2}}\left(1{\displaystyle \frac{2M}{r}}\right){\displaystyle \frac{4M}{r^3}}\left(2{\displaystyle \frac{5M}{r}}\right)\right]h_{1,l}.`$
Following Ipser ip71 , we make use of the fact that the near zone extends well into the nonrelativistic region (see Table 1): Eq. (63) is easily solved in the nonrelativistic limit ($`M/r0`$) and the solution will be valid both in the wave zone and in the outer part of the near zone. Thus, we will be able to impose outgoing-wave boundary conditions at $`r\mathrm{}`$, and still match to our near-zone solution, Eq. (45). In the nonrelativistic zone, Eq. (63) is simply

$$h_{1,l}^{\prime \prime }\frac{2}{r}h_{1,l}^{}+\left[\sigma ^2\frac{(l^2+l2)}{r^2}\right]h_{1,l}=0$$
(64)
with solution

$$h_{1,l}=\frac{i\sigma ^{l+2}R^l}{(l1)(2l1)!!}r^2\left[An_l(\sigma r)+Bj_l(\sigma r)\right]$$
(65)
where $`A`$ and $`B`$ are constants, $`n_l`$ and $`j_l`$ are spherical Bessel functions and the overall normalization has been chosen for later convenience. This solution for $`h_{1,l}`$ gives rise to a solution for $`h_{0,l}`$ using Eq. (62), which, in the nonrelativistic zone becomes simply,

$$i\sigma h_{0,l}=h_{1,l}^{}.$$
(66)
In the outer part of the near zone ($`\sigma r<<1`$) we have

$$h_{0,l}A\left(\frac{R}{r}\right)^l\left\{1+\frac{B}{A}\frac{(l+2)(\sigma r)^{2l+1}}{(l1)(2l+1)[(2l1)!!]^2}\right\}\left[1+O(\sigma ^2r^2)\right].$$
(67)
In the wave zone ($`\sigma r>>1`$) we have

$$h_{0,l}\frac{(\sigma R)^{l+1}}{(l1)(2l1)!!}\left(\frac{r}{R}\right)\left[A\mathrm{sin}(\sigma rl\pi /2)+B\mathrm{cos}(\sigma rl\pi /2)\right].$$
(68)
On this wave zone solution we now impose the boundary condition that the radiation be purely outgoing. With time dependence $`e^{i\sigma t}`$ this requires that $`B=iA`$. Our solution then becomes,

$`\text{wave zone:}h_{0,l}`$ $``$ $`{\displaystyle \frac{iA(\sigma R)^{l+1}}{(l1)(2l1)!!}}\left({\displaystyle \frac{r}{R}}\right)e^{i(\sigma rl\pi /2)}`$ (69)
$`\text{near zone:}h_{0,l}`$ $``$ $`A\left({\displaystyle \frac{R}{r}}\right)^l\left\{1+{\displaystyle \frac{i(l+2)(\sigma r)^{2l+1}}{(l1)(2l+1)[(2l1)!!]^2}}\right\}\left[1+O(\sigma ^2r^2)\right].`$ (70)
The normalization constant $`A`$ is fixed by matching Eq. (70) in the outer (nonrelativistic) part of the near zone to the solution (45), which is valid throughout the near zone. This gives $`A=\widehat{h}_{l,0}`$ and justifies our neglect in Sect. II.2 of the singular solution to Eq. (44) by the fact that the singular solution is indeed of order $`(\sigma r)^{2l+1}`$ relative to the regular solution. The constant $`\widehat{h}_{l,0}`$ is, in turn, fixed by the matching to the interior solution, Eq. (48), whose normalization ultimately is set by Eq. (50).

We now have a solution for our axial-parity metric perturbation valid in the entire domain $`r[0,\mathrm{})`$. To compute the rate of energy radiation we use the gauge-invariant expression

$$\frac{dE}{dt}=\frac{1}{32\pi }(l1)l(l+1)(l+2)\underset{r\mathrm{}}{lim}\left|\frac{k_0}{r}\right|^2,$$
(71)
where

$$k_0:=h_0\frac{1}{2}_th_2.$$
(72)
In Regge-Wheeler gauge we have $`h_2=0`$, $`k_0=h_0`$, and hence we obtain

$$\frac{dE}{dt}=\frac{l(l+1)(l+2)\left(\sigma R\right)^{2l+2}|\widehat{h}_{l,0}|^2}{32\pi (l1)[(2l1)!!]^2R^2}.$$
(73)
### IV.3 The gravitational radiation reaction timescale

We now combine our expressions (57) for the mode energy and (73) for the power radiated to find the gravitational radiation reaction timescale (54). Restoring factors of $`G`$ and $`c`$, we write the timescale as

$$\frac{1}{\tau _{\text{GR}}}=\underset{l2}{}\frac{1}{\stackrel{~}{\tau }_l}\left(\frac{\mathrm{\Omega }^2}{\pi G\overline{ϵ}}\right)^{l+1}$$
(74)
where

$$\frac{1}{\stackrel{~}{\tau }_l}=\frac{c}{R}\left(\frac{3GM}{4c^2R}\right)^l\frac{l(l+1)(l+2)\kappa (\kappa m)^{2l+1}|\widehat{h}_{l,0}|^2}{16(l1)[(2l1)!!]^2I}$$
(75)
with the integral $`I`$ defined by Eq. (58).

Having obtained this expression, we are in a position where we can evaluate the growth timescales associated with the unstable inertial modes of relativistic stars. As a quick check on our calculation, we plug into Eqs. (58) and (75) our post-Newtonian solution (II,5.33)-(II,5.40) for the corrections to the $`l=m`$ Newtonian r-modes of a uniform density star. The resulting timescale associated with the $`l=m=2`$ multipole agrees with previously published results ks99 ; lf ; ak :

$$\tau _{\text{GR}}=1.56\text{}\left(\frac{1.4M_{}}{M}\right)^4\left(\frac{R}{12.53\text{km}}\right)^5\left(\frac{\pi G\overline{ϵ}}{\mathrm{\Omega }^2}\right)^3.$$
(76)
(See, for example, Table 8 of Paper I.)

As a further check on the calculation we verify that, for weakly relativistic models, the timescale (74) exhibits the same scaling with the stellar parameters as the r-mode growth timescale estimated from the Newtonian calculations. The Newtonian timescales for the $`l=m`$ r-modes exhibit the following scaling with the neutron star mass and radius ks99 ; ak :

$$\tau _{\text{GR}}=T_l\frac{GM}{c^3}\left(\frac{GM}{c^2R}\right)^{(l+3)}\left(\frac{\pi G\overline{ϵ}}{\mathrm{\Omega }^2}\right)^{(l+1)}$$
(77)
where $`T_l`$ is a constant that depends only on $`l`$ and the star’s EOS. with To test our formula for the growth timescale we keep the baryon mass, $`M_B`$, fixed at $`1.4M_{}`$ and set $`\mathrm{\Omega }^2=\pi G\overline{ϵ}`$. One would then expect a log-log plot to clearly reveal that the timescale depends on the star’s compactness as $`(M/R)^{(l+3)}`$ for low $`M/R`$. That this is, indeed, the case can be seen from the data in Fig. 11, which compares the Newtonian and relativistic growth timescales of the modes whose Newtonian analogues are the first five $`l=m`$ r-modes. Fig. 12 illustrates our results in a different way by showing the dependence on the polytropic index, $`n`$, of the timescale for the $`l=m=2`$ r-mode and its relativistic counterpart.

These two figures suggest that for highly relativistic stars, the relativistic calculation tends to give a slightly longer growth timescale than that of a Newtonian star with the same EOS, baryon mass and compactness. This suggests that general relativity tends to stabilize the modes making the r-mode instability slightly weaker than previously expected (in accordance with the results for nonbarotropic stars yf01 ; rk02 ). That this should be the case is natural: All inertial modes have relatively low frequencies since $`\sigma \mathrm{\Omega }`$. This means that the associated gravitational waves will suffer significant backscattering by the spacetime curvature as they escape to infinity. As the star becomes increasingly compact more of the curvature potential, which can be approximated by

$$V\left(1\frac{2M}{r}\right)\frac{(l1)(l+2)}{r^2}$$
(78)
cf. Eq. (63), is unveiled. Hence, one would expect low-frequency modes of oscillation to radiate less efficiently as the star becomes more compact.

We should point out that it is somewhat misleading to do a direct comparison of Newtonian and relativistic models, because there is no one-to-one correspondence between the two. Although a weakly relativistic polytrope agrees well in all of its physical characteristics (mass, radius, compactness etc.) with the Newtonian model having the same polytropic index and central density, this is not true for strongly relativistic models. Newtonian and relativistic stars with the same polytropic index can be constructed so as to agree with respect to two of their physical characteristics, but in general they will differ on the others and this will affect the growth timescales of unstable modes.

We would nevertheless like to quantify how our fully relativistic radiation timescales differ from the Newtonian ones. To account for the lack of one-to-one correspondence between Newtonian and relativistic models, we have constructed a number of different relativistic stars that agree in two of their physical characteristics with a given Newtonian polytrope of the same index, $`n`$. The fiducial Newtonian models are chosen in such a way that a $`1.4M_{}`$ star has a radius of $`12.53\text{km}`$. This facilitates comparison with results from the literature, which made use of stars with these parameters lf ; yl ; lom98 ; aks99 ; owenetal . It is natural to choose one of the physical characteristics on which the relativistic and Newtonian stars agree to be the baryon mass, $`M_B`$, (i.e. the rest mass) of the star. (For relativistic stars the baryon mass is slightly higher than the total gravitational mass, whereas these quantities are the same for Newtonian stars.) For the other physical characteristic on which the stars agree we choose such quantities as the radius, central pressure, central energy density and so on. We list the characteristics of these various models in Table 2. In the table we compare a number of relativistic uniform density ($`n=0`$) stars with our fiducial Newtonian uniform density star, labelled N0. Similiarly, we compare a set of relativistic $`n=1`$ polytropes with our fiducial Newtonian $`n=1`$ polytrope, labelled N1. All of these stars are constructed so as to have a baryon mass of $`1.4M_{}`$ and also to agree with the fiducial Newtonian model on one of their other physical characteristics. The models are listed in order of increasing compactness and the final column of the table indicates the second physical characteristic on which the stars agree.

The effect of this ambiguity in comparing the relativistic and Newtonian models is indicated in Table 3, which presents the gravitational radiation reaction timescales for the fastest growing $`l=m`$ Newtonian r-modes and their relativistic hybrid counterparts. (Some of these Newtonian timescales have been computed in previous work lom98 ; aks99 ; owenetal ; lf ; yl .) We see that, depending on which physical characteristics one chooses to equate, the timescales for the analogous mode of the Newtonian and relativistic models can differ by as much as an order of magnitude, with the relativistic mode generally having the longer growth time. The fastest relativistic growth times are obtained by equating the compactness, $`M/R`$, of the relativistic and Newtonian models. In this case, the relativistic growth times are typically weaker than the corresponding Newtonian growth times by only a factor of a few.

In connection with the results listed in Table 3 it is relevant to make two observations. First of all, the tabulated data suggest that for each $`l`$ there is a variation of about a factor of two in $`\stackrel{~}{\tau }_l`$ between the various relativistic models (for the same equation of state). It is relevant to point out that had we instead tabulated the combination $`\stackrel{~}{\tau }_lM^{l+2}/R^{l+3}`$ (i.e. accounted for the expected scaling with mass and radius) then the variation between the models would have been much smaller. For $`l=2`$ we would have found a variation of about 6% between the uniform density models, while the result for the polytropes vary by about 11%. This indicates that the variation in the growth timescale between models 1-4 and models 5-9 is mainly due to the differences in mass and radius. The second feature worth noticing from the data in Table 3 is that, while the relativistic timescales differ from the Newtonian ones by only a factor of 2-3 for the quadrupole mode, the difference increases with $`l`$. For example, for $`l=6`$ the difference is at least an order of magnitude. This result can likely be explained in terms of backscattering from the curvature potential in the exterior spacetime. From (78) we see that the “height” of the potential increases as $`(l1)(l+2)`$. Thus we would expect the difference between our fully relativistic results and the Newtonian ones to increase with $`l`$ roughly as $`(l1)(l+2)/4`$ (after rescaling with the result for the quadrupole mode). Our numerical results are in reasonable agreement with this expectation.

Finally, we want to confirm the expectation that the analogue of the $`l=m=2`$ Newtonian r-mode remains the fastest growing unstable mode also when the growth times are estimated in full general relativity. That this is the case can be seen from Table 4 where we list the growth timescales for a number of unstable inertial modes. We compare the timescales of modes of the Newtonian models (N0 and N1) with those of the corresponding modes of relativistic models 4 and 9. (These are the models with the same compactness as their Newtonian counterparts, and which lead to the fastest growth times, cf. Table 3). For each mode considered, we list (in four consecutive rows) the data for the four different stellar models — all of which have baryon mass $`M_B=1.4M_{}`$ and compactness $`M/R=0.165`$. We list the frequency of the mode in each star as well as the growth timescales associated with the various current multipole moments of the mode. We also list (enclosed in parentheses) the timescales associated with some of the inertial mode mass multipoles of the Newtonian $`n=1`$ polytrope. These mass multipoles are higher order in $`\mathrm{\Omega }`$ than can be computed within our slow rotation formalism. However, they have been calculated by Yoshida and Lee using a self-consistent third order

Newtonian formalism (see Table 4 of Ref. yl ) and we list them here simply for ease of comparison with our new results. The point is to compare the new relativistic timescales with the previously published Newtonian timescales (see, in particular, Tables 7-9 of Paper I, Table 4 of yl  and Table 1 of owenetal ).

As with the timescales listed in Table 3, most of the relativistic growth times presented in Table 4 are basically unchanged compared with their Newtonian analogs. There are, however, cases where there is a dramatic difference between the Newtonian and the relativistic results. The best example of this is provided by the fourth mode listed in Table 4. This is the axial-led inertial mode presented in Figs. 5-7 in Sect. III, whose current quadrupole moment vanishes (or nearly vanishes) in the Newtonian models. In Paper I, it was argued based on the Newtonian slow-rotation calculation that the growth of this mode is dominated by its $`l=4`$ current multipole. However, by including rotational corrections to higher order in $`\mathrm{\Omega }`$, Yoshida and Lee yl  were able to compute the $`l=3`$ mass multipole and found that it drives the mode on an even shorter timescale. Now, with the inclusion of relativistic corrections to the mode, we see that the growth of the mode is, in fact, dominated by the $`l=2`$ current multipole — as one would have expected in the first place. This mode is significantly more unstable in general relativity than in the Newtonian calculations. However, it’s growth time is nevertheless still much longer than that of the mode whose Newtonian analogue is the $`l=m=2`$ r-mode.

To conclude: The data presented in Tables 3 and 4 suggest that the relativistic corrections to the timescales are not large enough to alter the standard picture of the gravitational-wave driven instability of sufficiently hot and rapidly rotating neutron stars. The fastest growing mode (by at least an order of magnitude) is the axial-led inertial mode corresponding to the $`l=m=2`$ Newtonian r-mode. And although the actual growth timescale of this mode is uncertain due to the uncertainty in the neutron star EOS, it is unlikely to be significantly shorter than has been estimated here.

## V Concluding remarks

In this paper we have studied the inertial modes of rotating relativistic stars. Numerical results were presented for the mode-eigenfrequencies and the associated eigenfunctions of barotropic models. These results were shown to be in good agreement with results in the literature (in particular in the post-Newtonian limit), and provide a significant improvement on previous studies as far as the strongly relativistic regime is concerned.

We also analyzed the rate at which these modes radiate gravitationally. In particular, we studied the growth timescale of various modes that are unstable due to the emission of gravitational waves fs78a ; fs78b . Our calculation was based on two ingredients: The energy associated with the mode oscillation was determined as the canonical energy defined by Friedman f78 , while the gravitational-wave luminosity followed from an analysis parallel to that of Ipser ip71 . By combining these two quantities we arrived at the required damping/growth timescale due to gravitational-wave emission. Our approach to the problem is novel, and differs from the methods previously used to estimate the corresponding timescales for the modes of non-barotropic stars rk02 ; yf01 . Nevertheless, the final results show the same qualitative behaviour. Most notably, the post-Newtonian estimates for the growth rate of an unstable inertial mode are found to be surprisingly accurate even for strongly elativistic models. Still, the results show a deviation from the post-Newtonian estimates as the star reaches compactness similar to that expected of a neutron star $`M/R0.1`$. Then the efficiency of radiation reaction tends to decrease. This result is likely due to the fact that the low-frequency waves from an inertial mode experience enhanced backscattering by the spacetime curvature as the star becomes increasingly compact (recall that the spacetime of a spherical star has a curvature potential barrier with a peak in the region $`R/M3`$). This means that general relativistic effects tend to stabilize the inertial modes. This is in contrast to the results for the instability associated with the acoustic f-modes. As was shown by Stergioulas and Friedman sf , the f-modes are significantly destabilized by relativistic effects.

Despite significant progress in the last few years, it is still not clear to what extent the gravitational-wave driven instabilities in rotating compact stars are of astrophysical relevance. In particular, we do not yet have a clear answer to the question of whether the unstable modes may lead to detectable graviational waves or whether they limit the spin of nascent neutron stars or of old neutron stars spun up by accretion. However, it is important to realize that serious theoretical challenges need to be overcome if we want to make further progress in this area of research. For the unstable r-modes, key questions include the role of complex (not well understood) interior physics, e.g. the strength of hyperon bulk viscosity pbj ; lo02 , and the saturation amplitude set by nonlinear coupling arras .

## Appendix A The spectral method used to solve the eigenvalue problem

To solve the set of equations (22)-(38) for the inertial modes of a barotropic relativistic star, we use a variant of the method developed in Paper I for the analogous Newtonian hybrid/inertial mode problem lf . We express our equilibrium and perturbation variables as a sum over a set of basis functions and substitute these series into our system of differential equations. This results in a system of algebraic equations for the expansion coefficients, which may then be solved using standard linear algebra techniques. In Paper I, we made use of power series expansions and were able to accurately compute the Newtonian hybrid mode eigenvalues and eigenfunctions. However, for the relativistic problem we have found it necessary to use a basis of orthogonal functions and we work instead with Chebychev polynomials (or, more specifically, “Type I” Chebychev polynomials arfken ). In other words, we use a spectral method to approach the problem.

Since the Chebychev polynomials are naturally defined on the domain $`[1,1]`$, we define a new coordinate

$$y=2\left(\frac{r}{R}\right)1$$
(79)
which maps the interior of the star, $`0rR`$, to the required domain. The Chebychev polynomial of degree $`i`$ is then given by

$$T_i(y)=\mathrm{cos}(i\mathrm{arccos}y).$$
(80)
We begin by expressing our equilibrium variables in terms of this Chebychev basis. We construct an equilibrium star using standard integration recipes and then use Chebychev approximation numrec  to find a Chebychev series that accurately fits each of our equilibrium variables. In other words, we take all of the background variables appearing in Eqs. (23)-(38), such as $`ϵ(r)`$, $`p(r)`$, $`\overline{\omega }(r)`$ etc., and represent each of them by a Chebychev series of the form,

$$B(r)=\underset{i=0}{\overset{\mathrm{}}{}}b_iT_i(y)\frac{1}{2}b_0,$$
(81)
where the coefficents, $`b_i`$, are determined from the Chebychev approximation algorithm (see Numerical Recipes numrec ).

Having constructed Chebychev series for our (known) background variables, we write each of our (unknown) perturbation variables in terms of an expansion in Chebychev polynomials of the form,

$$F_l(r)=\left(\frac{r}{R}\right)^{l+q}\left[\underset{i=0}{\overset{\mathrm{}}{}}f_{l,i}T_i(y)\frac{1}{2}f_{l,0}\right].$$
(82)
The factor $`(r/R)^{l+q}`$ provides for the condition of regularity at the origin, which requires the perturbation variables to vanish as an appropriate power of $`r`$ as $`r0`$. The axial parity variables $`U_l(r)`$ and $`h_l(r)`$ have $`q=1`$, while the polar parity variables $`W_l(r)`$ and $`V_l(r)`$ have $`q=0`$.

If we substitute the various Chebychev series represented by (81) and (82) into our perturbation equations, each term in these equations will take the form of a product of two Chebychev series. That is, the generic term in our perturbation equations will have the form $`B(r)F_l(r)`$ with $`B(r)`$ a known background function and $`F_l(r)`$ an unknown perturbation variable. We would like to be able to write such a product as a new expansion in Chebychev polynomials. This can be accomplished using the identity $`2T_iT_j=T_{i+j}+T_{|ij|}`$, which follows from Eq. (80) and standard cosine identities. After some careful rearrangement of terms we find for the product of the series (81) and (82),

$$B(r)F_l(r)=\frac{1}{2}\left(\frac{r}{R}\right)^{l+q}\left[\underset{i=0}{\overset{\mathrm{}}{}}\pi _{l,i}T_i(y)\frac{1}{2}\pi _{l,0}\right]$$
(83)
where

$$\pi _{l,i}=\underset{j=0}{\overset{\mathrm{}}{}}\left[b_{i+j}+\mathrm{\Theta }(j1)b_{|ij|}\right]f_{l,j}$$
(84)
with

$$\mathrm{\Theta }(k)=\{\begin{array}{c}0\text{for }k<0\hfill \\ 1\text{for }k0\hfill \end{array}.$$
(85)
We also need an expression for the derivatives of our perturbation variables in terms of the Chebychev expansions (82). If we define the Chebychev expansion for the derivative of $`F_l(r)`$ as follows:

$$R\frac{d}{dr}\left[\left(\frac{r}{R}\right)^{(l+q)}F_l\right]\underset{i=0}{\overset{\mathrm{}}{}}\stackrel{~}{f}_{l,i}T_i(y)\frac{1}{2}\stackrel{~}{f}_{l,0}$$
(86)
and then make use of standard identities involving Chebychev polynomials arfken  it is not too difficult to show that the coefficients $`\stackrel{~}{f}_{l,i}`$ of this series are related to the coefficients $`f_{l,i}`$ of (82) by

$$\stackrel{~}{f}_{l,i}\stackrel{~}{f}_{l,i+2}=4(i+1)f_{l,i+1}.$$
(87)
Our method of solution is as follows: We expand all of the quantities appearing in Eqs. (23)-(38) in Chebychev series, and substitute these expansions into the equations and into the boundary and matching conditions. We then use the formulas (83) and (87) to express the resulting equations as a linear algebraic system of the form

$$Ax=0$$
(88)
where $`A`$ is a known matrix that depends nonlinearly on the parameter $`\kappa `$ and $`x`$ is a vector whose components are the unknown coefficients in the Chebychev series for the variables $`h_l`$, $`U_l`$, $`V_l`$, $`W_l`$ and their derivatives.

To satisfy Eq. (88) we must search for those values of $`\kappa `$ for which the matrix $`A`$ is singular; that is, we must find the zeroes of the determinant of $`A(\kappa )`$. Since $`A`$ is infinite dimensional, we must truncate our spherical harmonic expansions (10) and (17) at some maximum index $`l_{\text{max}}`$ and also truncate our Chebychev expansions (81), (82) and (86) at some maximum index $`i_{\text{max}}`$. The resulting finite matrix is band diagonal. To find its zeroes we use standard linear algebra and root finding routines. We then check for convergence of these eigenvalues as we increase $`l_{\text{max}}`$ and $`i_{\text{max}}`$.

The eigenfunctions associated with these eigenvalues are determined by the perturbation equations only up to normalization. Given a particular eigenvalue, we find its associated eigenfunction by replacing one of the equations in the system (88) with the normalization condition (50). Since this eliminates one of the rows of the singular matrix $`A`$ in favor of the normalization equation, the result is an algebraic system of the form

$$\stackrel{~}{A}x=b,$$
(89)
where $`\stackrel{~}{A}`$ is now a non-singular matrix and $`b`$ is a known column vector. We solve this system for the vector $`x`$ and reconstruct the various Chebychev expansions from this solution vector of coefficients.

## Appendix B Gauge-invariant expression for the luminosity

We present here a brief derivation of the gauge-invariant expression (71) for the rate at which energy is radiated to future null infinity ($`^+`$) by an axial mode. One first writes the energy radiated in an asymptotically regular gauge and then observes that the gauge-dependent quantity in the expression can be replaced by a quantity that is gauge invariant.

We begin with standard expressions for $`dE/dt`$. Written in terms of the leading term $`\sigma ^0`$ in the asymptotic shear of outgoing null geodesics, in outgoing null coordinates $`(u,r,\theta ,\phi )`$ of a flat asymptotic metric $`\eta _{\alpha \beta }`$, the instantaneous power radiated is

$$\frac{dE}{dt}=\frac{1}{4\pi }_{\mathrm{}}\left|_t\sigma ^0(u,\theta ,\phi )\right|^2𝑑\mathrm{\Omega }.$$
(90)
where $`{\displaystyle _{\mathrm{}}}:=\underset{r\mathrm{}}{lim}{\displaystyle _{S_r}}`$, with $`S_r`$ a sphere of constant $`r`$ and $`u`$. In an asymptotically regular gauge ip71 ; f78 , the components $`h_{(\mu )(\nu )}`$ of the metric perturbation in an orthonormal basis $`\{𝐞_{(\mu )}\}`$ fall off like $`1/r`$ or faster in an outgoing null direction, and (for a real perturbation), Eq. (B1) takes the form,

$$\frac{dE}{dt}=\frac{1}{16\pi }_{\mathrm{}}\left[\left(_th_{(\theta )(\phi )}\right)^2+\frac{1}{4}\left(_th_{(\theta )(\theta )}_th_{(\phi )(\phi )}\right)^2\right]𝑑\mathrm{\Omega }.$$
(91)
Because this expression is identical to the flat-space expression in terms of the deviation $`h_{\alpha \beta }`$ of the metric from its Minkowski value, the radiation field can be treated as a perturbation of flat space. In particular, the parts of $`h_{\alpha \beta }`$ belonging to different representations of the rotation group decouple, and we can restrict consideration to perturbations $`h_{\alpha \beta }`$ belonging to an $`(l,m)`$ representation with axial symmetry. The gauge invariance of the vector $`k_\alpha `$ in (71) arises from this decoupling.

In Regge-Wheeler notation rw57 , a (complex) axial $`(l,m)`$ perturbation has the form

$`h_{tA}=h_0\mathrm{\Phi }_{lA}^m,h_{rA}=h_1\mathrm{\Phi }_{lA}^m,h_{AB}=h_2\chi _{lAB}^m`$ (92)
with $`A,B`$ indices on the sphere. Here, as in Eq. (17), the axial vector harmonics have the form $`\mathrm{\Phi }_{lA}^m=ϵ_A^B_BY_l^m`$, with $`ϵ_\theta ^\phi =\frac{1}{\mathrm{sin}\theta }`$, $`ϵ_\phi ^\theta =\mathrm{sin}\theta `$, while the axial tensor harmonics have components

$$\chi _{\theta \theta }=\frac{1}{\mathrm{sin}^2\theta }\chi _{\phi \phi }=\frac{1}{\mathrm{sin}\theta }(_\theta \mathrm{cot}\theta )_\phi Y_l^m\chi _{\theta \phi }=\frac{1}{2}\mathrm{sin}\theta (_\theta ^2\mathrm{cot}\theta _\phi \frac{1}{\mathrm{sin}^2\theta }_\phi ^2)Y_l^m.$$
(93)
A gauge transformation associated with a vector field $`\zeta ^\alpha `$ changes $`h_{\alpha \beta }`$ by $`_\alpha \zeta _\beta +_\beta \zeta _\alpha `$. The transformation preserves the axial $`(l,m)`$ representation to which the physical perturbation belongs if and only if $`\zeta _\alpha `$ is itself an axial vector on the sphere: $`\zeta _t=\zeta _r=0,`$

$$\zeta _A=\zeta \mathrm{\Phi }_{lA}^m$$
(94)
Then
xxxxxxxxxxxxxxx $`h_{tA}h_{tA}+_t\zeta \mathrm{\Phi }_{lA}^m,h_{rA}h_{rA}+r^2_r({\displaystyle \frac{1}{r^2}}\zeta )\mathrm{\Phi }_{lA}^m,h_{AB}h_{AB}+\zeta \chi _{lAB}^m,`$
and

$$h_0h_0_t\zeta ,h_1h_1r^2_r(\frac{1}{r^2}\zeta ),h_2h_2+2\zeta ,$$
(95)
implying gs ; gm  that the vector $`k_\alpha `$, whose nonzero components are

$$k_0=h_0+\frac{1}{2}_th_2,k_1=h_1+\frac{1}{2}r^2_r(r^2h_2),$$
(96)
is gauge invariant.

When an axial perturbation is written in an asymptotically regular gauge, expression (91) is valid. For a complex perturbation,

$`{\displaystyle \frac{dE}{dt}}(h)`$ $`:=`$ $`{\displaystyle \frac{dE}{dt}}(\mathrm{}h)+{\displaystyle \frac{dE}{dt}}(\mathrm{}h)`$ (97)
$`=`$ $`{\displaystyle \frac{1}{16\pi }}{\displaystyle _{\mathrm{}}}\left|_th_2\right|^2\left[\left|\chi _{(\theta )(\phi )}\right|^2+{\displaystyle \frac{1}{4}}\left|\chi _{(\theta )(\theta )}\chi _{(\phi )(\phi )}\right|^2\right]𝑑\mathrm{\Omega }`$
$`=`$ $`{\displaystyle \frac{1}{16\pi }}(l1)l(l+1)(l+2)\underset{r\mathrm{}}{lim}\left|{\displaystyle \frac{_th_2}{2r}}\right|^2.`$
Finally, noting that, in this gauge, $`h_2=O(r)`$, $`h_0=O(1)`$, we have $`{\displaystyle \frac{_th_2}{2r}}={\displaystyle \frac{k_0}{r}}+\text{O}(r^2).`$ Thus

$$\frac{dE}{dt}(h)=\frac{1}{16\pi }(l1)l(l+1)(l+2)\underset{r\mathrm{}}{lim}\left|\frac{k_0}{r}\right|^2,$$
(98)
for a complex perturbation. A real axial mode then radiates average power

$$\frac{dE}{dt}(\mathrm{}h)=\frac{dE}{dt}(\mathrm{}h)=\frac{1}{32\pi }(l1)l(l+1)(l+2)\underset{r\mathrm{}}{lim}\left|\frac{k_0}{r}\right|^2,$$
(99)
with the limit, as above, taken along a radially outgoing null direction. Although, for an unstable mode, the perturbed metric and its gauge invariant quantities blow up exponentially at spatial infinity, the power radiated to null infinity is well-defined and finite.

###### Acknowledgements.

We are grateful to the Institute for Theoretical Physics at the University of California - Santa Barbara for generous hospitality during the workshop “Spin and Magnetism in Young Neutron Stars” where part of this research was conducted. KHL acknowledges with thanks the support provided by a Fortner Research Fellowship at the University of Illinois and by the Eberly research funds of the Pennsylvania State University. This research was also supported in part by NSF grant AST00-96399 at Illinois and NSF grant PHY00-90091 at Pennsylvania. JLF was supported in part by NSF Grant PHY00-71044. NA acknowledges support from the Leverhulme Trust via a Prize Fellowship, PPARC grant PPA/G/1998/00606 as well as the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).