text:
# Packing-Limited Growth

## I Introduction

The dynamic packing of objects is an often overlooked variation on the theme of static packing. Given a mechanism for the creation, growth, movement, and interaction of like objects in a given dimension, what structures result? Here, we find universal features of a simple yet broad class of such mechanisms falling under the rubric “packing-limited growth” (PLG). A PLG mechanism entails objects being seeded randomly, growing according to a rule that may be specific to each object, and stopping when a part of their boundary hits that of another object. A motivating physical example of this kind of pattern formation may be found in the competition between tree crowns in dense forests Horn (1971); Takahashi (1996), the structure of porous media van der Marck (1995); Turcotte (1997); Hecht (2000), and the generalized problem of dense packings Conway and Sloane (1999).

A simple model of PLG has previously been studied by Andrienko, Brilliantov, and Krapivsky Andrienko et al. (1994); Brilliantov et al. (1994) (the ABK model). In their setting, spheres are seeded randomly in space and time. A sphere’s radius increases linearly with time, halting when another sphere is touched. This model is amenable to an approximate analysis for $`d>1`$ Brilliantov et al. (1994) and has an exact solution in $`d=1`$. In this present work, we are interested in the limiting distribution of radii, $`N(r)`$. For PLG models, we expect $`N(r)r^\alpha `$ for small $`r`$. Note that the fractal dimension $`D`$ of the set comprising all sphere centers, which is often measured instead of $`\alpha `$, is related to $`\alpha `$ as $`D=\alpha 1`$ Manna (1992). In $`d=1`$, the exact solution gives $`\alpha =1`$, which corresponds to $`D=0`$ (meaning the number of centers diverges logarithmically) Brilliantov et al. (1994). For $`d=2`$, Andrienko et al. Andrienko et al. (1994) report that $`D1.75`$ and hence $`\alpha 2.75`$ based on numerical evidence. Elsewhere, the same authors Brilliantov et al. (1994) determine numerically that $`\alpha 2.53`$.

We claim that the actual value of $`\alpha `$ is essentially independent of the specifics of the growth dynamics. To see why, we examine a related $`d=2`$ random packing model due to Manna Manna (1992). The packing process begins in a finite-sized volume (or alternatively with an initial population of fixed radius disks). New disks are added one at a time by randomly choosing a point in the packing’s pore space and centering there the largest possible non-overlapping disk. Manna finds $`\alpha 2.62`$ and $`\alpha 2.64`$ for two example packings.

We refer to the Manna model as “random Apollonian packing” (RAP) since it may be seen as a variation of the well known Apollonian packing Mandelbrot (1983). For the $`d=2`$ version of the latter, pore spaces are always formed by three disks each touching the other two and disks are added so as to fill the pore space fully (i.e., the inserted disk touches all three of its surrounding neighbors). Numerics give $`\alpha 2.31`$ for Apollonian packing in agreement with analytically derived bounds Manna and Hermann (1991). Aste Aste (1996) has shown for static polydisperse packings that are space filling, Apollonian packing provides the lowest bound on $`\alpha `$ while the upper bound is $`\alpha =d+1`$.

There are two important observations to make here. First, as noted by Brilliantov et al. Brilliantov et al. (1994), the RAP model is the ABK model in the case of infinitely fast growth: as soon as a sphere is nucleated, it instantaneously expands to hit the nearest sphere boundary. Second, for general PLG models, pore spaces evidently increase in number and decrease in size with time. This means that the likelihood of two spheres nucleating in the same pore space also decreases. Eventually the mechanism of the RAP model must take over, and all collisions will be between a growing sphere in a pore space and a “stuck” sphere forming part of the pore space boundary. Thus, the RAP model is not just a curious end point of the ABK model but, in fact, entails the sole mechanism describing how small radius spheres pack. We suggest then that measurements of $`\alpha `$ in the ABK model should coincide with revised estimates from simulations of the RAP model. In the remainder of the paper, we describe a general class of PLG models for which $`\alpha `$ is invariant, provide a simple solution for the $`d=1`$ problem, develop a scaling theory that describes how the radii distribution and volume of pore space evolve with time, and provide results from extensive numerical simulations.

## II General model of packing-limited growth

For our general conception of PLG we take $`d`$-dimensional spheres growing in a volume $`V`$. Spheres are nucleated randomly in space at a rate $`\kappa `$ per unit volume surviving only if they are injected into the unoccupied pore space. All spheres stop growth upon contact with a neighboring sphere. Figure 1 provides a visual context for the dynamics explained below. The $`i`$th sphere, which is initiated at $`t_i`$, grows at a rate $`G_i(tt_i)`$ giving the radius $`r_i(t)`$ as

$$r_i(t)=_0^{tt_i}G_i(u)\text{d}u,$$
(1)
for $`tt_i`$. We assume the weak requirement that each sphere grows in a strict monotonic fashion, i.e., for each $`i`$, $`G_i(tt_i)ϵ>0`$. The model is thus one of arbitrary individual growth limited solely by packing. The dynamics continue in the above fashion until the volume $`V`$ is completely filled and a final jammed state is reached.

Consider an individual pore $`\mathrm{\Gamma }`$ of diameter $`\lambda `$ (i.e., the maximum separation of two boundary points that can be joined by a line contained entirely within the pore). The rate $`\nu `$ at which spheres nucleate in $`\mathrm{\Gamma }`$ (excluding growth for the moment) is given by

$$\nu \kappa V_d\lambda ^d,$$
(2)
where $`V_d`$ is the volume of $`d`$-dimensional sphere of unit radius. The typical time $`\tau `$ between nucleations in $`\mathrm{\Gamma }`$ is therefore

$$\tau =1/\nu \kappa ^1V_d^1\lambda ^d.$$
(3)
On the other hand, if sphere $`i`$ nucleates in $`\mathrm{\Gamma }`$, it will jam in a time $`\tau _{\mathrm{jam}}`$ that can be no greater than the time it takes for its radius to reach $`\lambda /2`$. Together with the assumption that $`G_i(tt_i)ϵ>0`$, this gives

$$\tau _{\mathrm{jam}}\lambda /2ϵ.$$
(4)
So when $`\tau _{\mathrm{jam}}\tau `$, i.e., when

$$V_d\lambda ^{d+1}\kappa /2ϵ1,$$
(5)
we expect the packing mechanism to be the same as the RAP model—all spheres will be stopped by an existing stopped sphere and never by another moving sphere. Growth rate therefore becomes irrelevant as far as the final packing is concerned. From Eq. (5), we see that the general model reduces to the RAP model when the maximum pore size $`\lambda _{\mathrm{max}}`$ satisfies

$$\lambda _{\mathrm{max}}\left(2ϵ/V_d\kappa \right)^{1/(d+1)}.$$
(6)
Since $`\lambda _{\mathrm{max}}`$ steadily decreases with $`t`$ (i.e., space fills up) this condition will always be eventually satisfied.

## III Exact $`d=1`$ solution

As we have noted, the ABK model has been solved exactly in the $`d=1`$ case Brilliantov et al. (1994) with the outcome being $`\alpha =1`$. Here, we achieve the same result with a considerably simpler calculation. Since we have posited that $`\alpha `$ is a universal exponent for a general class of packing-limited growth models, we may choose the straightforward example of the $`d=1`$ RAP model. We take the unit interval $`[0,1]`$ as the initial vacant pore space. Spheres are now line segments and are limited either on the left or right as they are added with the points 0 and 1 providing initial boundaries. Thus, the unit interval is filled in from each end with the one solitary pore space diminishing in the middle. Writing the length of the $`n`$th line segment as $`l_n`$ we have

$$l_n=z_{n1}\left(1\underset{i=1}{\overset{n1}{}}l_i\right),$$
(7)
where $`z_i`$ is a random number uniformly distributed on the unit interval and $`l_1=z_0`$. In other words, each new line segment is a random fraction of the current pore space. Iteratively solving Eq. (7) gives

$$l_n=z_{n1}\underset{i=0}{\overset{n2}{}}(1z_i).$$
(8)
We can in principle find the distributions of the $`l_i`$ but for our present purposes their means are sufficient and we have

$$l_n=2^n.$$
(9)
We therefore expect the typical number of line segments with $`l^{}l=2^n`$ to be

$$N(l^{}l=2^n)=n=\frac{\mathrm{ln}l}{\mathrm{ln}2}.$$
(10)
Since this is the cumulative frequency distribution, i.e., $`N(l^{}l)=_l^{\mathrm{}}N(l^{})\text{d}l^{}`$, we find the frequency distribution $`N(l)`$ must behave as

$$N(l)\frac{1}{\mathrm{ln}2}l^1,$$
(11)
yielding, as expected, $`\alpha =1`$.

## IV Description of insertion probability and pore space volume

We next investigate the form of $`P_{\mathrm{ins}}(r;n)`$, the probability distribution of the $`(n+1)`$th sphere’s radius to be inserted into a packing. In addition, we characterize $`P(r;n)`$, the probability distribution of sphere radii after $`n`$ spheres have been packed. Note that since we have shown that a general class of PLG models reduce to the RAP model, we now use the number of spheres $`n`$ rather than time $`t`$ to index our quantities.

We first note that the distribution $`P(r;n)`$ fills in from the right with increasing $`n`$. As the packing fills in space, the maximum pore size either decreases or does not change and so does the maximum size of any sphere that may be added. We write the radius of this largest sphere as $`r_c`$. We observe that to a first approximation, $`P(r;n)`$ above this cutoff scale $`r_c`$ follows its limiting power law form while below the distribution it is essentially flat (see Fig. 2). For the purposes of estimation, we assume this form exactly as

$$P(r;n)=\{\begin{array}{cc}\frac{\alpha 1}{\alpha }r_c^1\hfill & \text{for }r<r_c\hfill \\ \frac{\alpha 1}{\alpha }r_c^1\left(\frac{r}{r_c}\right)^\alpha \hfill & \text{for }rr_c.\hfill \end{array}$$
(12)
The corresponding frequency distribution is $`N(r;n)P(r;n)`$. Since the tail of $`N(r;n)`$ remains fixed as $`n`$ increases, we can obtain the scaling of $`r_c`$ with $`n`$. From Eq. (12), the tail of $`N(r;n)`$ behaves as

$$N(r;n)=\frac{\alpha 1}{\alpha }nr_c^{\alpha 1}r^\alpha =kr^\alpha ,$$
(13)
for $`r>r_c`$. Since the prefactor $`k`$ must be constant, we have

$$r_c=\left(\frac{\alpha 1}{\alpha k}n\right)^{1/(\alpha 1)}n^{1/(\alpha 1)}.$$
(14)
The uniformity of $`P(r;n)`$ for $`r<r_c`$ closely relates to the form of $`P_{\mathrm{ins}}(r;n)`$. It is possible to write down an exact expression for $`P_{\mathrm{ins}}(r;n)`$, that is,

$$P_{\mathrm{ins}}(r;n)=\frac{dV\delta (D_n(\stackrel{}{x})r)}{dV},$$
(15)
where the integrals are over the pore space, and $`D_n(\stackrel{}{x})`$ is the distance from the point $`\stackrel{}{x}`$ to the closest pore space boundary after $`n`$ spheres have been inserted. This integral may be solved exactly in the limit of very small radii,

$$\underset{r0}{lim}P_{\mathrm{ins}}(r;n)=\frac{S(n)}{\mathrm{\Phi }(n)},$$
(16)
where $`S(n)`$ and $`\mathrm{\Phi }(n)`$ are the surface area and available pore space of the existing $`n`$ packed spheres. This means that the region of pore space available for insertion of an infinitesimal sphere is proportional to the surface area of the extant packing. Assuming that this holds approximately for all spheres below the cutoff, the full distribution $`P_{\mathrm{ins}}(r;n)`$ may be modeled as a purely flat distribution,

$$P_{\mathrm{ins}}(r;n)=\{\begin{array}{cc}r_c^1\hfill & \text{for }r<r_c\hfill \\ 0\hfill & \text{for }rr_c.\hfill \end{array}$$
(17)
Good agreement with this approximation is shown in the inset of Fig. 2. Comparing Eq. (16) with Eq. (17), we see that $`r_c=\mathrm{\Phi }(n)/S(n)`$. We calculate $`S(n)`$ and $`\mathrm{\Phi }(n)`$ using $`N(r;n)`$,

$$S(n)=\frac{k\alpha V_d}{\alpha d}r_c^{(\alpha d)},$$
(18)
and

$$\mathrm{\Phi }(n)=\frac{k\alpha V_d}{(d+1)(d+1\alpha )}r_c^{d+1\alpha },$$
(19)
where as before $`V_d`$ is the volume of a unit radius sphere in $`d`$ dimensions. Note that $`S(n)\mathrm{}`$ and $`\mathrm{\Phi }(n)0`$ as $`n\mathrm{}`$.

Using Eqs. (18) and (19), and the result $`r_c=\mathrm{\Phi }(n)/S(n)`$, we find an estimate of $`\alpha `$ as

$$\widehat{\alpha }=d+\frac{d+1}{d+2}.$$
(20)
In contrast, a modified calculation for the ABK model finds Brilliantov et al. (1994)

$$\alpha _{\mathrm{mf}}=1+d\left(1+\mathrm{exp}\left[2\frac{2^{d+2}2}{d+2}\right]\right).$$
(21)
Note that $`\widehat{\alpha }`$ must be an upper bound on the true value of $`\alpha `$ for normal Apollonian packing. Furthermore, due to the nature of the approximations made in forming $`P_{\mathrm{ins}}(r;n)`$ and $`P(r;n)`$, $`\widehat{\alpha }`$ is also an upper bound for the exponent of the RAP model. As we show in the following section when we consider our numerical results, both of these predictions of $`\alpha `$ appear to hold for certain (mutually exclusive) ranges of $`d`$.

Using Eqs. (14), (18), and (19), we are able to find the scalings of surface area and pore space volume with $`n`$,

$$S(n)n^\gamma n^{(\alpha d)/(\alpha 1)},$$
(22)
and

$$\mathrm{\Phi }(n)n^\beta n^{(d+1\alpha )/(\alpha 1)}.$$
(23)
These relations provide us with further methods for determining $`\alpha `$ and for testing scaling relations between the iterative and structural nature of the packing. Equation 23, in particular, affords a more robust measurement than does obtaining $`\alpha `$ directly from $`P(r;n)`$, and we employ this fact in our numerical investigations.

Note that for comparison between the current theory and that of ABK we must convert scaling predictions that depend on time to those that depend on sphere number. Specifically, the ABK theory predicts $`\mathrm{\Phi }(t)t^A`$ where $`A=\mathrm{exp}[2(2^{d+2}2)(d+2)]`$. Because the ABK model can be seen as an RAP model with infinite growth velocity, the time between events is inversely proportional to the pore space available (given a uniform rate for attempted nucleation $`\kappa `$). So the result for the decay of pore space from ABK may be rewritten as $`\mathrm{\Phi }_{ABK}(n)n^\beta ^{}`$, where

$$t_n=\frac{1}{\kappa }\underset{i=1}{\overset{n}{}}\frac{1}{\mathrm{\Phi }_{ABK}(i)}.$$
(24)
This implies that $`t_nn^{1+\beta ^{}}`$. So if ABK predicts $`\mathrm{\Phi }(t)t^A`$, this is equivalent to writing $`\mathrm{\Phi }_{ABK}(n)n^{A(1+\beta ^{})}`$. Equating the two forms leads to the conclusion that

$$\beta ^{}=\frac{A}{1A}.$$
(25)
This will be useful in the following section when the scaling of pore space decay is examined. Finally, key scaling relations are summarized in Table 1.

## V Numerical results

The numerical procedure for all variants of PLG schemes obey the same basic algorithm. Spheres are seeded according to a Poisson distribution, with rate $`\kappa `$ and only grow within the pore space. Once injected, spheres grow according to model rules (linear velocity, heterogeneous, exponential, etc.) until they collide and are jammed. In the case of the RAP model, only one sphere is allowed to nucleate and expand at any given point in the simulation. Using this event driven procedure we are able to calculate limiting states and determine improved estimates of the universal value of $`\alpha `$ for $`d=2`$, $`3`$, and $`4`$. All packings are simulated on a unit hypercube with periodic boundary conditions.

Figure 3 shows an example packing using the RAP model in $`d=2`$. It is clear that the number of spheres will increase without bound and that the pore space ultimately vanishes. The associations with the Apollonian case are visually evident.

To demonstrate the universality of $`\alpha `$, we consider four different PLG models in $`d=2`$ with a variety of initial conditions, the details of which are given in Table 2. The frequency distribution $`N(r)`$ for these four PLG models is shown in Figure 4. The main plot shows the distributions recentered using their respective means. The recentered distributions are indistinguishable to the eye clearly indicating that the specifics of the growth mechanisms are irrelevant and that the exponent $`\alpha `$ is a universal one.

The scaling theory developed in the previous section relates the structure of a packing (its fractal dimension) to the change of basic elements (critical radius, surface area, volume, etc.) as a function of iteration. In order to test the scaling theory, as well as the specific prediction $`\widehat{\alpha }=d+(d+1)/(d+2)`$ in higher dimensions, we now examine iterative and structural scalings of the RAP model in $`d=2`$, 3, and 4.

We estimate $`\alpha `$ using a variety of means and the results are summarized in Table 3. The agreement between the predicted value of $`\alpha `$ calculated using $`\mathrm{\Phi }(n)`$ and that calculated directly from the geometry of the packing provides further proof that the scaling theory outlined in the previous section is valid, regardless of the specific value of $`\alpha `$. The current theory appears to improve with increasing dimension, whereas the approximate theory of ABK only holds in $`d=2`$.

Both the current theory and that of ABK predict the scaling of pore space as a function of iteration or time. Evaluating the pore space decay $`\mathrm{\Phi }(n)`$ in $`d=2`$ suggests that the agreement with ABK is most likely coincidental. As evidenced by the data in Fig. 5 we are able to establish $`95\%`$ confidence intervals that exclude the predictions of ABK in $`d=2`$. This is not altogether surprising as ABK is an approximate theory describing the collision of moving spheres whereas in the RAP procedure all collisions are between a sphere and the presently jammed state. In higher dimensions, the theory of ABK fails, and the current theory becomes increasingly appropriate, falling within the numerically defined confidence intervals for $`\beta `$ in $`d=4`$. The results are summarized in Table 4.

## VI Conclusion

The iterative nature of the random Apollonian packing model and its simple prescription for packing spheres suggests that at sufficiently small scales, the structure of packing-limited growth models is essentially unaffected by initial conditions or dynamics. The scaling of pore space, surface area, number, and critical radius $`r_c`$ are all interrelated and can be expressed in terms of simple power laws using only the dimension $`d`$ and the exponent $`\alpha `$. Extensive numeric simulations demonstrate the validity of the predicted upper bound $`\widehat{\alpha }`$. Although $`\widehat{\alpha }`$ appears to be an overestimate of the actual value of $`\alpha `$ in $`d=2`$ and $`d=3`$, it agrees with simulations in $`d=4`$, and presumably beyond. A refined approximation of the insertion probability $`P_{\mathrm{ins}}(r;n)`$ seems to be the best approach by which to improve the estimate of $`\alpha `$. The applicability of predictions of the present model and its relevance to physical and biological problems may lie in the process of balancing the idea of packing-limited growth with other dynamical possibilities such as aggregation, competition, and death.

###### Acknowledgements.

This work was supported in part by NSF Grant No. EAR-9706220.


format: text

text:
# Inertial modes of slowly rotating relativistic stars in the Cowling approximation

## 1 INTRODUCTION

Non-radial oscillations of relativistic stars gained a lot of interest in the last decades because of the possible detection of their associated gravitational waves. Especially after the discovery that the $`r`$ modes of rotating neutron stars are generically unstable to gravitational radiation (Andersson 1998, Friedman & Morsink 1998), the interest in the rotational instabilities renewed. Great effort has been put into understanding this instability and its implications on the spin evolution of both newly born and old accreting neutron stars. If not damped by viscous effects, the $`r`$ mode instability should be a very efficient way to generate gravitational radiation. For a recent account of the $`r`$-mode instability we refer to the review by Andersson and Kokkotas (2001).

Up to now, most of the studies of the $`r`$-mode instability have been done within the Newtonian framework. In a first attempt to study the relativistic analogue, Kojima (1997, 1998) considered the axial equations for slowly rotating stars in the so-called “low-frequency approximation”. He found a situation quite different from Newtonian theory, in that to first order the relativistic frame dragging in rotating stars gives rise to a continuous spectrum of frequencies instead of a single $`r`$-mode frequency. In a subsequent work, Beyer & Kokkotas (1999) proved in a mathematically rigorous way the existence of this continuous spectrum in the low-frequency approximation.

Lockitch, Andersson & Friedman (2001) showed that for uniform density stars, in addition to the continuous spectrum, Kojima’s equations also admits discrete mode solutions, the relativistic $`r`$ modes. Ruoff & Kokkotas (2001) and Yoshida (2001) extended this study to polytropic and realistic models and showed that for certain stellar models, no $`r`$ modes can be found, even when the low-frequency assumption is dropped and the radiation reaction is included (Ruoff & Kokkotas 2002; Yoshida & Futamase 2001). In all these studies, however, only the axial perturbations have been considered, and any coupling to the polar equations has been neglected. Lockitch et al. (2001) have argued that this assumption is only justified for non-barotropic stars, as the relativistic $`r`$ modes in barotropic stars should be hybrids, i.e. they retain a non-vanishing polar contribution in the non-rotating limit.

By including higher order rotational corrections and the coupling to the polar equations, Kojima & Hosonuma (2000) derived a fourth order equation, an extension of Kojima’s master equation, which as a consequence does not have the singular structure of the latter. Although not rigorously proved, it was argued that this could possibly resolve the problem of the non-existence of the $`r`$ modes in certain stellar models (Lockitch & Andersson 2001).

Recently Yoshida & Lee (2002) announced that they can find relativistic $`r`$ modes in the Cowling approximation for non-barotropic stellar models. They obtained their results by only including first order rotational correction terms. However, they did not restrict themselves to merely considering the coupling between the axial equation of order $`l`$ to the polar equations of order $`l+1`$, but instead they solved the complete system starting from $`l=|m|`$ up to some value $`l_{\mathrm{max}}`$, which they chose sufficiently high enough to ensure the convergence of the modes. Lockitch et al. (2001) had already formulated this problem in a somewhat different way, by neglecting only the perturbations which vanish in the stationary limit, such as the energy and pressure perturbations for instance, and some of the metric variables. So far they have presented numerical results only in a post Newtonian treatment of relativistic $`r`$ modes for uniform density stars. Their numerical eigenfunctions were used by Stergioulas & Font (2001) as approximate $`r`$-mode initial data for a 3-D general-relativistic hydrodynamical time evolution in the Cowling approximation. They found that the $`r`$ mode oscillates at a predominately discrete frequency with an indication of the kinematical drift, as predicted by Rezzolla, Lamb & Shapiro (2000).

In this paper, we follow the approach of Yoshida & Lee (2002) and extend it also to the time domain. Furthermore, we consider both barotropic and non-barotropic stellar models and discuss in detail how the continuous spectrum appears and how it is depends on the cut-off $`l_{\mathrm{max}}`$. By computing a sequence of stellar models with different $`M/R`$ ratios, we study how the continuous spectrum affects their inertial modes. In this paper, we call all the modes “inertial modes”, whose restoring force is mainly the Coriolis force. Thus we consider the $`r`$ modes as a subclass of inertial modes whose motion is dominated by the axial velocity field.

The paper is organized as follows: In Section 2 we present both the time dependent and the eigenvalue forms of the coupled perturbation equations. We first specialize on the axisymmetric case, where we describe how to compute the continuous spectrum and show its dependence on $`l_{\mathrm{max}}`$. The numerical methods and results are presented in Section 3 and we conclude with the discussion in Section 4. Throughout the paper we use geometrical units $`G=c=1`$. The over-dots and primes denote differentiation with respect to the time and the radial coordinates, respectively. All perturbation variables should carry indices $`l`$ and $`m`$, however, we conveniently omit them. Only in cases, where the $`l`$ index is important, we will include it. Finally, we assume $`m0`$, as one can always transform a $`m<0`$ mode with frequency $`\sigma `$ into a $`m>0`$ mode with frequency $`\sigma `$ by taking its complex conjugate.

## 2 FORMULATION OF THE PROBLEM

### 2.1 Background model

We consider a uniformly rotating relativistic star with angular velocity $`\mathrm{\Omega }`$, which is assumed to sufficiently small that any deviation from sphericity can be neglected. The metric can then be written as (Hartle 1967)

$$ds^2=e^{2\nu }dt^2+e^{2\lambda }dr^2+r^2\left(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2\right)2\omega \mathrm{sin}^2\theta dtd\varphi .$$
(1)
Assuming a perfect fluid star, the energy momentum tensor

$$T_{\mu \nu }=\left(p+ϵ\right)u_\mu u_\nu +pg_{\mu \nu },$$
(2)
with pressure $`p`$, energy density $`ϵ`$ and four-velocity components

$`u^t`$ $`=`$ $`e^\nu ,`$ (3)
$`u^\varphi `$ $`=`$ $`\mathrm{\Omega }e^\nu .`$ (4)
We will use a relativistic polytropic equation of state given by

$`p`$ $`=`$ $`\kappa \rho ^\mathrm{\Gamma },`$ (5)
$`ϵ`$ $`=`$ $`\rho +{\displaystyle \frac{p}{\mathrm{\Gamma }1}},`$ (6)
where $`\kappa `$ denotes the polytropic constant, $`\mathrm{\Gamma }`$ is the polytropic index and $`\rho `$ the rest mass density. With this form of the equation of state, the polytropic index $`\mathrm{\Gamma }`$ is equivalent to the adiabatic index and obeys the relation

$$\mathrm{\Gamma }=\frac{p+ϵ}{p}\frac{dp}{dϵ}.$$
(7)
### 2.2 Perturbed equations of motion

We assume the oscillations to be adiabatic, so that the relation between the Eulerian pressure perturbation $`\delta p`$ and energy density perturbation $`\delta ϵ`$ is given by

$$\delta p=\frac{\mathrm{\Gamma }_1p}{p+ϵ}\delta ϵ+p^{}\xi ^r\left(\frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}1\right),$$
(8)
where $`\mathrm{\Gamma }_1`$ represents the adiabatic index of the perturbed configuration and $`\xi ^r`$ is the radial component of the fluid displacement vector $`\xi ^\mu `$. The sound speed $`C_s`$ is

$`C_s^2`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}{\displaystyle \frac{dp}{dϵ}}.`$ (9)
The complete set of the perturbed Einstein equations has been derived using the BCL gauge (Battiston, Cazzola & Lucaroni 1971) in Ruoff, Stavridis & Kokkotas (2001), which we will refer to as RSK. As we are working in the Cowling approximation, we only need to consider the five fluid perturbations, which are the three components of the velocity perturbations $`\delta u_i`$, the enthalpy perturbation $`H`$ and the radial component of the displacement vector $`\xi ^\mu `$. Following RSK, we expand these quantities as

$`\delta u_r`$ $`=`$ $`e^\nu {\displaystyle \underset{l,m}{}}u_1^{lm}Y_{lm},`$ (10)
$`\delta u_\theta `$ $`=`$ $`e^\nu {\displaystyle \underset{l,m}{}}\left[\stackrel{~}{u}_2^{lm}_\theta Y_{lm}\stackrel{~}{u}_3^{lm}{\displaystyle \frac{_\varphi Y_{lm}}{\mathrm{sin}\theta }}\right],`$ (11)
$`\delta u_\varphi `$ $`=`$ $`e^\nu {\displaystyle \underset{l,m}{}}\left[\stackrel{~}{u}_2^{lm}_\varphi Y_{lm}+\stackrel{~}{u}_3^{lm}\mathrm{sin}\theta _\theta Y_{lm}\right],`$ (12)
$`\delta ϵ`$ $`=`$ $`{\displaystyle \underset{l,m}{}}{\displaystyle \frac{\left(p+ϵ\right)^2}{\mathrm{\Gamma }_1p}}\left(H^{lm}\xi ^{lm}\right)Y_{lm},`$ (13)
$`\xi ^r`$ $`=`$ $`\left[\nu ^{}\left(1{\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}\right)\right]^1{\displaystyle \underset{l,m}{}}\xi ^{lm}Y_{lm}.`$ (14)
Later we will introduce new variables $`u_2`$ and $`u_3`$ instead of $`\stackrel{~}{u}_2`$ and $`\stackrel{~}{u}_3`$. These include rotational corrections and therefore differ from the definitions of $`u_2`$ and $`u_3`$ in RSK. The time dependent perturbation equations follow directly from Eqs. (68)–(72) of RSK with all the metric perturbations set to zero:

$`\left(_t+\text{i}m\mathrm{\Omega }\right)H`$ $`=`$ $`C_s^2\left\{e^{2\nu 2\lambda }\left[u_1^{}+\left(2\nu ^{}\lambda ^{}+{\displaystyle \frac{2}{r}}\right)u_1e^{2\lambda }{\displaystyle \frac{\mathrm{\Lambda }}{r^2}}\stackrel{~}{u}_2\right]+\text{i}m\varpi H\right\}\nu ^{}e^{2\nu 2\lambda }u_1,`$ (15)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)u_1`$ $`=`$ $`H^{}+{\displaystyle \frac{p^{}}{\mathrm{\Gamma }_1p}}\left[\left({\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}1\right)H+\xi \right]\left[\omega ^{}+2\varpi \left(\nu ^{}{\displaystyle \frac{1}{r}}\right)\right]\left(\text{i}m\stackrel{~}{u}_2+_1^{\pm 1}\stackrel{~}{u}_3\right),`$ (16)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)\stackrel{~}{u}_2`$ $`=`$ $`H+{\displaystyle \frac{2\varpi }{\mathrm{\Lambda }}}\left(\text{i}m\stackrel{~}{u}_2+_3^{\pm 1}\stackrel{~}{u}_3\right){\displaystyle \frac{\text{i}mr^2}{\mathrm{\Lambda }}}A,`$ (17)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)\stackrel{~}{u}_3`$ $`=`$ $`{\displaystyle \frac{2\varpi }{\mathrm{\Lambda }}}\left(\text{i}m\stackrel{~}{u}_3_3^{\pm 1}\stackrel{~}{u}_2\right)+{\displaystyle \frac{r^2}{\mathrm{\Lambda }}}_2^{\pm 1}A,`$ (18)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)\xi `$ $`=`$ $`\nu ^{}\left({\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}1\right)e^{2\nu 2\lambda }u_1,`$ (19)
where

$`A`$ $`=`$ $`\varpi C_s^2e^{2\lambda }\left[u_1^{}+\left(2\nu ^{}\lambda ^{}+{\displaystyle \frac{2}{r}}\right)u_1e^{2\lambda }{\displaystyle \frac{\mathrm{\Lambda }}{r^2}}\stackrel{~}{u}_2\right]+\left[\varpi \left(\nu ^{}{\displaystyle \frac{2}{r}}\right)+\omega ^{}\right]e^{2\lambda }u_1`$ (20)
and

$`\mathrm{\Lambda }`$ $`=`$ $`l(l+1).`$ (21)
The operators $`_1^{\pm 1}`$, $`_2^{\pm 1}`$, and $`_3^{\pm 1}`$ are the same as in RSK and are defined by their action on a perturbation variable $`P^{lm}`$

$`_1^{\pm 1}P^{lm}`$ $`=`$ $`(l1)Q_{lm}P^{l1m}(l+2)Q_{l+1m}P^{l+1m},`$ (22)
$`_2^{\pm 1}P^{lm}`$ $`=`$ $`(l+1)Q_{lm}P^{l1m}+lQ_{l+1m}P^{l+1m},`$ (23)
$`_3^{\pm 1}P^{lm}`$ $`=`$ $`(l1)(l+1)Q_{lm}P^{l1m}+l(l+2)Q_{l+1m}P^{l+1m},`$ (24)
with

$`Q_{lm}`$ $`:=`$ $`\sqrt{{\displaystyle \frac{(lm)(l+m)}{(2l1)(2l+1)}}}.`$ (25)
By changing to new variables, we can simplify the above set of evolution equations. To this end we use Eq. (15) and rewrite Eq. (20) as

$`A`$ $`=`$ $`\varpi e^{2\nu }_tH+e^{2\lambda }Bu_1+O(\mathrm{\Omega }^2)`$ (26)
with

$`B`$ $`=`$ $`\omega ^{}+2\varpi \left(\nu ^{}{\displaystyle \frac{1}{r}}\right).`$ (27)
Neglecting second order terms and moving the terms including $`_tH`$ to the lefthand side, we can write Eqs. (17) and (18) as

$`_t\left(\stackrel{~}{u}_2+{\displaystyle \frac{\text{i}m\varpi r^2}{\mathrm{\Lambda }}}e^{2\nu }H\right)+\text{i}m\mathrm{\Omega }\stackrel{~}{u}_2`$ $`=`$ $`H+2{\displaystyle \frac{\varpi }{\mathrm{\Lambda }}}\left(\text{i}m\stackrel{~}{u}_2+_3^{\pm 1}\stackrel{~}{u}_3\right){\displaystyle \frac{\text{i}mr^2}{\mathrm{\Lambda }}}e^{2\lambda }Bu_1,`$ (28)
$`_t\left(\stackrel{~}{u}_3{\displaystyle \frac{\varpi r^2}{\mathrm{\Lambda }}}e^{2\nu }_2^{\pm 1}H\right)+\text{i}m\mathrm{\Omega }\stackrel{~}{u}_3`$ $`=`$ $`2{\displaystyle \frac{\varpi }{\mathrm{\Lambda }}}\left(\text{i}m\stackrel{~}{u}_3_3^{\pm 1}\stackrel{~}{u}_2\right)+{\displaystyle \frac{r^2}{\mathrm{\Lambda }}}e^{2\lambda }B_2^{\pm 1}u_1.`$ (29)
This suggests to define new variables

$`u_2`$ $`:=`$ $`\stackrel{~}{u}_2+\text{i}m{\displaystyle \frac{\varpi r^2}{\mathrm{\Lambda }}}e^{2\nu }H,`$ (30)
$`u_3`$ $`:=`$ $`\stackrel{~}{u}_3{\displaystyle \frac{\varpi r^2}{\mathrm{\Lambda }}}e^{2\nu }_2^{\pm 1}H,`$ (31)
and rewrite Eqs. (15)–(19) in terms of $`u_2`$ and $`u_3`$. After discarding any second order terms introduced by this replacement the perturbations equations now read

$`\left(_t+\text{i}m\mathrm{\Omega }\right)H`$ $`=`$ $`e^{2\nu 2\lambda }\left\{C_s^2\left[u_1^{}+\left(2\nu ^{}\lambda ^{}+{\displaystyle \frac{2}{r}}\right)u_1e^{2\lambda }{\displaystyle \frac{\mathrm{\Lambda }}{r^2}}u_2+2\text{i}m\varpi e^{2\lambda 2\nu }H\right]\nu ^{}u_1\right\},`$ (32)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)u_1`$ $`=`$ $`H^{}+{\displaystyle \frac{p^{}}{\mathrm{\Gamma }_1p}}\left[\left({\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}1\right)H+\xi \right]B\left(\text{i}mu_2+_1^{\pm 1}u_3\right),`$ (33)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)u_2`$ $`=`$ $`H+2{\displaystyle \frac{\varpi }{\mathrm{\Lambda }}}\left(\text{i}mu_2+_3^{\pm 1}u_3\right){\displaystyle \frac{\text{i}mr^2}{\mathrm{\Lambda }}}e^{2\lambda }Bu_1,`$ (34)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)u_3`$ $`=`$ $`2{\displaystyle \frac{\varpi }{\mathrm{\Lambda }}}\left(\text{i}mu_3_3^{\pm 1}u_2\right)+{\displaystyle \frac{r^2}{\mathrm{\Lambda }}}e^{2\lambda }B_2^{\pm 1}u_1,`$ (35)
$`\left(_t+\text{i}m\mathrm{\Omega }\right)\xi `$ $`=`$ $`\nu ^{}\left({\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}1\right)e^{2\nu 2\lambda }u_1.`$ (36)
In order to find the eigenmodes of this coupled system of equations, we assume our perturbations variables to have a harmonic time dependence $`\mathrm{exp}(\text{i}\sigma t)`$. By replacing all time derivatives by $`\text{i}\sigma `$ and letting $`H\text{i}H`$, $`u_3\text{i}u_3`$ and $`\xi \text{i}\xi `$ in order to obtain a purely real valued system of equations, we obtain two ODEs for $`H`$ and $`u_1`$ and three algebraic relations for $`u_2`$, $`u_3`$ and $`\xi `$. The relation for $`\xi `$, which follows from Eq. (36) is rather simple and can be immediately used to eliminate $`\xi `$, giving

$`H^{}`$ $`=`$ $`\left(\sigma +m\mathrm{\Omega }\right)u_1{\displaystyle \frac{p^{}}{\mathrm{\Gamma }_1p}}\left({\displaystyle \frac{\mathrm{\Gamma }_1}{\mathrm{\Gamma }}}1\right)\left[H\left(\sigma +m\mathrm{\Omega }\right)^1\nu ^{}e^{2\nu 2\lambda }u_1\right]+B\left(mu_2+_1^{\pm 1}u_3\right),`$ (37)
$`u_1^{}`$ $`=`$ $`\left(2\nu ^{}\lambda ^{}+{\displaystyle \frac{2}{r}}\right)u_1+2m\varpi e^{2\lambda 2\nu }H+e^{2\lambda }{\displaystyle \frac{\mathrm{\Lambda }}{r^2}}u_2C_s^2\left[\left(\sigma +m\mathrm{\Omega }\right)e^{2\lambda 2\nu }H\nu ^{}u_1\right],`$ (38)
$`u_2`$ $`=`$ $`\mathrm{\Sigma }^1\left[H{\displaystyle \frac{mr^2}{\mathrm{\Lambda }}}e^{2\lambda }Bu_1+{\displaystyle \frac{2\varpi }{\mathrm{\Lambda }}}_3^{\pm 1}u_3\right],`$ (39)
$`u_3`$ $`=`$ $`\mathrm{\Sigma }^1\left[{\displaystyle \frac{r^2}{\mathrm{\Lambda }}}e^{2\lambda }B_2^{\pm 1}u_1{\displaystyle \frac{2\varpi }{\mathrm{\Lambda }}}_3^{\pm 1}u_2\right],`$ (40)
where we have defined

$`\mathrm{\Sigma }:=\sigma +m\mathrm{\Omega }{\displaystyle \frac{2m\varpi }{\mathrm{\Lambda }}}.`$ (41)
The relations (39) and (40) cannot be solved directly for $`u_2`$ and $`u_3`$ as they involve the coupling operators $`_i^{\pm 1}`$. We should stress again that because of this coupling, both sets of equations (32)–(36) and (37)–(40) have to be considered as infinite systems with $`l`$ running from $`m`$ to infinity. There is no coupling to equations with $`l<m`$, since from its definition (25) it follows that the coupling coefficient $`Q_{mm}`$, which would couple to $`l=m1`$, is zero. Furthermore, the equations form two independent sets, each belonging to a different parity. In the non-rotating case, one usually distinguishes between polar and axial perturbations. Under parity transformation polar perturbations change sign as $`(1)^l`$, axial perturbations as $`(1)^{l+1}`$. In the rotating case, however, the polar and axial equations are coupled, and therefore the distinction between polar and axial modes cannot be upheld any more. However, the equations do not mix the overall parity as can be seen as follows. A polar equation with even $`l=m`$ has even parity. It is coupled to an axial equation with $`l+1`$, whose parity is also even. The next coupling is again to a polar equation with $`l+2`$, thus having even parity. As this continues in the same manner, this means that for even $`m`$, the complete coupled system with a leading polar equation has even parity. Conversely, the other system starting with a leading axial equation has odd parity. For odd $`m`$, we obtain the reversed situation. This implies that for any given $`m`$, it makes sense to distinguish the modes according to their overall parity. Lockitch & Friedman (2000) introduced the notion of polar or axial led modes, depending on whether the leading equations with $`l=m`$ are axial or polar. For even $`m`$, the polar led modes have even parity and the axial modes odd parity, and vice versa for odd $`m`$. The axisymmetric case $`m=0`$ is somewhat special in that only the even parity modes start with $`l=0`$ whereas the odd parity modes have to start with $`l=1`$, as there are no axial $`l=0`$ equations. In the following section, we will restrict ourselves to positive $`m`$.

### 2.3 The Continuous Spectrum

For simplicity, we now focus on axisymmetric perturbations with $`m=0`$, which leads to considerable simplifications of the time independent equations (37)–(40). With the expansion of the operators $`_i^{\pm 1}`$ according their definitions (22) – (24), the $`m=0`$ equations read (we now explicitly include the index $`l`$):

$`(H^l)^{}`$ $`=`$ $`\sigma u_1^l+A\left(l(l1)Q_lu_3^{l1}(l+1)(l+2)Q_{l+1}u_3^{l+1}\right),`$ (42)
$`(u_1^l)^{}`$ $`=`$ $`\left(2\nu ^{}\lambda ^{}+{\displaystyle \frac{2}{r}}\right)u_1^l+e^{2\lambda }{\displaystyle \frac{\mathrm{\Lambda }}{r^2}}u_2^lC_s^2\left[\sigma e^{2\lambda 2\nu }H^l\nu ^{}u_1^l\right],`$ (43)
$`u_2^l`$ $`=`$ $`\sigma ^1\left[H^l+2\varpi \left((l1)Q_lu_3^{l1}+(l+2)Q_{l+1}u_3^{l+1}\right)\right],`$ (44)
$`u_3^l`$ $`=`$ $`\sigma ^1\left[r^2e^{2\lambda }A\left(Q_lu_1^{l1}Q_{l+1}u_1^{l+1}\right)+2\varpi \left((l1)Q_lu_2^{l1}+(l+2)Q_{l+1}u_2^{l+1}\right)\right],`$ (45)
with

$`Q_l`$ $`:=`$ $`\left(4l^21\right)^{1/2}.`$ (46)
For $`l=0`$, we have $`u_2^0=u_3^0=0`$. To solve the ODEs (42) and (43), we have to compute $`u_2^l`$ and $`u_3^l`$ in terms of $`H^l`$ and $`u_1^l`$. This can be done the most easily by combining $`u_2^l`$ and $`u_3^l`$ as well as $`H^l`$ and $`u_1^l`$ into vectors and rewriting the two algebraic relations (44) and (45) as

$`\sigma \left(\begin{array}{c}u_2^l\\ u_3^l\end{array}\right)`$ $`=`$ $`2\varpi (l1)Q_l\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}u_2^{l1}\\ u_3^{l1}\end{array}\right)+2\varpi (l+2)Q_{l+1}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)\left(\begin{array}{c}u_2^{l+1}\\ u_3^{l+1}\end{array}\right)`$ (70)
$`+r^2e^{2\lambda }AQ_l\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}H^{l1}\\ u_1^{l1}\end{array}\right)+\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right)\left(\begin{array}{c}H^l\\ u_1^l\end{array}\right)+r^2e^{2\lambda }AQ_{l+1}\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right)\left(\begin{array}{c}H^{l+1}\\ u_1^{l+1}\end{array}\right).`$
Defining

$`u^l`$ $`=`$ $`\left(\begin{array}{c}u_2^l\\ u_3^l\end{array}\right),s^l=\left(\begin{array}{c}H^l\\ u_1^l\end{array}\right),`$ (75)
$`U_{}^l`$ $`=`$ $`2\varpi (l1)Q_l\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),U_\sigma =\sigma \left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),U_+^l=2\varpi (l+2)Q_{l+1}\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),`$ (82)
$`S_{}^l`$ $`=`$ $`r^2e^{2\lambda }AQ_l\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right),S=\left(\begin{array}{cc}1& 0\\ 0& 0\end{array}\right),S_+^l=r^2e^{2\lambda }AQ_{l+1}\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right),`$ (89)
we can cast Eq. (70) into a more compact form

$`U_{}^lu^{l1}+U_\sigma u^l+U_+^lu^{l+1}`$ $`=`$ $`S_{}^ls^{l1}+Ss^l+S_+^ls^{l+1},l=1,\mathrm{},\mathrm{}`$ (90)
This, again, can be viewed as a matrix equation for two infinitely dimensional matrices U and S acting on the vectors

$`u`$ $`=`$ $`(u^1,u^2,\mathrm{},u^l,\mathrm{})^\mathrm{T},`$ (91)
$`s`$ $`=`$ $`(s^1,s^2,\mathrm{},s^l,\mathrm{})^\mathrm{T},`$ (92)
whose respective elements are the 2-vectors $`u^l`$ and $`s^l`$. Both U and S are tridiagonal block matrices, with each block given by the above $`2\times 2`$ matrices. Explicitly, U is given by

$`𝖴`$ $`=`$ $`\left(\begin{array}{cccccccc}U_\sigma & U_+^1& 0& \mathrm{}& & & & \\ U_{}^2& U_\sigma & U_+^2& 0& \mathrm{}& & & \\ 0& U_{}^3& U_\sigma & U_+^3& 0& \mathrm{}& & \\ & & \mathrm{}& \mathrm{}& \mathrm{}& & & \\ & \mathrm{}& 0& U_{}^l& U_\sigma & U_+^l& 0& \mathrm{}\\ & & & & \mathrm{}& \mathrm{}& \mathrm{}& \end{array}\right).`$ (99)
A similar structure holds for S. Now we can write Eq. (90) as

$`{\displaystyle \underset{l^{}}{}}\text{U}^{ll^{}}u^l^{}`$ $`=`$ $`{\displaystyle \underset{l^{}}{}}\text{S}^{ll^{}}s^l^{},l=1\mathrm{}\mathrm{},`$ (100)
which can be solved for $`u^l`$ by multiplying both sides with $`𝖴^1`$

$`u^l`$ $`=`$ $`{\displaystyle \underset{l^{}l^{\prime \prime }}{}}\left(\text{U}^1\right)^{ll^{}}\text{S}^{l^{}l^{\prime \prime }}s^{l^{\prime \prime }},l=1\mathrm{}\mathrm{}.`$ (101)
It is important to note that the matrix U is $`r`$-dependent since its elements $`U_\pm ^l`$ contain the function $`\varpi `$. This raises the question whether one can assure that U will be invertible for any frequency $`\sigma `$ at each value of $`r`$ in the stellar interior. For U to be invertible its determinant must not vanish.

In the Newtonian limit $`\varpi \mathrm{\Omega }`$, the matrix U does not depend on $`r`$ any more, and one can easily show that U becomes singular only for a discrete set of values of $`\sigma `$. These frequencies represent the solutions of the homogeneous part $`𝖴u=0`$ of Eq. (100). In the relativistic case, however, we have to replace $`\mathrm{\Omega }`$ by the $`r`$-dependent effective angular velocity $`\varpi `$. Thus the zeroes of $`det𝖴`$ now depend on the position inside the star. This means that each single Newtonian frequency will be spread out in the relativistic case into a continuous band of frequencies, determined by the values of $`\varpi `$ at the centre and at the surface of the star, which we will denote by $`\varpi _0`$ and $`\varpi _R`$, respectively. The total range of the continuous spectrum is then the sum of all the individual bands originating from their Newtonian discrete values.

In order to numerically perform the inversion and solve equation (101), we have to truncate the system at some value $`l_{\mathrm{max}}>1`$. Let us consider the simple case, where we truncate the system at $`l_{\mathrm{max}}=2`$, i.e. we couple only $`l=1`$ with $`l=2`$ and neglect all $`l3`$. Then the matrix U is given by

$`𝖴`$ $`=`$ $`\left(\begin{array}{cc}U_\sigma & U_+^1\\ U_{}^2& U_\sigma \end{array}\right)=\left(\begin{array}{cccc}\sigma & 0& 0& 6\varpi Q_2\\ 0& \sigma & 6\varpi Q_2& 0\\ 0& 2\varpi Q_2& \sigma & 0\\ 2\varpi Q_2& 0& 0& \sigma \end{array}\right)`$ (108)
with $`Q_2=1/\sqrt{15}`$ and $`det𝖴=\mathrm{𝟢}`$ yields

$`\sigma ^2`$ $`=`$ $`{\displaystyle \frac{4}{5}}\varpi ^2=\mathrm{\hspace{0.33em}0.8}\varpi ^2.`$ (109)
Hence, the range of the continuous spectrum is given by

$`0.894\varpi _0\sigma _c0.894\varpi _R,`$ (110)
if we restrict ourselves to positive frequencies. If we also include the $`l=3`$ terms, we obtain

$`\sigma ^2`$ $`=`$ $`{\displaystyle \frac{84}{49}}\varpi ^2\mathrm{\hspace{0.33em}1.714}\varpi ^2.`$ (111)
and the range is shifted to

$`1.309\varpi _0\sigma _c1.309\varpi _R.`$ (112)
Going one step further and including $`l=4`$, we obtain two solutions

$`\sigma ^2`$ $`=`$ $`{\displaystyle \frac{4}{441}}\left(147\pm 42\sqrt{7}\right)\varpi ^2,`$ (113)
corresponding to the frequency intervals

$`1.530\varpi _0\sigma _c1.530\varpi _R`$ (114)
for the plus sign and

$`0.5704\varpi _0\sigma _c0.5704\varpi _R`$ (115)
for the minus sign. Adding $`l=5`$ will still yield two ranges and for $`l=6`$ we obtain three. The pattern is as follows: When a coupling to a higher even value of $`l`$ is included, one additional range of the continuous spectrum will appear. Depending on the respective values of $`\varpi _0`$ and $`\varpi _R`$, the various ranges might as well overlap. This actually happens for the more relativistic stars since the variation in $`\varpi `$ is much greater.

A similar picture also holds in the case $`m0`$. Here, each coupling to a higher value of $`l`$ results in one additional band of the continuous spectrum. We have already stated that in the case of the matrix U becoming singular at some point inside the star, we cannot invert it anymore and the integration of the eigenvalue system (37)–(40) fails. In the purely axial case, Ruoff & Kokkotas (2001) showed that by a series expansion one could still obtain mathematically valid mode solutions inside the continuous spectrum, however, only at the price of a having a divergence in the associated fluid perturbation at the singular point. Thus they concluded that these modes are unphysical and should be discarded. The question is whether we can adopt this point of view and claim that inside the continuous bands modes cannot exist.

It is clear that the widths of the continuous bands strongly depend on the stellar parameters and go to zero in the Newtonian limit. For weakly relativistic stellar models, each individual band is fairly small, but the total range, which is the sum of these small bands, can become quite large. If it is true that modes should not exist inside the continuous bands, then even for weakly relativistic stellar models, there should be almost no true mode solutions. As we cannot perform the integration of (37)–(40) inside the continuous bands, we cannot tell whether there still might be valid mode solutions. Trying series expansion of the coupled equations around the singular point will be quite involved. This is why we rather rely on the evolution of the time dependent equations, which are free of singularities. By examining the Fourier spectra we can assess whether or not we can find modes inside the continuous spectrum. As we shall show and discuss, there are indeed modes which lie inside the continuous bands. However, as the stellar models become more compact, the continuous spectrum has the ability to destroy more and more of the inertial modes.

## 3 Numerical methods and results

The numerical evolution of the coupled equations does not present any severe numerical problems. Although one might encounter long term instabilities, if a large number of $`l`$ is included, these can usually be overcome by increasing the spatial resolution or shifting the numerical origin some grid points away from zero. Up to values of about $`l_{\mathrm{max}}=10`$, evolutions with $`200`$ grid points inside the star yield quite accurate results for the inertial modes. As the inertial modes are located at at the lower end of the frequency spectrum, we have to perform long time evolutions in order to obtain a sufficient frequency resolution. Usually, we have to evolve up to the order of $`500`$ms, corresponding to a total number of time steps well above $`10^6`$.

The numerical mode calculation is somewhat more troublesome. As explained in the previous section, we cannot compute the modes once they are inside the continuous spectrum, because the matrix U becomes singular and inversion is no longer possible. To find the modes outside the continuous spectrum, we proceed as follows. We first rescale the variables according to

$`H^l`$ $``$ $`r^lH^l,`$ (116)
$`u_1^l`$ $``$ $`r^{l1}u_1^l,`$ (117)
$`u_2^l`$ $``$ $`r^lu_2^l,`$ (118)
$`u_3^l`$ $``$ $`r^{l+1}u_3^l.`$ (119)
This ensures that all our new variables have Taylor expansions around the origin starting with the constant term. To initiate the integration from the centre toward the stellar surface, we have to prescribe initial data for each $`H^l`$ and $`u_1^l`$. At each integration step, we have to invert the matrix U in order to compute $`u_2^l`$ and $`u_3^l`$. At the origin, $`H^l`$ and $`u_1^l`$ are not independent, but are related to leading order by

$`H^l(0)`$ $`=`$ $`\left(\sigma +m\mathrm{\Omega }{\displaystyle \frac{2m\varpi _0}{l}}4\varpi _0^2\mathrm{\Sigma }_{l1}^1{\displaystyle \frac{(2l+1)(l1)}{l^2}}Q_{lm}^2\right){\displaystyle \frac{u_1^l(0)}{l}}.`$ (120)
To leading order, however, there is no relation between $`H^l`$ and any other $`H^l^{}`$, so that we have $`n=l_{\mathrm{max}}m+1`$ degrees of freedom for the initial values, as $`l`$ runs from $`m`$ to $`l_{\mathrm{max}}`$. The frequency $`\sigma `$ represents a further degree of freedom, which we have to choose deliberately. The missing boundary conditions, which fix these degrees of freedom, come from the requirement of a vanishing Lagrangian pressure perturbation at the stellar surface. Since this requirement translates into the simultaneous vanishing of the quantities

$`\mathrm{\Delta }p^l`$ $`:=`$ $`\left[\left(\sigma +m\mathrm{\Omega }\right)e^{2\lambda 2\nu }H^l\nu ^{}u_1^l\right]_{r=R},mll_{\mathrm{max}},`$ (121)
we obtain $`l_{\mathrm{max}}m+1`$ conditions. We have thus reduced the number of degrees of freedom to one, which is the value one expects for linear systems, as the overall amplitude of the perturbations remains arbitrary. Let us denote a set of initial data for the $`H^l`$ by $`_k:=(H_k^m,\mathrm{},H_k^{l_{\mathrm{max}}})`$. In order to compute the modes, we choose a frequency $`\sigma `$ and perform $`n`$ integrations for the following $`n`$ different sets of initial data: $`_1=(1,0,\mathrm{},0),_2=(0,1,0,\mathrm{},0),\mathrm{},_n=(0,\mathrm{},1)`$. The respective initial data for $`u_1`$ follow from (120). For each set $`_k`$, we compute the $`n`$ Langragian pressure perturbations $`\mathrm{\Delta }p_k^l`$, which in general will not be zero. After having performed the $`n`$ integrations, we end up with $`n\times n`$ values for the Langrangian pressure perturbations $`\mathrm{\Delta }p_k^l`$, which can be arranged in a $`n\times n`$ “pressure matrix” P. As the $`n`$ solutions pertaining to the $`n`$ sets of initial data $`_k`$ are linearly independent, we can try to use them to construct another solution by a suitable linear combination, which satisfies $`\mathrm{\Delta }p^l=0`$ for all $`l`$ simultaneously. In other words we have to find $`n`$ coefficients $`a_k`$, such that the linear combinations $`a_k\mathrm{\Delta }p_k^l`$ vanish for all $`l`$. This is possible if and only if the determinant of the pressure matrix P vanishes. The algorithm for finding the eigenmode frequency is then obvious. For each $`\sigma `$, we construct the pressure matrix P and compute its determinante. If zero, we have found an eigenmode.

To also compute the eigenfunction, we have to solve the homogeneous system of equations

$`{\displaystyle \underset{k=m}{\overset{l_{\mathrm{max}}}{}}}a_k\mathrm{\Delta }p_k^l=0,l=m\mathrm{}l_{\mathrm{max}}`$ (122)
for the coefficients $`a_k`$, which can be easily done by a “singular value decomposition” of P (see e.g. Press et al. 1992). The eigenfunction can be obtained by the linear combination of the $`n`$ solutions corresponding to the $`n`$ sets of initial data $`_k`$, with the appropriate weighing coefficients $`a_k`$.

Although this method is quite easy to implement, it suffers from the drawback of breaking down when $`l_{\mathrm{max}}m+7`$. This failure results from the insufficient prescription of the initial data, where we should include higher order terms in the Taylor expansion, as the leading orders are set to zero for all $`l`$ but one. For the purpose of computing the inertial modes for small values of $`l`$, however, our method is robust enough, and we obtain perfect agreement with the frequencies obtained from the time evolutions. Different methods are used by Lockitch et al. (2001) and Yoshida & Lee (2002).

We have checked the consistency of our codes by feeding the evolution code with the eigenfunctions obtained from the eigenvalue code. As expected the time evolution yields a standing wave with one single frequency. Arbitrary initial data excite the $`f`$ and $`p`$ modes together with the inertial and in the non-barotropic case also the $`g`$ modes. The frequencies found by Fourier transforming the time sequences agree with those found from the mode calculation within less than one per cent. This difference depends only on the numerical resolution of the time evolutions and tends to zero for increasing resolution.

### 3.1 The axisymmetric case

Let us first assess the influence of the rotation on the $`f`$ and $`p`$ modes. We take the same stellar model as Font et al. (2001), which is a $`\mathrm{\Gamma }=2`$ polytrope with $`\kappa =217.86`$km<sup>2</sup>. (Note their value $`\kappa =100`$ is in dimensionless units where $`G=c=M_{}=1`$). With a central rest mass density $`\rho _0=7.914\times 10^{14}`$g/cm<sup>3</sup>, the total mass is $`M=1.4M_{}`$ and the radius is $`R=14.15`$km. The rotation rate is given by $`\mathrm{\Omega }=4230`$s<sup>-1</sup> corresponding to $`\epsilon =\mathrm{\Omega }/\mathrm{\Omega }_K=0.5223`$, where $`\mathrm{\Omega }_K=\sqrt{M/R^3}`$. In the following we consider only barotropic perturbations, i.e. we use $`\mathrm{\Gamma }_1=2`$, thus suppressing the $`g`$ modes. We will study non-barotropic perturbations in the non-axisymmetric case.

In Table 1, we show the change in frequencies of the (quasi)radial $`f^0`$ and $`p^0`$ modes, when the coupling to higher values of $`l`$ is included. As we are considering axisymmetric perturbations, there is no $`m`$-splitting and rotational effects only come from the coupling of the equations with different $`l`$. For $`l_{\mathrm{max}}=0`$, there are no rotational corrections, i.e. the equations describe the non-rotating case. For $`l_{\mathrm{max}}=1`$, the polar $`l=0`$ equations couple to the axial $`l=1`$ equation for $`u_3`$, and we find a shift in the frequencies towards higher values. This shift is about 3.7 per cent for the $`f^0`$ mode, 1.3 per cent for the $`p_1^0`$ mode and even less for higher modes. When further couplings to higher $`l`$ are included, the additional frequency corrections become smaller and the frequencies converge rapidly. Beyond $`l_{\mathrm{max}}=3`$, the frequency corrections are less than 0.01 per cent. Hence, we see that the quasi-radial $`f`$ and $`p`$ modes are quite insensitive to rotation. The effects of rotation on the axisymmetric modes have been investigated using the fully nonlinear equations in the Cowling approximation by Font et al. (2001). However, a comparison between our results and theirs is not very meaningful as they study rapidly rotating models, which undergo quite large radius and mass corrections with respect to the non-rotating case. In the perturbation formalism, these corrections appear only at the $`\mathrm{\Omega }^2`$ level, so that we cannot account for them when only first order corrections are included. That these corrections are important can be seen from the fact that their frequencies decrease with increasing rotation rate while ours actually increase. By including some of the second order corrections we have convinced ourselves that the frequencies get pushed down to smaller values, which then agree much better with the ones of Font et al. (2001). Therefore, we compare our values in Table 1 with theirs only in the non-rotating limit. For the first four radial modes, they find the frequencies $`2708`$Hz, $`4547`$Hz, $`6320`$Hz and $`8153`$Hz, which agree very well with the first row of Table 1. For $`m0`$, we find the well-known linear splitting of the $`f`$ and $`p`$ modes into a lower and a higher frequency part.

Let us now turn our attention towards the inertial modes. In the axisymmetric case, the first inertial modes appear, when coupling to $`l=2`$ is included. As described in the previous section, a continuous band of frequencies appears at the same time. In the time evolution, the modes can be distinguished from the continuous band as follows. Place some observers at different locations inside the star. Take the FFT of each observer’s time series and plot the various spectra on top of each other. Those peaks that coincide for all observers belong to modes, while those peaks which are different for each observer correspond to frequencies from the continuous spectrum.

As discussed above, the location of the continuous spectrum depends strongly on how many $`l`$ are taken into account. Depending on $`l_{\mathrm{max}}`$, it can consist of several patches. The width of each patch is determined by the compactness of the stellar models and shrinks to a set of discrete values in the Newtonian limit. For very compact stars, the patches can become quite broad and in general overlap. In the following, we will investigate the behaviour of both the continuous spectrum and the inertial modes as functions of the stellar compactness $`M/R`$ and the value of $`l_{\mathrm{max}}`$. As a first step, however, it is helpful to understand the behaviour of the inertial modes for an almost Newtonian stellar model, where the ranges of the continuous spectrum are very narrow.

Lindblom & Ipser (1999) have shown that for the Newtonian Maclaurin spheroids, the eigenvalue problem can be made separable by choosing appropriate spheroidal coordinates. Hence, $`l`$ and $`m`$ are “good” quantum numbers, which can be used to label the modes. For each value of $`m>0`$ and $`l>m`$ there exist $`lm`$ inertial modes. For $`m=0`$, one finds $`l1`$ mode solutions. In this case, however, the eigenvalue problem can be written in terms of $`\sigma ^2`$, hence both $`\sigma `$ and $`\sigma `$ are valid eigenfrequencies. For even $`l`$, there is an odd number of eigenvalues, hence one eigenvalue has to be $`\sigma =0`$ (c.f. Table 1 of Lindblom & Ipser 1999). If one restricts oneself to positive frequencies, there is no mode for $`l=2`$, one mode for $`l=3`$ and 4, two modes for $`l=5`$ and 6 and so on.

In our case, $`l`$ is not a good quantum number any more, but we still can assign a certain value $`l`$ to each mode as we shall demonstrate. To this end we now study for an almost Newtonian stellar model with compactness $`M/R=0.01`$ how both the continuous spectrum and the inertial modes are affected as $`l_{\mathrm{max}}`$ is successively incremented.

Figure 1 shows the locations of the continuous bands (grey shaded areas) together with the rotational modes one can find for different values of $`l_{\mathrm{max}}`$, ranging from $`l_{\mathrm{max}}=2`$ to 5. For $`l_{\mathrm{max}}=2`$, we have a single continuous band in whose vicinity many modes can be found. These become more densely spaced as they approach the continuous spectrum. The modes to the left of the continuous spectrum have even parity and the modes to the right have odd parity. When we go to $`l_{\mathrm{max}}=3`$, we obtain a similar picture. The continuous spectrum has moved to the right and again is surrounded by a large number of modes. In addition, we find a more isolated mode closer to the position, where the continuous spectrum for $`l_{\mathrm{max}}=2`$ was located. This is actually the first “true” inertial mode, as it does not change if we further increment $`l_{\mathrm{max}}`$. All the other modes which we can find for $`l_{\mathrm{max}}=2`$ and $`l_{\mathrm{max}}=3`$ are “fake” modes, i.e. they are artefacts due to the truncation of the equations and dissappear again when the coupling to higher $`l`$ is included. We will label the “true” inertial modes with $`i_n^{l_0}`$, where $`l_0`$ is the value of $`l_{\mathrm{max}}`$ for which it appears for the first time and does not change its frequency significantly when $`l_{\mathrm{max}}`$ is further increased. Since for a given $`l_{\mathrm{max}}`$, more than just one new inertial mode can appear, we order these modes through the lower index $`n`$ with ascending frequency, i.e. $`\sigma _n^{l_0}<\sigma _{n+1}^{l_0}`$. According to this convention, the inertial mode appearing for $`l_{\mathrm{max}}=3`$ is the $`i_1^3`$ mode, which has odd parity. Our upper index $`l_0`$ agrees with the definition of Yoshida & Lee (2000a) and differs from the one of Lockitch & Friedman, who use $`l_01`$ to denote the same mode. Yoshida & Lee (2000a), however, use the lower index $`n`$ in a somewhat different form, as they compute the modes in the rotating frame, which then form a set with equal fractions of negative and positive frequency modes. Yoshida & Lee attach negative indices $`n`$ to the negative frequencies modes and positive $`n`$ to the positive frequency modes.

For $`l_{\mathrm{max}}=4`$, the continuous spectrum consists of two patches, and we find the $`i_1^4`$ mode, which now has even parity. For $`l_{\mathrm{max}}=5`$, these two patches shift to the right and two new inertial modes, $`i_1^5`$ and $`i_2^5`$ appear. The emerging pattern is as follows: For each even $`l_{\mathrm{max}}`$, a new patch of the continuous spectrum appears. If $`l_{\mathrm{max}}>2`$ is odd, $`(l_{\mathrm{max}}1)/2`$ new odd parity inertial modes appear, if $`l_{\mathrm{max}}`$ is even, we find $`(l_{\mathrm{max}}2)/2`$ new even parity modes. Besides the inertial modes for a given $`l_{\mathrm{max}}`$ appear close to the positions where the patches of the continuous spectrum were located for $`l_{\mathrm{max}}1`$. Examination of the eigenfunctions reveals that a mode $`i_n^l`$ has no node in the coefficient $`H^l`$, one node in the coefficient $`H^{l2}`$, and $`n`$ nodes in the coefficient $`H^{l2n}`$.

It can be noticed that for $`l_{\mathrm{max}}=4`$ and 5, there are only very few modes in the vicinity of the left patches of the continuous spectrum, which is due to the fact that these patches are somewhat broader than the other ones. If one considers less relativistic models, the width of the patches shrink and more modes can be found around them. In the Newtonian case, where all the patches degenerate to single points, we expect an infinite number of modes around these points. As the star becomes more relativistic and the continuous bands expand, they will swallow more and more modes. In principle this is of no relevance as these modes are not physical, anyway. The remaining question is, whether or not the continuous spectra are able to swallow the “true” inertial modes, if they grow broad enough.

The answer is depicted in Fig. 2, where we plot the inertial modes together with the continuous spectra as functions of both the compactness $`M/R`$ and $`l_{\mathrm{max}}`$. Here, we plot only the “true” inertial modes and not those resulting from the truncation of the equations at $`l_{\mathrm{max}}`$ as in Fig. 1. The four panels illustrate the situation for the different cut-off values $`l_{\mathrm{max}}`$ ranging from 3 to 6. The compactness $`M/R`$ is from 0.01 (as in Fig. 1) to 0.22. The stability limit with respect to radial collapse is reached at $`M/R=0.215`$. All stellar models have the same $`\epsilon =0.5223`$, which means that the less compact models rotate much more slowly than the highly relativistic ones. Since the frequencies of the inertial modes depend essentially linearly on the rotation rate, they would become very small for less compact models. This is the reason why we express them in units of $`\mathrm{\Omega }_K`$. We thus see that the inertial modes lie in the interval between 0 and 1 for the whole compactness range. (Actually they lie in the interval between -1 and 1, but as negative and positive frequencies have the same magnitude, we plot only the positive ones).

Starting with $`l_{\mathrm{max}}=3`$, we have one single continuous spectrum, which broadens considerably as the compactness $`M/R`$ is increased. For small $`M/R`$, we also find the $`i_1^3`$ mode, which moves to lower frequencies as $`M/R`$ increases. However, the lower border of the continuous spectrum decreases much faster and for $`M/R0.17`$ it reaches the mode. For larger $`M/R`$ the mode cannot be found any more. We stress again that once the mode would be inside the continuous spectrum, the eigenvalue code breaks down because the matrix $`𝖴`$ becomes singular and inversion is no longer possible. In this case, we take the results of the time evolutions to decide whether or not the mode exists inside the continuous spectrum. And indeed, the time evolution shows that for $`M/R>0.17`$ the mode does not penetrate the continuous spectrum but instead disappears.

When we go to $`l_{\mathrm{max}}=4`$, the same happens to the $`i_1^4`$ mode which reaches the upper patch of the continuous spectrum at about $`M/R0.11`$ and disappears for larger $`M/R`$. Surprisingly, the $`i_1^3`$ mode now exists over the complete compactness range. Indeed, the time evolution shows that it can even exist inside the continuous spectrum for $`M/R>0.2`$. This feature is even more apparent when the coupling up to $`l_{\mathrm{max}}=5`$ is included. In this case, the $`i_1^3`$ mode is always inside the continuous spectrum over the whole compactness range. However, the $`i_1^4`$ mode still does not penetrate the continuous spectrum and the same is true for the $`i_1^5`$ and $`i_2^5`$ modes. Notice also that for $`M/R>0.17`$, the two patches of the continuous spectrum begin to overlap.

For $`l_{\mathrm{max}}=6`$, both the $`i_1^3`$ and $`i_1^3`$ modes exist for the complete compactness range. Again, the $`i_1^5`$ is always inside the continuous spectrum, whereas the $`i_1^5`$ mode enters it only for $`M/R>0.2`$. The $`i_2^5`$ mode, in contrast, ceases to to exist for $`M/R>0.1`$. For the $`l_{\mathrm{max}}=7`$ case, which is not shown, only the $`i_1^3`$, $`i_1^5`$ and $`i_1^7`$ modes exist for the complete compactness range, most of the others exist only outside the continuous spectrum. Some actually do penetrate the continuous spectrum, but if the star becomes too compact, they eventually vanish.

The emerging picture is as follows. Only the modes with low frequencies, such as the $`i_1^3`$, $`i_1^5`$ and $`i_1^7`$ modes for instance are able to exist for the complete compactness range and even inside the continuous spectrum. Those with higher frequencies exist only for a much smaller range. This might be attributed to the fact that the patches of the continuous spectrum are much broader for larger frequencies. As a consequence, the total number of inertial modes decreases as the stellar models become more relativistic. However, this leaves the question why some of the modes can survive inside the continuous spectrum and some cannot. We conjecture that this has to do with the value of $`l_{\mathrm{max}}`$, where the mode occurs for the first time. If for a certain value of $`l_{\mathrm{max}}`$ it appears well outside the continuous spectrum, then it will survive if $`l_{\mathrm{max}}`$ is further increased and the mode suddenly finds itself inside the continuous spectrum. In other words, it is not the mode which moves inside the continuous spectrum as $`l_{\mathrm{max}}`$ is increased, it is the continuous spectrum which shifts itself to the position where the mode is located. If, however, for a certain $`l_{\mathrm{max}}`$ the continuous spectrum is already sufficiently broad to occupy the position where an inertial mode would appear, then this does not seem to be able to establish itself, even when higher values of $`l`$ are included.

### 3.2 The non-axisymmetric case

For $`m0`$, we find essentially the same picture as in the axisymmetric case, i.e. some of the modes can exist for the whole compactness range while some cannot. Here, the most interesting question, of course, is what happens to the relativistic $`r`$-modes and in which way are they affected by the continuous spectrum. One can see that the relativistic case differs qualitatively from the Newtonian one, as in the latter, one can immediately show that the leading order of the $`r`$-mode frequency (in the inertial frame) is given by

$`\sigma `$ $`=`$ $`m\mathrm{\Omega }\left(1{\displaystyle \frac{2}{l(l+1)}}\right),`$ (123)
whereas the analogous relativistic calculation yields

$`\sigma (r)`$ $`=`$ $`m\mathrm{\Omega }\left[1{\displaystyle \frac{2}{l(l+1)}}\left(1{\displaystyle \frac{\omega (r)}{\mathrm{\Omega }}}\right)\right].`$ (124)
Instead of a single frequency, we have a continuous spectrum of frequencies, determined by the values of the frame dragging $`\omega `$ at the centre and surface of the star. Herein, we have chosen the sign of $`\sigma `$ such that it yields negative frequencies for positive values of $`m`$.

The above values for $`\sigma `$ were obtained by assuming that the $`r`$ mode is purely axial to leading order, and the coupling to the polar equations is just a higher order correction. However, it was shown by Lockitch et al. (2001) that this is not true for barotropic perturbations, as this assumption leads to some incompatibilities within the polar equations. Instead, the barotropic $`r`$ mode necessarily has to include polar contributions, which do not vanish in the non-rotating limit. For a “pure” $`r`$ mode, the polar contribution has to vanish. In the terminology of Lockitch et al. (2001), the relativistic barotropic $`r`$ mode is therefore an “axial-led hybrid” mode. In the non-barotropic case, things are different, as the $`r`$ mode cannot have non-vanishing polar contributions in the non-rotating limit. This is because for non-rotating non-barotropic stellar models, there are no stationary polar perturbations, since they would immediately excite $`g`$ modes. For non-rotating barotropic stars on the other hand, the $`g`$ modes form a degenerate space of zero-frequency modes, which is why the axial inertial modes can have polar contributions. Thus, for non-barotropic stars, the $`r`$ mode should not be an axial-led hybrid, but it should reduce in the non-rotating limit to a purely axial stationary perturbation. Nevertheless, we will call the inertial mode $`i_1^l`$ with $`l=m+1`$ an $`r`$ mode, regardless whether the star is barotropic or non-barotropic. (Note that Lindblom & Ipser call any inertial mode a generalized $`r`$ mode).

Let us start with the barotropic case for $`m=2`$. As already mentioned above, the axial equation for $`l=2`$ yields the band of frequencies given by Eq. (124). This is depicted in the left graph of Fig. 3, where we show the continuous spectrum as a function of $`M/R`$. No inertial modes are present. We have also included the Newtonian value for the $`r`$ mode as given by Eq. (123). Note that this value is completely independent of any stellar structure parameters. If we now include the coupling to the polar equations with $`l=3`$, the situation changes drastically. We now obtain two continuous frequency bands and a mode which lies below the Newtonian $`r`$ mode value, the further below the more compact the star is. Translating to positive frequencies, this means that this mode always has a higher frequency than the Newtonian $`r`$ mode. Inspection of the eigenfunction reveals that the axial part fits nicely to a $`r^{l+1}`$ power law, i.e.

$`u_3^{l=2}`$ $``$ $`r^3,`$ (125)
which is exactly what one expects for the $`r`$ mode. Inclusion of higher $`l`$ only gives minor frequency corrections, which tells us that this is indeed a physical inertial mode, the relativistic version of an $`m=2`$ $`r`$ mode. According to our nomenclature for inertial modes,the $`r`$ mode is actually an $`i_1^3`$ inertial mode. In Fig. 3, we have also included the coupling up to $`l_{\mathrm{max}}=4`$ and 5. Here, each inclusion of one additional $`l`$ leads to another patch in the continuous spectrum. Notice that the width of the uppermost and lowermost patches is very broad, whereas the middle patches around $`\sigma /\mathrm{\Omega }_L1`$ remain very narrow, even for large $`M/R`$. For each $`l_{\mathrm{max}}`$, we find $`l_{\mathrm{max}}m`$ new inertial modes, appearing very close to those frequencies which were occupied by the continuous spectrum for $`l_{\mathrm{max}}1`$. Again, some modes exist throughout the whole compactness range, while some are destroyed by the continuous spectrum above some $`M/R`$. In particular this is the case for those modes residing at both ends of the frequency range, as the continuous bands become very broad for large $`M/R`$.

To assess the accuracy of our results, we have compared our $`r`$ mode frequencies with non-linear time evolutions. For the $`1.4`$M neutron star model with $`\mathrm{\Omega }=2180`$s<sup>-1</sup>, the non-linear evolution in the Cowling approximation yields an $`r`$-mode frequency of 500 Hz (Stergioulas, 2002). Our linear calculation gives a frequency of $`512`$ Hz, which differs only by 2.5 per cent from the non-linear value. For larger rotation rates the difference increases, and for $`\mathrm{\Omega }=4986`$s<sup>-1</sup>, the non-linear value is $`1030`$Hz, whereas we find $`1172`$Hz, which is by 14 per cent off.

When turning to the non-barotropic case the picture becomes more complicated because of the presence of the $`g`$ modes. When the star is set into rotation, the $`g`$ and inertial modes start to interact. Yoshida & Lee (2000b) have studied in depth the behaviour of the $`g`$ and inertial modes as functions of the rotation rate and $`\mathrm{\Gamma }_1`$ for Newtonian stars. In this paper, we would like to focus only on the influence of the non-barotropicity on the $`r`$ mode.

In accordance with the results Yoshida & Lee (2000b) we find that the $`r`$ mode frequency remains almost unaffected as the star starts to deviate from being barotropic. However, we find that for small rotation rates, the eigenfunctions change in a quite drastic way if $`\mathrm{\Gamma }_1`$ starts to differ from $`\mathrm{\Gamma }`$. The larger the rotation, the more similar they get. We have already mentioned above that the relativistic $`r`$ mode of a barotropic star has to be an axial-led hybrid mode, which must retain its polar contribution in the non-rotating limit, whereas for a non-barotropic star the polar contribution must vanish. This can be readily checked by computing a sequence of eigenfunctions for decreasing values of $`\mathrm{\Omega }`$. These results are shown in Fig. 4, were we contrast the $`m=2`$ $`r`$ modes for a barotropic model with those of the corresponding non-barotropic model with $`\mathrm{\Gamma }_1=3`$. We only show the dominant parts of the eigenfunctions, which are the axial contribution $`u_3`$ with $`l=2`$ and the polar parts $`H`$, $`u_1`$ and $`u_2`$ with $`l=3`$. Furthermore, we rescale the eigenfunctions in such a way that $`u_3(R)=1`$.

For the barotropic sequence (left panels), we can see that the shape of $`u_3`$ does not depend on the rotation rate and is always given by $`r^3`$. The polar parts $`H`$, $`u_1`$ and $`u_2`$ actually do change along the sequence. Whereas $`H`$ goes to zero, $`u_1`$ and $`u_2`$ approach some finite values. This is in full agreement with Lockitch et al. (2001), who have shown that the non-trivial stationary perturbations must have vanishing $`H`$ and non-vanishing $`u_1`$ and $`u_2`$. From Eq. (32) is clear that for vanishing $`H`$ (and therefore vanishing $`\dot{H}`$), $`u_1`$ and $`u_2`$ must satisfy

$`u_1^{}+\left(2\nu ^{}\lambda ^{}+{\displaystyle \frac{2}{r}}\right)u_1e^{2\lambda }{\displaystyle \frac{l(l+1)}{r^2}}u_2`$ $`=`$ $`0,`$ (126)
which is equivalent to Eq. (3.23) of Lockitch et al. (2001).

The non-barotropic sequence (right panels) is completely different. For a large rotation rate $`\epsilon =0.2`$, the eigenfunction is very similar to the barotropic one. However, for smaller rotation rates, they begin to deviate considerably. The shape of the axial part $`u_3`$ is no longer described by $`r^3`$, but now becomes zero throughout most parts of the star and only peaks close to the surface. Our chosen value of $`\mathrm{\Gamma }_1=3`$ seems to be quite large, however, the same effect can be observed for $`\mathrm{\Gamma }_1`$ much close to 2, in this case, it appear for much smaller rotation rates. In principle, we should not expect this behaviour, instead the shape of $`u_3`$ should remain the same and only the polar perturbations should go to zero. From Fig. 4 we can see that latter is actually the case and $`H`$ as well as $`u_1`$ and $`u_2`$ all go to zero in the non-rotating limit. We can track this singular behaviour of $`u_3`$ back to the frame dragging $`\omega `$, for if we set it to zero, thus simulating a Newtonian stellar model, the eigenfunction of $`u_3`$ is well behaved; i.e. it remains proportional to $`r^3`$ for any rotation rate, whereas all the polar perturbations go to zero in the non-rotating limit. We suppose that this singular behaviour of $`u_3`$ is an artefact of the Cowling approximation, which should be cured by taking into account the metric perturbations.

## 4 Summary

We have investigated the inertial modes of slowly rotating relativistic stars in the Cowling approximation. Considering only first order rotational corrections, we have performed our computations in both the time and the frequency domain. Comparison with non-linear results for the $`f`$, $`p`$ and $`r`$ modes showed very good agreement, at least for slowly rotating stars. For larger rotation, second order effects become important and would have to be included for accurate mode calculations.

For the numerical calculations, the system of equations has to be truncated at some cut-off value $`l_{\mathrm{max}}`$. This truncation, however, introduces many spurious modes, which either vanish or shift in frequencies when $`l_{\mathrm{max}}`$ is increased. The “true” inertial modes, which appear for a given $`l_{\mathrm{max}}`$ are those which do not change significantly when $`l_{\mathrm{max}}`$ becomes larger. For each $`l>m>0`$, one can find $`lm`$ inertial modes, which can be labeled $`i_n^l`$ ($`1nlm`$). For the axisymmetric case $`m=0`$, there are $`l1`$ different modes, although in this case one mode has zero frequencies and the others always come in pairs with positive and negative frequencies.

We have shown that, as a result of the relativistic frame dragging, there always exists a continuous spectrum, which, however, depends very strongly on the number of coupled equations. For $`l_{\mathrm{max}}=2`$, the continuous spectrum is confined to a single connected region, with an increasing width for more relativistic stellar models. If $`l_{\mathrm{max}}`$ is increased the continuous spectrum splits into an increasing number of patches, which are disconnected for weakly relativistic stars but tend to overlap for more compact stellar models. For $`m0`$, each successive inclusion of a higher $`l`$ creates a new patch, whereas for $`m=0`$, only every second $`l`$ does so. In the limit $`l_{\mathrm{max}}\mathrm{}`$ the continuous spectrum should completely cover the frequency range of the inertial modes. In previous work, which neglected the coupling between the equations, it was stated that modes should not be able to exist inside the continuous spectrum, as the eigenfunctions would become singular. Through the time evolutions, however, we were able to show that there are modes, which can exist inside the continuous spectrum.

We provide the following (phenomenological) explanation. If an inertial mode can exist inside the continuous spectrum for a given $`l_{\mathrm{max}}`$, then it had to appear for the first time outside the continuous spectrum for some smaller $`l_{\mathrm{max}}`$. If a mode with label $`i_n^{l_0}`$ appears outside the continuous spectrum for $`l_{\mathrm{max}}=l_0`$, then it will not be affected if the continuous spectrum shifts itself on top of the mode for some $`l_{\mathrm{max}}>l_0`$. On the other hand, if the continuous spectrum for $`l_{\mathrm{max}}=l_0`$ is already so broad that it covers the position where the $`i_n^{l_0}`$ mode would appear, then this mode cannot come into existence.

Hence, it seems that for very relativistic stars, the continuous spectrum is able to destroy some of the inertial modes. We should, however, be still very cautious to draw any conclusions about the (non)existence of these modes for a rapidly rotating neutron star. We have seen that the continuous spectrum is very sensitive to the number of equations that are coupled. As we are only taking into account the first order rotational corrections, there is only coupling from $`l`$ to $`l\pm 1`$. Inclusion of second order corrections leads to coupling to $`l\pm 2`$; however, this should not affect the continuous spectrum, since in the relevant equations for $`u_2`$ and $`u_3`$, only terms coupling to $`H`$ and $`u_1`$ arise. Third order corrections might do the job and it could well be that they will affect the continuous spectrum in such a way that modes which do not exist in the first order analysis can now exist because the responsible patch could have moved somewhere else.

So far, we have completely neglected any metric perturbations, as we worked in the Cowling approximation. Lockitch et al. (2001) included some (non-radiative) metric perturbations, but restricted their studies to a post Newtonian treatment. Nevertheless it is clear that we should expect some quantitative difference when including these terms, i.e. when we move from the Cowling approximation to the low-frequency approximation. Still, the continuous spectrum will remain unaffected when metric variables are included, for it is only the fluid equations which are responsible for its existence. By including the metric perturbations, the modes might be affected in such a way, however, that they get pushed out of the continuous spectrum and we might find some of the modes which we cannot find in the Cowling approximation. We think that this scenario seems to be quite likely as the study of the purely axial $`l=m`$ case has shown. There, the Cowling approximation leads to a purely continuous spectrum without any mode solutions. In the low-frequency approximation, where one retains a certain metric component, one is lead to Kojima’s master equation, which in certain cases admits $`r`$ mode solutions, whose frequencies lie outside the continuous spectrum (Ruoff & Kokkotas 2001; Yoshida 2001). But the range of the continuous spectrum was still the same as in the Cowling approximation and even when all the metric perturbations were included (Ruoff & Kokkotas 2002; Yoshida & Futamase 2001), the qualitative picture remained the same. Hence, we conjecture that the inclusion of the metric perturbations in the fully coupled equations will not affect the behaviour of the continuous spectrum, but it certainly will shift the modes, and possibly in such a way that some of them might be pushed outside the patches of the continuous spectrum and therefore reappear.

As far as the $`r`$ modes are concerned, our results show that their existence does not seem to called into question by the continuous spectrum. The previous results by Ruoff & Kokkotas (2001, 2002), Yoshida (2001) and Yoshida & Futamase (2002), which suggested that in certain stellar models with rather soft equations of state, the (purely axial) relativistic $`r`$ modes cannot exist due to the existence of the continuous spectrum, seem to be artefacts arising from neglecting the coupling to the polar equations. The inclusion of this coupling shifts the location of the continuous spectrum and the $`r`$ mode can exist well outside it for a large range of stellar models. In this paper, we have only presented results for a single polytropic equation of state with $`\mathrm{\Gamma }=2`$, but from our previous results (Ruoff & Kokkotas 2001), we know that for $`\mathrm{\Gamma }>2`$, the width of the continuous spectrum actually shrinks and thus favours the presence of inertial modes. As most of the current realistic equations of state are stiffer than a $`\mathrm{\Gamma }=2`$ polytrope, we expect that the existence of $`r`$ modes should not be questioned. The same should also hold for strange stars.

Nevertheless it is clear that the relativistic case is quite different from the Newtonian one, and further studies are necessary to gain a full understanding of the inertial modes of (rapidly) rotating relativistic stars. In particular the metric perturbations have to be included to obtain the correct $`r`$-mode frequencies and the associated growth times. We cannot rule out that, in particular the latter will deviate significantly from the results obtained from applying the quadrupole formula to Newtonian stellar models.

## ACKNOWLEDGMENTS

We wish to thank Nikolaos Stergioulas and Ulrich Sperhake for useful comments. J.R. is supported by the Marie Curie Fellowship No. HPMF-CT-1999-00364. A.S. was supported by the Greek National Scholarship foundation (I.K.Y.) during this work.This work has been supported by the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).


text:
# The near-infrared Hubble diagram for sub-mm galaxies

## 1. Exploring the nature of the sub-mm galaxy population

There are various routes by which to address the question of whether or not sub-mm selected sources are the progenitors of massive elliptical galaxies.

One approach is to determine the extent to which known high-redshift populations are capable of producing strong mm–far-infrared emission. This (relatively efficient) experiment has been performed for powerful radio galaxies (Archibald et al. 2001), quasars (e.g. Isaak et al. 2002), Lyman-break galaxies (e.g. Chapman et al. 2000), and X-ray selected AGN (e.g. Page et al. 2001). The results of these studies indicate that the sort of massive starburst (SFR $`5001000\mathrm{M}_{}\mathrm{yr}^1`$) required to produce a bright ($`S_{850\mu m}>4`$ mJy) sub-mm source is found only in the most massive objects, and generally only at high redshift ($`z>1.52`$). However, it is clearly possible that the sub-mm selected population could be dominated by a somewhat different class of source.

A second approach is to compare the estimated comoving number density of bright sub-mm sources with those of other galaxy populations. As discussed by Scott et al. (2002) the comoving number density of sources with $`S_{850\mu m}>8`$ mJy is $`1\times 10^5\mathrm{Mpc}^3`$, comparable to that of extremely red objects at $`z1.5`$, and present-day ellipticals with $`L>23L^{}`$. However, the interpretation of such approximate coincidences depends on the assumed typical duration (and frequency) of the starbursts which power the sub-mm emission.

A third approach is to determine how the clustering properties of the sub-mm population compare with those displayed by other galaxy populations. Unfortunately current sub-mm and mm surveys are too small to yield a statistically significant detection of source clustering. However, existing data are certainly consistent with clustering as strong as that displayed by, for example, extremely red objects (Scott et al. 2002), and marginally significant evidence of clustering has been gleaned from cross correlation with Chandra sources (Almaini et al. 2002) and Lyman-break galaxies (Webb et al. 2002).

A fourth approach, which I explore in this brief article, is to investigate how the host galaxies of sub-mm sources compare with other known high-redshift galaxies. Here I have chosen to make this comparison in the near-infrared (i.e. $`K`$-band) both because many of the sub-mm source host galaxies have only been detected at $`K`$, and because observed $`K`$-band brightness currently offers the best means to estimate the masses of high-redshift galaxies.

## 2. The $`Kz`$ diagram

The $`Kz`$ diagram for radio galaxies has been studied and augmented ever since Lilly & Longair (1984) first demonstrated the existence of a tight relation for the 3CR galaxies out to $`z1`$. Recently Eales et al. (1997), van Breugel et al. (1998) and Jarvis et al. (2002) have extended the observed relation to $`z>4`$. As shown in Fig. 1, and as discussed by Jarvis et al. (2002), the radio-galaxy data are consistent with the $`Kz`$ relation expected for a passively-evolving galaxy of (present-day) luminosity $`45L^{}`$ out to the highest redshifts. Whether pure passive evolution is the true explanation of this relation is not really important for the present purpose. What is clear is that powerful radio galaxies are the most massive known galaxies at high redshift and are undoubtedly destined to be massive ellipticals (given they contain very massive black holes). They thus provide a clear benchmark against which to test the hypothesis that sub-mm selected sources are the progenitors of today’s massive ellipticals.

Before attempting to place sub-mm sources on this diagram it is important to realize that the $`K`$ magnitudes of the radio galaxies have all been measured through large apertures; in Fig. 1 the 3CR and 6C values have been corrected to a metric aperture of 50 kpc, while the $`K`$ magnitudes given by van Breugel et al. (1998) have been measured through an (in effect very similar) aperture of diameter 4 arcsec. This is important because the continued low scatter in the radio galaxy relation beyond $`z2.5`$ would not be produced if, for example, small aperture magnitudes were simply measured for the brightest visible clumps in these often complex sources. In considering sub-mm sources I have therefore included only objects with measured or reasonably solid estimated redshifts, and for which 4-arcsec aperture or near-total $`K`$ magnitudes have been measured (or at least attempted) armed with sub-arcsec positional accuracy.

The first three sub-mm sources I have plotted in Fig. 1 are therefore the three galaxies with spectroscopic redshifts from the cluster lens survey, namely SMMJ02399-0134 ($`z=1.06`$), 14011+0252 ($`z=2.55`$), and 02399-0136 ($`z=2.80`$) (Ivison et al. 2000, Smail et al. 2000). The $`K`$ magnitudes of these three lensed sources have been de-magnified by factors of 2.5, 2.8 and 2.5 respectively. To explore the relation further I have added the 4 best-studied bright sub-mm sources from the CUDSS and 8-mJy SCUBA surveys (Eales et al. 1999; Scott et al. 2002), all of which have very accurate positions (from IRAM PdB and VLA interferometry) yielding reliable $`K`$-band identifications in deep near-infrared images. These are CUDSS14A (Gear et al. 2000), Lockman850.1 (Lutz et al. 2001), ELAISN2850.1 and ELAISN2850.2 (Ivison et al. in prep). The redshifts of these objects have been estimated from the mean relation between 350GHz/1.4GHz flux-density ratio and $`z`$ given by Carilli & Yun (2000) which, in the redshift range of interest, is effectively identical to the extrapolated relation for Arp 220. Finally, I have included a lower limit at $`z=4.5`$ to indicate the average position of the 3 millimetre sources recently discussed by Dannerbauer et al. (2002). These sources merit inclusion because, while they still lack $`K`$-band identifications, they again have positions from mm-wave interferometry of sufficient accuracy to be sure that their counterparts lie below the detection threshold of the existing $`K`$-band imaging (to convert the published 2-arcsec diameter detection limits to 4-arcsec values I have simply doubled the noise).

Despite the obvious uncertainties it can be seen from Fig. 1 that the sub-mm galaxies display a plausible relationship between $`K`$ and $`z`$. Moreover this $`Kz`$ relation appears indistinguishable to that displayed by the radio galaxies, at least out to $`z3`$. At higher redshifts there is a suggestion that sub-mm galaxies are signficantly fainter at $`K`$ than radio galaxies, but at present it is unclear whether this indicates a significant difference in stellar mass or the increasing impact of dust obscuration on the rest-frame light from the sub-mm hosts.

## References

Almaini, O. et al., 2002, MNRAS, in press (astro-ph/0108400)

Archibald, E.N. et al. 2001, MNRAS, 323, 417

Carilli, C. & Yun, M. 2000, ApJ, 530, 618

Chapman, S.C. et al. 2000, MNRAS, 319, 318

Dannerbauer, H. et al. 2002, ApJ, submitted (astro-ph/0201104)

Eales, S.A. et al. 1997, MNRAS, 291, 593

Eales, S.A. et al. 1999, ApJ, 515, 518

Gear, W.K. et al. 2000, MNRAS, 316, L51

Isaak, K.G. et al. 2002, MNRAS, 329, 149

Ivison, R.J. et al. 2000, MNRAS, 315, 209

Jarvis, M.J. et al. 2002, in press (astro-ph/0112341)

Lilly, S.J. & Longair, M.S. 1984, MNRAS, 211, 833

Lutz, D. et al. 2001, A&A, 378, 70

Page, M.J. et al. 2001, Sci, 294, 2516

Scott, S. et al. 2002, MNRAS, in press (astro-ph/0107446)

Smail, I. et al. 2000, ApJ, 528, 612

van Breugel, W.J.M., et al. 1998, ApJ, 502, 614

Webb, T.M.A. et al. 2002, ApJ, in press (astro-ph/0201181)


text:
# LIFTS: an Imaging Fourier Transform Spectrograph for Astronomy

## 1. Introduction

We present the first astronomical results of LIFTS (see Wurtz et al. 2001 for details). The double output Michelson interferometer allow the observation of a spectrum for each of the 1k x 1k pixels of the combined images.

In a collimated beam, the incoming light hits a beam splitter and follows two possible optical paths before recombining into two output CCD cameras. All the photons are therefore contained in the sum of the two images while the relative phase shift information is stored in their difference. When stepping through a range of path length differences, by moving one arm of the interferometer, one ends up with the Fourier transform of the spectrum for each element of the images.

### 1.1. The Unique Advantages of an IFTS

By first principles, an FTS has a highly flexible resolution from $`=1`$ to several thousands, only limited by the maximum travel of the moving arm ( 1 cm in our case). Every pixel contains a spectrum, and all the images can be co-added to obtain a very deep image, called the “panchromatic” image. Follow-up spectroscopy of faint objects discovered in the panchromatic image can be directly obtained from the datacube. The observer can even decide to increase the spectral resolution in real-time by increasing the number of steps to a higher path length difference. We stress that this technique does not require preselecting objects or subregions in a few-arcminute field to be fed to a spectrograph.

Compared to any dispersive system, where all the spatial and spectral elements have to fit on a single CCD, the IFTS has a substantial multiplex advantage. Not only is the full CCD available for spatial information, but the spectral information is just limited by the number of channels one is willing to take on a given field. This multiplex advantage is proportional to the density of relevant pixels in the image (including some amount of sky and calibration objects).

The IFTS must be distinguished from filter systems including imaging Fabry-Perot (FP) spectrographs or tunable filters. An IFTS accept all the light at all wavelengths. It only creates one interference between the two “branches” of a photon and records the phase difference. On the other hand, an FP creates multiple “branches” by multiple reflections inside a cavity creating constructive interferences only for a narrow band of wavelengths, rejecting most of the light. Also, the passband is periodic and one order has to be selected, limiting the spectral range to a few angstroms. F-P systems are better suited for narrow band observations of faint extended objects as they cut all the out-of-band noise to which IFTS would be sensitive.

## 2. Observations with the Livermore IFTS

We have had a total of seven clear bright-time half-nights on the APO 3.5-meter with the completed instrument. Seeing ranged from 0.7 to 1.5 arcsec, so we binned to $`512`$ pixels so that our pixel size was 0.28 arcsec.

The Ring Nebula datacube, at $`430`$, is presented in Figure 1.

## References

Bennett, C. L., Carter, M. R., & Fields, D. J. 1995, Proc. SPIE, 2252, 274

Wurtz, R., Cook, K. H., Bennett, C. L., Bixler, J., Carr, D.& Wishnow, E. H. 2000, ASP Conf. Series, 207., 203