format: text

text:
# Memory Efficient Arithmetic

## 1 Introduction

Suppose that $`\alpha `$ is a positive real number. Then, a standard fact regarding base-$`b`$ representations is that given any integer $`b2`$, there exists and integer $`J`$, and a sequence of integers $`r_J,r_{J1},\mathrm{}`$, such that

$$\alpha =\underset{jJ}{}r_jb^j,0r_jb1,r_J0,$$
(1)
and we write

$$\alpha =(r_Jr_{J1}\mathrm{}r_0.r_1\mathrm{})_b$$

to denote this expansion (when $`J>0`$). If we further disallow $`r_n=b1`$ for all $`nN`$, for some $`N`$, then the sequence of $`r_i`$’s is unique. If we make no such restrictions on $`r_n`$, then the sequence of $`r_i`$’s may not be unique, as in the base-$`10`$ expansion

$$1=0.9999999\mathrm{}$$

A question which has attracted recent attention (see ) is the following: Given a base $`b`$, and an integer $`m0`$, efficiently determine a good approximation to

$$\nu =\{b^m\pi \},$$

where for a real number $`\theta `$, the notation $`\{\theta \}`$ means

$$\{\theta \}=\theta \theta (\mathrm{Note}:\{\theta \}\theta (mod1),0\{\theta \}<1).$$

We will say that such an approximation $`\gamma `$ is a level-$`p`$ approximation if and only if

$$|\gamma \nu |<\frac{1}{b^p}.$$

Now, if the approximation $`\gamma `$ is sufficiently good, then one can use it to determine the digits of $`\pi `$: For example, if we take $`m=3`$ and $`b=10`$, then we get

$$\nu b^m\pi 10^3\pi 0.59265358979\mathrm{}(mod1).$$

Suppose we had an approximation $`\gamma `$ to $`\nu `$ satisfying

$$|\gamma \nu |<\frac{1}{1000};$$

so, $`\gamma =0.59\mathrm{}`$. Then, the leading base-$`10`$ digit of $`\gamma `$ is the same as the fifth base-$`10`$ digit of $`\pi `$ (from the left), which is $`5`$.

More generally, a good approximation to $`\nu `$ gives us the $`(m+2)`$nd digit of $`\pi `$; however, depending on the value of $`m`$ selected, this approximation may need to be extremely close to $`\nu `$, in order to determine this digit. For example, in the above instance with $`m=3`$, if $`\gamma =0.6`$, then

$$|\gamma \nu |<\frac{1}{100},$$

and we note that the leading digit of $`\gamma `$ is not the same as the leading digit of $`\nu `$.

In this paper we will describe a method for determining the $`m`$th base-$`b`$ digit of an integer, where $`m=1`$ corresponds to the least significant digit (note that for the digits-of-$`\pi `$ problem above, $`m=1`$ corresponded to the leading digit; so, the $`m`$th digit is defined differently in this context). As in the problem concerning digits of $`\pi `$, this method produces arbitrarily precise approximations to

$$\nu =\left\{\frac{n}{b^m}\right\},$$

where $`n`$ is this integer, which will be defined by a certain type of expansion that we will describe below. Note that the value of $`m`$ is at least $`1`$ (if $`m0`$, then $`\nu `$ is trivially $`0`$). Now, if

$$\nu =(0.r_1r_2\mathrm{})_b,$$

then the $`m`$th digit of $`n`$ equals $`r_1`$; and, if $`\lambda `$ is a sufficiently good approximation to $`\nu `$, then the leading digit of $`\lambda `$ will also equal $`r_1`$.

The type of expansion for $`n`$ we will use is defined as follows: In , Smale and Shub say that a computation of length $`L`$ of a positive integer $`n`$ is a sequence of integers $`s_1,s_2,\mathrm{},s_L`$, where $`s_1=0`$, $`s_2=1`$, and for $`i3`$,

$$s_i=s_js_k,\mathrm{where}j,k<i,$$

where $``$ is either addition, subtraction, or multiplication, and where $`s_L=n`$. If one expresses such a computation as a string indices $`(j,k)`$ and operations $`+,`$ and $`\times `$, then given an integer $`N2`$, one can compute $`n(modN)`$ using only $`O(L(\mathrm{log}L)(\mathrm{log}N)(\mathrm{log}\mathrm{log}N)^2)`$ bit operations, by just computing the sequence $`s_1,\mathrm{},s_L`$ modulo $`N`$. The factor $`(\mathrm{log}N)(\mathrm{log}\mathrm{log}N)^2`$ in this big-O appears because a product of integers modulo $`N`$ can be computed using Fast Fourier Transforms using only $`O((\mathrm{log}N)(\mathrm{log}\mathrm{log}N)^2)`$ bit operations.

We will say that an integer $`n`$ has computational complexity $`L`$ if and only if there exists a computation of lenth $`L`$ for computing $`n`$.

In the next section we will prove a general result, which we will use to prove the following theorem:

###### Theorem 1

Given a positive integer $`n`$ having computational complexity $`L`$, an approximation $`A`$ to the number of base-$`b`$ digits of $`n`$ (see the input specs below), an integer $`m0`$, and a level $`y`$, there exists an algorithm for computing a level-$`y`$ approximation to

$$\nu =\left\{\frac{n}{b^m}\right\}.$$

This algorithm requires only

$$O(yL(\mathrm{log}^{2/3}n)\mathrm{log}^C(b+m+y+L+\mathrm{log}n))\mathrm{bits}\mathrm{of}\mathrm{memory},$$

and

$$O(yL(\mathrm{log}n)\mathrm{log}^C(b+m+y+L+\mathrm{log}n))\mathrm{bit}\mathrm{operations}.$$

The input and output requirements of this algorithm are as follows:

> Input: The integers $`m`$, $`y`$, and a string representing the length-$`L`$ computation of $`n`$. Also, the algorithm requires as input an integer $`A`$, which is an approximation to the number $`d`$ of base-$`b`$ digits of $`n`$. This approximation need only satisfy
>
> $$\frac{1}{2}<\frac{A}{d}<2.$$
>
> We further restrict $`m`$ so that $`m<2A+y+1`$, since otherwise $`\gamma =0`$ satisfies the conclusion of our Theorem.

> Output: The level-$`y`$ approximation to $`\nu `$, encoded as a string of $`y+1`$ base-$`b`$ digits.

This theorem requires a little more explanation. First of all, the input to the algorithm will be a string of $`O(L\mathrm{log}L+\mathrm{log}m+\mathrm{log}y)`$ bits, which is smaller than the space requirement listed above (when $`C>1`$). The $`L\mathrm{log}L`$, $`\mathrm{log}m`$, and $`\mathrm{log}y`$ terms here account for the number of bits needed to specify the length $`L`$ computation of $`n`$, the index $`m`$, and the level $`y`$, respectively. The approximation $`A`$ to $`d`$ requires only $`O(\mathrm{log}\mathrm{log}n)`$ bits of space, and this turns out to be $`O(L\mathrm{log}L)`$: To see why this is so, we note that any length $`L`$ computation produces an integer $`n<2^{2^{L1}}`$, which can be proved by induction. It follows then that $`\mathrm{log}\mathrm{log}n=O(L)`$.

The output of the algorithm will be a string of $`O(y)`$ bits, representing the base-$`b`$ approximation $`\gamma `$ to $`\nu `$. The number $`\gamma `$ will have only $`y+1`$ base-$`b`$ digits, and will satisfy

$$|\gamma \nu |<\frac{1}{b^y}.$$

Perhaps the most surprising aspect of the above theorem is that the indicated algorithm can require significantly less memory to find the approximation to $`\nu `$ than it does to write down the number $`n`$, which will have $`O(\mathrm{log}n)`$ base-$`b`$ digits. The following corollary of the above theorem will make this point clear:

###### Corollary 1

Suppose that $`a,b2`$ and $`m,t,y1`$ are all integers. There exists an algorithm which computes a level-$`y`$ base-$`b`$ approximation $`\gamma `$ to $`\nu `$, where

$$\nu =\left\{\frac{a^t}{b^m}\right\}.$$

This algorithm requires only

$$O(yt^{2/3}\mathrm{log}^C(a+b+y+m+t))\mathrm{bits}\mathrm{of}\mathrm{memory},$$

where $`C>0`$, and performs

$$O(yt\mathrm{log}^C(a+b+y+m+t))\mathrm{bit}\mathrm{operations}.$$

Now, the number of bits needed to write down the number $`n=a^t`$ is clearly $`O(t\mathrm{log}a)`$; and yet, if, say, we take $`y=1`$, this algorithm requires only $`t^{2/3+o(1)}`$ bits of space.

This corollary follows since $`a^t`$ has computational complexity

$$L=O(\mathrm{log}^2(a+t)).$$

To see this, we note that $`a^t`$ can be generated by repeated squaring: If $`t=2^{t_1}+\mathrm{}+2^{t_s}`$, then $`a^t=a^{2^{t_1}}\mathrm{}a^{2^{t_s}}`$. These numbers $`a^{2^h}`$ can be computed by starting with $`a`$; then squaring to get $`a^2`$; then squaring again to get $`a^4`$; then continuing, this produces the list $`a,a^2,a^4,\mathrm{},a^{2^h}`$ after only $`h`$ multiplications.

The rest of this paper is organized as follows: In the next section we will state the Main Theorem (Theorem 2) and then use it to deduce Theorem 1. In section 3 we give a proof of the Main Theorem. Finally, in section 4 we give a proof of a proposition (Proposition 1), which is an auxillary result needed for the proof of the Main Theorem.

## 2 Main Theorem and Proof of Theorem 1

Theorem 1 is actually a corollary of a more general result concerning approximations to $`\nu `$. In this section we will state this result, which will henceforth be called the Main Theorem, and then show how to apply it to prove Theorem 1. The proof of this theorem, as well as a brief description of the ideas used to prove it, can be found in section 3; also, at the end of Section 3.1, we will give a brief statement on how to parallelize the algorithm.

We suppose that $`b2`$ is an integer, which is to be the base used; that $`0<n/a^t<1`$ is some rational number where $`n,a1`$, $`t0`$ are integers, and where $`n`$ has computational complexity $`L`$; that $`y1`$ is some level of precision to be used; and finally, that $`\mu 0`$ is some integer. Then, given any pair of integers $`S,T`$ satisfying

$$ST>\frac{3(\mathrm{log}n+(\mu +y+2)\mathrm{log}b)}{\mathrm{log}a},a^S>T^2,$$
(2)
we have the following

###### Theorem 2 (Main Theorem)

Let

$$\nu =\left\{b^\mu \frac{n}{a^t}\right\}.$$
(3)
There exists an algorithm which computes a level-$`y`$ approximation $`\gamma `$ to this number $`\nu `$, where the space and time requirements of the algorithm are as follows:

> Space: $`O(yL(S+T)\mathrm{log}^C(y+L+S+T+a+b+\mu +\mathrm{log}t))`$ bits of memory, where $`C>0`$.

> Time: $`O(yL(ST+T^3)\mathrm{log}^C(y+L+S+T+a+b+\mu +\mathrm{log}t))`$ bit operations.

Here we give more precise information about the input and output specifications of the algorithm:

> Input: The positive integers $`a,t,\mu ,y,S`$ and $`T`$, as well as a string of $`O(L\mathrm{log}L)`$ characters representing the length-$`L`$ computation needed to produce $`n`$.

> Output: The algorithm will give an approximation $`\gamma `$ to $`\nu `$. This approximation will have $`y+1`$ base-$`b`$ digits, and will satisfy
>
> $$|\gamma \nu |<\frac{1}{b^y}.$$

To prove Theorem 1, using this result, we let $`a=b`$, $`t=2A`$, and $`\mu =tm`$. We note that this gives

$$\frac{b^\mu n}{a^t}=\frac{n}{b^m}.$$

We also let

$`S`$ $`=`$ $`\left({\displaystyle \frac{3(\mathrm{log}b)(3A+\mu +y+2)}{\mathrm{log}a}}\right)^{2/3}+1,`$
$`T`$ $`=`$ $`\left({\displaystyle \frac{3(\mathrm{log}b)(3A+\mu +y+2)}{\mathrm{log}a}}\right)^{1/3}+1.`$
We note that this choice of $`S`$ and $`T`$ satisfies (2). 111To show this, one needs the fact that $`3A\mathrm{log}b>\mathrm{log}n`$, which follows since $`n`$ has $`2A`$ base-$`b`$ digits.

Now, applying the algorithm described in Theorem 2 with the parameters indicated above, we get the same output as described in Theorem 1. The running time and space requirements to run this algorithm are also as stated in Theorem 1 for our particular choices of $`S`$ and $`T`$.

## 3 Proof of Theorem 2

Let $`S`$ and $`T`$ be as in (2), and let $`r`$ and $`k`$ be integers such that

$$t=Skr,0rS1.$$

Then, we have that

$$\alpha =\frac{n}{a^t}=\frac{na^r}{a^{Sk}}.$$

The idea of the proof of Theorem 2 is to approximate $`b^\mu \alpha `$ (and therefore $`\nu `$) as follows:

$$b^\mu \alpha =\gamma _1+\mathrm{}+\gamma _T+E,$$
(4)
where

$$\gamma _j=b^\mu \frac{A_j}{a^Sj},$$

for some rationals $`A_1,\mathrm{},A_T`$, and where

$$|E|<\frac{1}{b^{y+2}}\mathrm{for}ST\mathrm{sufficiently}\mathrm{large}.$$
(5)
Then, we will find approximations $`\gamma _1^{},\mathrm{},\gamma _T^{}`$ to $`\{\gamma _1\},\mathrm{},\{\gamma _T\}`$. Now, if the precision of these approximations $`\gamma _1^{},\mathrm{},\gamma _T^{}`$ is high enough, and if we we let $`\mathrm{\Sigma }`$ satisfy

$$\mathrm{\Sigma }=\{\gamma _1^{}+\gamma _2^{}+\mathrm{}+\gamma _T^{}\},$$

then $`\mathrm{\Sigma }`$ will be an approximation to $`\nu `$; and, if we then take $`\gamma `$ to be the closest number to $`\mathrm{\Sigma }`$ having $`y+1`$ base-$`b`$ digits, then $`\gamma `$ will be a level-$`y`$ approximation to $`\nu `$.

We claim that the approximations $`\gamma _1^{},\mathrm{},\gamma _T^{}`$ to $`\{\gamma _1\},\mathrm{},\{\gamma _T\}`$ need only have

$$w=y+\frac{\mathrm{log}T}{\mathrm{log}b}+3,$$

base-$`b`$ digits (and be level-$`w`$ approximations), in order to guarantee that $`\mathrm{\Sigma }`$ is a level-$`y+1`$ approximation to $`\nu `$. Note that this would imply that

$$|\{\gamma _j\}\gamma _j^{}|<\frac{1}{Tb^{y+2}}.$$

To see only $`w`$ base-$`b`$ digits are needed, we note that if these numbers $`\gamma _j^{}`$ satisfy this last inequality, then by the triangle inequality,

$`|\gamma \nu |`$ $``$ $`{\displaystyle \frac{1}{b^{y+1}}}+|\mathrm{\Sigma }\nu |{\displaystyle \frac{1}{b^{y+1}}}+{\displaystyle \underset{j=1}{\overset{T}{}}}|\gamma _j^{}\{\gamma _j\}|+|E|`$
$`<`$ $`{\displaystyle \frac{1}{b^{y+1}}}+{\displaystyle \frac{1}{b^{y+2}}}+{\displaystyle \frac{1}{b^{y+2}}}{\displaystyle \frac{1}{b^y}},`$
as claimed.

Let us now find a set of values for $`A_1,\mathrm{},A_T`$ which make (4) hold: Using the geometric series identity, we have that

$$\underset{i=1}{\overset{T}{}}\frac{A_i}{a^Si}=\underset{j=1}{\overset{\mathrm{}}{}}\frac{B_j}{a^{Sj}},$$
(6)
where

$$B_j=A_1+A_22^{j1}+A_33^{j1}+\mathrm{}+A_TT^{j1}.$$

We seek values for $`A_1,\mathrm{},A_T`$ so that

$$B_j=\{\begin{array}{cc}0,\hfill & \text{for }1jT,jk\text{;}\hfill \\ na^r,\hfill & \text{for }j=k\text{.}\hfill \end{array}$$

The following Proposition gives the solution we seek

###### Proposition 1

We have that

$$A_j=\frac{na^r(\mathrm{Coef}.\mathrm{of}x^{k1}\mathrm{in}\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}xh)}{_{\genfrac{}{}{0pt}{}{h=1}{hj}}^Tjh};$$
(7)
and,

$$|A_j|nTa^S4^T.$$
(8)
Note that this implies

$$\gamma _j=b^\mu \frac{(1)^{Tj}na^r(\mathrm{Coef}.\mathrm{of}x^{k1}\mathrm{in}\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}xh)}{(j1)!(Tj)!(a^Sj)}.$$

From this proposition we deduce that

$`|E|`$ $`=`$ $`b^\mu \left|{\displaystyle \underset{j=T+1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{A_1+A_22^{j1}+\mathrm{}+A_TT^{j1}}{a^{Sj}}}\right|`$
$``$ $`b^\mu {\displaystyle \underset{j=T+1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(nTa^S4^T)T^j}{a^{Sj}}}`$
$`=`$ $`{\displaystyle \frac{na^S4^TT^{T+2}b^\mu }{a^{S(T+1)}}}{\displaystyle \underset{j=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{T^j}{a^{Sj}}}`$
$`=`$ $`{\displaystyle \frac{na^S4^TT^{T+2}b^\mu }{a^{ST}(a^ST)}}<{\displaystyle \frac{nb^\mu }{a^{ST/3}}}<{\displaystyle \frac{1}{b^{y+2}}},`$
for $`ST`$ large enough; and so, (5) follows.

We now have all the ingredients necessary to prove Theorem 2, which we will give as the following algorithm:

### 3.1 Algorithm 1

The input, output, and requirements of this algorithm are as stated in Theorem 2. Here are the steps of the algorithm:

1. Let

$$w=y+\frac{\mathrm{log}T}{\mathrm{log}b}+3.$$

Note that this choice of $`w`$ satisfies

$$\frac{T}{b^w}\frac{1}{b^{y+2}}.$$

2. Set $`\mathrm{\Sigma }=0`$, and let $`r`$, $`S`$, $`T`$ and $`k`$ be as described at the beginning of this section.

3. For $`j`$ from $`1`$ to $`T`$ do steps 4 through 8.

4. Compute

$$Qa^Sj.$$

5. Set

$$v(j1)!(Tj)!Q.$$

This number can be computed using $`O((T+S)\mathrm{log}^C(T+S))`$ bit operations, and just as much memory (for some $`C>0`$).

6. Apply Algorithm 2 (given in the next subsection of the paper) to compute

$$H\mathrm{Coef}.\mathrm{of}x^{k1}\mathrm{in}\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}xh.$$

This step requires $`O(T^2\mathrm{log}^CT)`$ bit operations and $`O(T\mathrm{log}^CT)`$ bits of memory.

7. Compute

$$u(1)^{Tj}na^rH(modv),0uv1.$$

Since $`n`$ has computational complexity $`L`$, this step requires only $`O(L(S+T)\mathrm{log}^C(L+S+T+a))`$ bit operations, and just as much memory.

8. Find a number $`\tau `$ having $`w+1`$ base-$`b`$ digits satisfying

$$|\tau \varphi |<\frac{1}{b^w},$$

where

$$\varphi =\left\{b^\mu \frac{u}{v}\right\}.$$

(So, $`\tau `$ will be a level-$`w`$ approximation to $`\varphi `$.)

We note that this number $`\tau `$ can be easily computed by first letting

$$u_0b^\mu u(modv),0u_0v1,$$

and then noting that

$$\left\{\frac{b^\mu u}{v}\right\}=\frac{u_0}{v}=(0.r_1r_2\mathrm{})_b.$$

Then, by finding the first $`w+1`$ significant digits of $`u_0/v`$, and letting $`\tau =(0.r_1\mathrm{}r_{w1})_b`$ one see that the above inequalities are satisfied.

9. Set

$$\mathrm{\Sigma }\{\mathrm{\Sigma }+\tau \}.$$

We only need to do level-$`w`$ arithmetic in base-$`b`$ here.

(If $`j<T`$, then increment $`j`$ and loop back to step 4.)

10. (We assume $`j=T`$.) Let $`\gamma `$ be the number having $`y+1`$ base-$`b`$ digits which comes nearest to $`\mathrm{\Sigma }`$, and then OUTPUT $`\gamma `$.

We note that we can perform the operations in steps 3 through 8, with different values of $`j`$, in parallel. For example, given two processors, we can assign processor 1 to perform steps 3 through 8, with values of $`jT/2`$, and then assign processor 2 to do the same, but with $`T/2<jT`$. This would result in an two-fold reduction in the running time, as long as $`\mu `$ is sufficiently large. Of course, the memory requirements would double, because each of the two processors would require their own seperate memories.

More generally, we have that, given $`M`$ processors, for $`\mu `$ sufficiently large, Algorithm 1 can be computed in parallel, resulting in an $`M`$-fold reduction in running time, but an $`M`$-fold increase in memory requirements.

### 3.2 Algorithm 2

> Input: $`T,k,j`$.

> Output: Coef. of $`x^{k1}`$ in $`_{\genfrac{}{}{0pt}{}{h=1}{hj}}^Txh`$.

> Requirements: The algorithm performs $`O(T^2\mathrm{log}^CT)`$ bit operations (for some $`C>0`$), but requires only $`O(T\mathrm{log}^CT)`$ bits of memory.

1. Let $`P`$ be the least integer such that

$$\mathrm{\Delta }=\underset{\genfrac{}{}{0pt}{}{pP}{p\mathrm{prime}}}{}p2^{T+1}T!$$

Note: $`P=O(T\mathrm{log}T)`$, and can be computed using $`O(T\mathrm{log}^DT)`$ bit operations (for some $`D>0`$); and so, we can compute and store $`P`$ within the time and space requirements listed above for the algorithm. We also note that every coefficient of the polynomial in the output specifications is less than $`\mathrm{\Delta }/2`$ in absolute value.

2. Set $`\mathrm{\Sigma }=0`$.

3. For each prime $`pP`$ do steps 4 through 8.

4. Compute the polynomial

$$f(x)\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}xh(modp).$$

Note: This polynomial can be stored as a length-$`T`$ coefficient vector, and the number of bits required to store such a vector is $`O(T\mathrm{log}p)=O(T\mathrm{log}T)`$; also, this polynomial can be computed using $`O(T\mathrm{log}^DT)`$ bit operations by making use of FFT’s and a divide-and-conquer strategy for polynomial multiplication. The divide-and-conquer part of the algorithm can probably best be described as the following recursive procedure: First, we suppose that $`L`$ is a set of polynomials to be producted together modulo $`p`$, and Product$`(L)`$ denotes the procedure for computing this product. The pseudocode for this procedure is given as follows:

If $`|L|=1`$ (i.e. $`L`$ has only one polynomial), then

RETURN the contents of $`L(modp)`$;

Else, if $`|L|2`$, say $`L=\{f_1,\mathrm{},f_t\}`$, then

RETURN

Product$`(\{f_1,\mathrm{},f_{t/2}\})`$Product$`(\{f_{t/2+1},\mathrm{},f_t\})(modp)`$

Now, using FFT’s to perform the polynomial multiplication in this second step (the ‘Else’ step), we see that if the two polynomials being multiplied together have degrees $`\mathrm{}_1`$ and $`\mathrm{}_2`$, respectively, then the multiplication should take no more than $`O((\mathrm{}_1+\mathrm{}_2)\mathrm{log}^D(\mathrm{}_1+\mathrm{}_2+p))`$ bit operations. Now, if we run Product$`(L)`$ starting with $`L`$ consisting of all linear factors $`xh`$, $`1hT`$, $`hj`$, then if $`T1`$ is a power of $`2`$, the procedure products together $`(T1)/2`$ pairs of degree $`1`$ polynomials; $`(T1)/4`$ pairs of degree $`2`$ polynomials; and so on, all the way down to two polynomials of degree $`(T1)/2`$. So, the total number of bit operations required to run this producedure is

$$(\mathrm{log}^DT)\underset{j(\mathrm{log}T)/\mathrm{log}2+1}{}\frac{T}{2^j}2^j=O(T\mathrm{log}^DT).$$

The memory requirements (in bits) are likewise of the same order.

5. Set

$$H\mathrm{Coef}.\mathrm{of}x^{k1}\mathrm{in}f(x)(modp).$$

6. Set

$$N(\mathrm{\Delta }/p)^1H(modp),\mathrm{where}0Np1.$$

7. Set

$$\mathrm{\Sigma }\mathrm{\Sigma }+\frac{N\mathrm{\Delta }}{p}.$$

8. Increment the value of $`p`$, and return to step 4, unless $`p>P`$, in which case we proceed to step 9.

9. Let $`r`$ be the least residue in absolute value of $`\mathrm{\Sigma }(mod\mathrm{\Delta })`$.

10. Return the value of $`r`$, and STOP.

It is relatively easy to see that the algorithm requires no more than the indicated space and time requirements.

The idea behind the algorithm is that we use the Chinese Remainder Theorem to compute the $`x^{k1}`$ coefficient of our polynomial, and the computation in step 7 is just an “on the fly” CRT calculation. This calculation is based on the following fact: If $`q_1,\mathrm{},q_h`$ are coprime, and if $`a_1,\mathrm{},a_h`$ are any integers, then if we set

$$\mathrm{\Delta }^{}=\underset{i=1}{\overset{h}{}}q_i,$$

and

$$\mathrm{\Sigma }^{}=\underset{i=1}{\overset{h}{}}b_i\frac{\mathrm{\Delta }^{}}{q_i},\mathrm{where}b_ia_i(\mathrm{\Delta }^{}/q_i)^1(modq_i),$$

then

$$\mathrm{\Sigma }^{}a_i(modq_i),\mathrm{for}\mathrm{every}i=1,2,\mathrm{},h.$$

One might guess that the coefficient of our polynomial can be computed using less resources by using a “Fourier Series” method; that is,

$$\mathrm{Coef}.x^{k1}\mathrm{in}f(x)=\frac{1}{T}\underset{\mathrm{}=0}{\overset{T1}{}}e^{2\pi i\mathrm{}(k1)/T}f(e^{2\pi i\mathrm{}/T}).$$

It is not obvious (to me) how to do this without using the special form of the polynomial $`f(x)`$: First of all, we would need to maintain $`T`$ digits of precision for each term in the sum, since any particular coefficient of the polynomial $`f(x)`$ can have size $`2^{cT\mathrm{log}T}`$, for some $`c>0`$. Thus, $`T^2`$ bit operations would be needed to compute each term $`f(e^{2\pi i\mathrm{}/T})`$. In total, $`T^3`$ bit operations would be needed to evaluate all the terms in the sum. If one tries to use FFT’s to evaluate all the terms in the sum at the same time, this reduces the running time to $`O(T^2\mathrm{log}^DT)`$ bit operations; however, the memory requirements then increase to $`T^2`$ bits of storage, which is the amount needed to store the all numbers $`f(e^{2\pi i\mathrm{}/T})`$, $`0\mathrm{}T`$ to $`T`$ bits of precision. Even if we try a discrete version of this method, where the polynomials are computed, say, modulo $`2^k`$ for $`kT`$, and the roots of unity are roots of unity modulo $`2^k`$, we would run into the same difficulties.

## 4 Proof of Proposition 1

The $`A_i`$’s can be computed by solving the equation

$$\left(\begin{array}{ccccc}1& 1& 1& \mathrm{}& 1\\ 1& 2& 3& \mathrm{}& T\\ 1& 2^2& 3^2& \mathrm{}& T^2\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 1& 2^{T1}& 3^{T1}& \mathrm{}& T^{T1}\end{array}\right)\left(\begin{array}{c}A_1\\ A_2\\ A_3\\ \mathrm{}\\ A_T\end{array}\right)=\left(\begin{array}{c}0\\ \mathrm{}\\ 0\\ na^r\\ 0\\ \mathrm{}\\ 0\end{array}\right).$$
(9)
If we call the matrix on the left-hand-side $`M`$, then

$$A_j=na^rM_{j,k}^1,$$
(10)
where $`M_{j,k}^1`$ is the entry in the $`j`$th row, $`k`$th column of $`M^1`$.

We will calculate $`M_{j,k}^1`$ via polynomial interpolation: We have that for any set of ordered pairs

$$(1,b_1),(2,b_2),\mathrm{},(T,b_T),$$

where $`b_1,\mathrm{},b_T`$, there exists a unique degree $`T1`$ polynomial $`f(x)[x]`$ such that

$$f(i)=b_i,\mathrm{for}\mathrm{all}i=1,2,\mathrm{},T;$$

moreover, if we write

$$f(x)=c_Tx^{T1}+c_{T1}x^{T2}+\mathrm{}+c_2x+c_1,$$

then these coefficients $`c_i`$ can be calculated in two different ways: The first way is through basic linear algebra, since

$$\left(\begin{array}{ccccc}1& 1& 1& \mathrm{}& 1\\ 1& 2& 2^2& \mathrm{}& 2^{T1}\\ 1& 3& 3^2& \mathrm{}& 3^{T1}\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 1& T& T^2& \mathrm{}& T^{T1}\end{array}\right)\left(\begin{array}{c}c_1\\ c_2\\ c_3\\ \mathrm{}\\ c_T\end{array}\right)=\left(\begin{array}{c}b_1\\ b_2\\ b_3\\ \mathrm{}\\ b_T\end{array}\right).$$
(11)
We notice that the matrix on the left-hand-side is $`M^{}`$, the transpose of our matrix $`M`$.

The second way of calculating the $`c_i`$’s is by Lagrange interpolation, which gives

$$f(x)=\underset{i=1}{\overset{T}{}}b_i\underset{\genfrac{}{}{0pt}{}{h=1}{hi}}{\overset{T}{}}\frac{xh}{ih}.$$
(12)
Now, if we suppose that

$$b_i=\{\begin{array}{cc}0,\hfill & \text{if }ij\text{,}\hfill \\ 1,\hfill & \text{if }i=j\text{,}\hfill \end{array}$$

then for this choice of $`b_i`$’s, one sees from (11) that

$$c_k=(M^{})_{k,j}^1=M_{j,k}^1.$$

On the other hand, from (12) we see that

$$c_k=\mathrm{Coef}.\mathrm{of}x^{k1}\mathrm{in}\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}\frac{xh}{jh}.$$

Thus,

$$M_{j,k}^1=\mathrm{Coef}.\mathrm{of}x^{k1}\mathrm{in}\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}\frac{xh}{jh},$$

and we conclude from this and (10) that (7) holds.

Finally, to prove (8), we note that

$$\left|\underset{\genfrac{}{}{0pt}{}{h=1}{hj}}{\overset{T}{}}jh\right|=(Tj)!(j1)!,$$

The coefficient of $`x^{k1}`$ in the above polynomial is clearly less than

$$T!\left(\genfrac{}{}{0pt}{}{T}{k1}\right)<T!2^T.$$

So,

$$|A_j|na^S2^T\frac{T!}{(Tj)!(j1)!}=jna^S2^T\left(\genfrac{}{}{0pt}{}{T}{j}\right)<nTa^S4^T,$$

which proves (8).

## 5 Acknowledgements

I would like to thank Richard Hudson for an email he sent to me, which got me interested in these digit calculation questions, which eventually lead me to prove the theorems listed above. I would also like to thank Kevin Hare for pointing out to me that my algorithm above is highly parallelizable.


text:
# Inertial modes of neutron stars with the superfluid core

## 1 Introduction

Since Andersson (1998) and Friedman & Morsink (1998) found that the $`r`$-modes of a rotating neutron star are unstable against gravitational radiation reaction even for infinitesimally small rotation speeds, many papers on the $`r`$-mode instability have appeared. As a result of these investigations, the importance of the instability in astrophysics has been pointed out mainly in the following two ways, that is, the instability is strong enough to limits both the spin of newly born neutron stars and that of old accreting neutron stars, and neutron stars in which the $`r`$-mode oscillation is excited will be a potential source of gravitational waves that are detectable by the LIGO and other detectors (Lindblom, Owen, & Morsink 1998; Andersson, Kokkotas, & Schutz 1999; Owen et al. 1998; Bildsten 1998; Andersson, Kokkotas, & Stergioulas 1999). For a recent review, see, e.g., Andersson & Kokkotas (2001).

It has long been suggested that neutrons in the inner crust and neutrons and protons in the core of neutron stars are in superfluid states when the interior temperatures cool down below $`T10^9\mathrm{K}`$ (e.g., Shapiro & Teukolsky 1983). Since the interior temperature of neutron stars is believed to cool down quickly by emitting neutrinos (e.g., Baym & Pethick 1979), we may expect that many observable neutron stars, except newly born ones, are likely to have a core that contains superfluids. Although one of the important ingredients that affect the oscillations of neutron stars is superfluidity in the interior, it was rather recently that Epstein (1988) pointed out the importance of superfluidity for the oscillations, suggesting the existence of a new family of acoustic modes associated with superfluidity by means of the local dispersion relations derived for the case of non-rotating neutron stars. Later on, Lindblom & Mendell (1994) studied the $`f`$-modes of rotating neutron stars that have a core filled with neutron and proton superfluids, taking account of entrainment effects between the superfluids, and showed that the ordinary fluid-like $`f`$\- modes are not strongly affected by the superfluidity in the core. In the same paper, Lindblom & Mendell (1994) suggested the existence of superfluid-like acoustic modes, which they called $`s_0`$-modes, but they could not find the $`s_0`$-modes numerically. First numerical examples of the superfluid-like acoustic modes of neutron stars with the superfluid core were obtained by Lee (1995), who also showed that $`g`$-modes are not propagating in the core filled with superfluids. Andersson & Comer (2001) discussed the dynamics of superfluid neutron star cores examining local dispersion relations in detail, and they confirmed that $`g`$-modes are excluded from the superfluid core (Lee 1995). Applying their argument to the $`r`$-modes, they suggested that the $`r`$-modes of a superfluid neutron star core are split into two distinct families, that is, ordinary fluid-like $`r`$-modes and superfluid-like $`r`$-modes. Most recently, Prix & Rieutord (2002) have computed in great detail acoustic modes of non-rotating neutron stars composed of neutron and proton superfluids for a wide range of the entrainment parameter. Although their results are basically consistent with the results obtained by Lee (1995) and Andersson & Comer (2001), they have also shown that the assumption that $`|\delta 𝒗_p\delta 𝒗_n|0`$ for the ordinary fluid-like acoustic modes and $`|\rho _p\delta 𝒗_p+\rho _n\delta 𝒗_n|0`$ for the superfluid-like acoustic modes is not necessarily helpful to distinguish the two families, where $`\delta 𝒗_p`$ and $`\delta 𝒗_n`$ are the Eulerian velocity perturbations of the proton and neutron superfluids, the mass densities of which are denoted as $`\rho _p`$ and $`\rho _n`$. We also note that Sedrakian & Wasserman (2000) examined rotational effects on the fundamental modes in uniform density superfluid stars, applying the so-called tensor virial technique to normal modes. As for relativistic non-rotating superfluid stars, Comer, Langlois, & Lin (1999) calculated superfluid-like acoustic modes and showed that the basic features of the superfluid modes are quiet similar to those of the corresponding Newtonian superfluid modes. More recently, Andersson, Comer & Langlois (2002) have found avoided crossings between the superfluid- and the ordinary fluid-like acoustic modes in a relativistic superfluid neutron star as the entrainment is varied.

To give an answer to the question whether the $`r`$-mode instability will survive the damping due to mutual friction in the superfluid core of cold neutron stars, Lindblom & Mendell (2000) computed ordinary fluid-like $`r`$-modes, taking account of the entrainment effects in the core. They showed that the mutual friction damping in the core is not effective to suppress the instability except for a few domains of the entrainment parameter $`\eta `$. Recently, Lee & Yoshida (2002) have studied the $`r`$-modes in neutron stars with the superfluid core with a different method, confirming most of the results by Lindblom & Mendell (2000). Lee & Yoshida (2002) have also found that, as suggested by Andersson & Comer (2001), the $`r`$-modes are split into ordinary fluid-like $`r`$-modes and superfluid-like $`r`$-modes, and that the instability caused by the superfluid-like $`r`$-modes is extremely weak and easily damped by dissipation processes in the interior.

The purpose of this paper is to study modal properties of inertial modes in rotating neutron stars containing superfluids in the core. As shown by Lockitch & Friedman (1999), the existence of inertial modes in the lowest order of angular frequency $`\mathrm{\Omega }`$ of rotation depends on the structure of the so-called zero-frequency subspace of eigensolutions to linearized hydrodynamic equations. Naively speaking, we may say that, for inertial modes to exist in the limit of $`\mathrm{\Omega }0`$, $`g`$-modes of a non-rotating star with $`\mathrm{\Omega }=0`$ need to be degenerate to the zero-frequency subspace. Investigating the zero-frequency subspace of pulsation equations for the superfluid core, Andersson & Comer (2001) showed the existence of the zero-frequency $`g`$-modes, and hence the existence of inertial modes (see, also, Comer 2002). In this paper, we calculate global inertial mode oscillations in neutron stars with a core filled with neutron and proton superfluids, taking account of the entrainment effects between the two superfluids. For the dynamics of superfluids in the core, we employ a two-fluid formalism similar to that given in Lindblom & Mendel (1994). To obtain the inertial modes in the limit of $`\mathrm{\Omega }0`$, we generalize the method devised by Lockitch & Friedman (1999) to stars composed of a superfluid core and a normal fluid envelope. In §2 we present the basic equations employed in this paper for the dynamics of superfluids in the core, and in §3 dissipation processes considered in this paper are described. §4 gives numerical results, and §5 and §6 are for discussions and conclusions.

## 2 oscillation Equations

We derive basic equations governing superfluid motions in the neutron star core in the Newtonian dynamics, assuming uniform rotation of the star. The core is assumed to be filled with neutron and proton superfluids and a normal fluid of electron. We also assume perfect charge neutrality between the protons and electrons because the plasma frequency is much higher than the oscillation frequencies considered in this paper (see, e.g., Mendell 1991a). Since the transition temperatures $`T_c10^9`$ K to neutron and proton superfluids are much higher than the interior temperatures of old neutron stars (see, e.g., Epstein 1988), we assume that all the neutrons and protons in the core are in superfluid states and the excited normal fluid components of the fluids can be ignored.

In this study, we assume that the neutron and proton superfluids and the electron normal fluid in the equilibrium state are in the same rotational motion with the angular velocity $`\mathrm{\Omega }`$ around the axis of rotation. The velocity of the fluid $`v^a`$ in the equilibrium state is therefore given by $`v^a=\mathrm{\Omega }\phi ^a`$, where $`\phi ^a`$ is the rotational Killing vector. In a perturbed state, however, the neutron and the proton superfluids can move differently from each other in the core, obeying their own governing equations. In order to describe the perturbations in the superfluid core, thus, we can choose the perturbed neutron density $`\delta \rho _n`$, the perturbed proton density $`\delta \rho _p`$, the perturbed velocity of the neutron $`\delta v_n^a`$, and the perturbed velocity of the proton $`\delta v_p^a`$ as dynamical variables. Since the equilibrium star is axisymmetric about the rotation axis, we assume that the time and azimuthal dependence of the perturbations is given by $`\mathrm{exp}(i\sigma t+im\phi )`$ with $`\sigma `$ being the oscillation frequency observed in an inertial frame, and $`m`$ is an integer representing the azimuthal wave-number. The basic perturbation equations employed in this paper for the neutron and proton superfluids in rotating neutron stars are essentially the same as those given in Mendell (1991a) and Lindblom & Mendell (1994) and are given as follows:

$$i\omega \delta \rho _n+_a(\rho _n\delta V_n^a)=0,$$
(1)
$$i\omega \delta \rho _p+_a(\rho _p\delta V_p^a)=0,$$
(2)
$`i\omega \delta v_n^a+2\mathrm{\Omega }\delta v_n^b_b\phi ^a+^a(\delta \mu _n+\delta \mathrm{\Phi })2\mathrm{\Omega }{\displaystyle \frac{\rho _{np}}{\rho _n}}(\delta v_p^b\delta v_n^b)^a\phi _b=0,`$ (3)
$`i\omega \left(\delta v_p^a+{\displaystyle \frac{m_e}{m_p}}\delta v_e^a\right)+2\mathrm{\Omega }\left(\delta v_p^b+{\displaystyle \frac{m_e}{m_p}}\delta v_e^b\right)_b\phi ^a`$
$`+^a\left(\delta \mu _p+{\displaystyle \frac{m_e}{m_p}}\delta \mu _e+\left(1+{\displaystyle \frac{m_e}{m_p}}\right)\delta \mathrm{\Phi }\right)2\mathrm{\Omega }{\displaystyle \frac{\rho _{np}}{\rho _n}}(\delta v_n^b\delta v_p^b)^a\phi _b=0,`$ (4)
where $`\omega `$ denotes the corotating frequency, defined by $`\omega =\sigma +m\mathrm{\Omega }`$, $`_a`$ means the covariant derivative in flat three dimensional space, and $`\delta Q`$ stands for the Eulerian change in the physical quantity $`Q`$. Here, $`m_p`$ and $`m_e`$ are the proton and the electron masses. In the mass conservation equations (1) and (2), $`\rho _n\delta V_n^a`$ and $`\rho _p\delta V_p^a`$ mean the mass current vector of the neutron and of the proton, respectively, and can be expanded in terms of the perturbations of the superfluid velocities as follows:

$$\rho _n\delta V_n^a=\rho _{nn}\delta v_n^a+\rho _{np}\delta v_p^a,$$
(5)
$$\rho _p\delta V_p^a=\rho _{pn}\delta v_n^a+\rho _{pp}\delta v_p^a,$$
(6)
where the relationships $`\rho _n=\rho _{nn}+\rho _{np}`$, $`\rho _p=\rho _{pn}+\rho _{pp}`$, and $`\rho _{np}=\rho _{pn}`$ must be satisfied due to the Galilean invariance of equations (1) and (2). The perturbed superfluid velocities are then given in terms of $`\delta V_n^a`$ and $`\delta V_p^a`$ as

$$\delta v_n^a=\frac{\rho _{11}}{\rho _n}\delta V_n^a+\frac{\rho _{12}}{\rho _n}\delta V_p^a,$$
(7)
$$\delta v_p^a=\frac{\rho _{21}}{\rho _p}\delta V_n^a+\frac{\rho _{22}}{\rho _p}\delta V_p^a,$$
(8)
where

$$\rho _{11}=\frac{\rho _{pp}\rho _n^2}{\stackrel{~}{\rho }^2},\rho _{22}=\frac{\rho _{nn}\rho _p^2}{\stackrel{~}{\rho }^2},\rho _{12}=\rho _{21}=\frac{\rho _{np}\rho _n\rho _p}{\stackrel{~}{\rho }^2},$$
(9)
and

$$\stackrel{~}{\rho }^2=\rho _{nn}\rho _{pp}\rho _{np}\rho _{pn}.$$
(10)
Note that $`\rho _{11}+\rho _{12}=\rho _n`$ and $`\rho _{22}+\rho _{21}=\rho _p`$. In Euler equations (3) and (4), $`\delta v_e^a`$ and $`\delta \mathrm{\Phi }`$ are the perturbed velocity of the electron and the perturbed gravitational potential, respectively, and $`\delta \mu _n`$, $`\delta \mu _p`$, and $`\delta \mu _e`$ stand for the perturbed chemical potentials per unit mass for the neutron, proton and electron, respectively. Because of the assumption of the perfect charge neutrality in the proton and electron plasma, we may have

$$\delta V_p^a=\delta v_e^a,\mathrm{and}\delta \rho _p/m_p=\delta \rho _e/m_e,$$
(11)
where $`\rho _e`$ is the electron density. The perturbed gravitational potential is determined by the Poisson equation

$$^2\delta \mathrm{\Phi }=4\pi G\delta \rho ,$$
(12)
where $`\rho =\rho _n+\rho _p+\rho _e`$ and $`G`$ is the gravitational constant.

Using a variant of Gibbs-Duhem relation, the pressure perturbation is given by

$$\delta p=\rho _n\delta \mu _n+\rho _p\delta \stackrel{~}{\mu },$$
(13)
where $`\delta \stackrel{~}{\mu }=\delta \mu _p+\frac{m_e}{m_p}\delta \mu _e`$. Here, we have ignored the entropy carried by the electron normal fluid. This treatment is justified for a star whose temperature is low enough. Note also that the superfluids carry no entropy.

To obtain a relation between the densities $`\rho _n`$ and $`\rho _p`$ and the chemical potentials $`\mu _n`$ and $`\stackrel{~}{\mu }`$, we begin with writing the energy density $`e`$ as

$$e=e(\rho _n,\rho _p),$$
(14)
with which the chemical potentials are defined as

$$\mu _n(\rho _n,\rho _p)=\left(e/\rho _n\right)_{\rho _p},\mu _p(\rho _n,\rho _p)=\left(e/\rho _p\right)_{\rho _n}.$$
(15)
If the chemical potential of the electron is given by

$$\mu _e(\rho _e)=c^2\sqrt{1+(3\pi ^2\mathrm{}^3\rho _e/m_e^4c^3)^{2/3}},$$
(16)
we have, due to charge neutrality condition $`\rho _e=\rho _pm_e/m_p`$,

$$\left(\begin{array}{c}\delta \mu _n\\ \delta \stackrel{~}{\mu }\end{array}\right)=\left(\begin{array}{cc}𝒫_{11}& 𝒫_{12}\\ 𝒫_{21}& 𝒫_{22}\end{array}\right)\left(\begin{array}{c}\delta \rho _n\\ \delta \rho _p\end{array}\right),$$
(17)
where

$$𝒫_{11}=\left(\frac{\mu _n}{\rho _n}\right)_{\rho _p},𝒫_{12}=\left(\frac{\mu _n}{\rho _p}\right)_{\rho _n}=\left(\frac{\stackrel{~}{\mu }}{\rho _n}\right)_{\rho _p}=𝒫_{21},𝒫_{22}=\left(\frac{\stackrel{~}{\mu }}{\rho _p}\right)_{\rho _n}.$$
(18)
We write the inverse of equation (17) as

$$\left(\begin{array}{c}\delta \rho _n\\ \delta \rho _p\end{array}\right)=\left(\begin{array}{cc}𝒬_{11}& 𝒬_{12}\\ 𝒬_{21}& 𝒬_{22}\end{array}\right)\left(\begin{array}{c}\delta \mu _n\\ \delta \stackrel{~}{\mu }\end{array}\right),$$
(19)
where

$$\left(\begin{array}{cc}𝒬_{11}& 𝒬_{12}\\ 𝒬_{21}& 𝒬_{22}\end{array}\right)=\left(\begin{array}{cc}𝒫_{11}& 𝒫_{12}\\ 𝒫_{21}& 𝒫_{22}\end{array}\right)^1.$$
(20)
The curl of the Euler equations, which is the vorticity equation, is useful for perturbation analysis of rotating stars, and is given in terms of $`\delta V_n^a`$ and $`\delta V_p^a`$ by

$`ϵ^{abc}_b\left[\left(i\omega {\displaystyle \frac{\rho _{11}}{\rho _n}}g_{cd}+2\mathrm{\Omega }_d\phi _c\right)\delta V_n^d+i\omega {\displaystyle \frac{\rho _{12}}{\rho _n}}g_{cd}\delta V_p^d\right]=0,`$ (21)
$`ϵ^{abc}_b\left[\left(i\omega \left({\displaystyle \frac{\rho _{22}}{\rho _p}}+{\displaystyle \frac{m_e}{m_p}}\right)g_{cd}+2\mathrm{\Omega }\left(1+{\displaystyle \frac{m_e}{m_p}}\right)_d\phi _c\right)\delta V_p^d+i\omega {\displaystyle \frac{\rho _{12}}{\rho _p}}g_{cd}\delta V_n^d\right]=0,`$ (22)
where $`g_{ab}`$ and $`ϵ^{abc}`$ mean the metric tensor and the Levi-Civita tensor in the flat space.

Although we have made no assumptions for rotation velocity up to now, let us restrict our consideration to slowly rotating stars in order to obtain inertial modes to the lowest order in the stellar rotation frequency $`\mathrm{\Omega }`$. In the limit of slow rotation, then, we assume perturbations to obey the following ordering in $`\mathrm{\Omega }`$:

$$\delta V_n^a=O(1),\delta V_p^a=O(1),\omega /\mathrm{\Omega }=O(1),\delta \rho _n=O(\mathrm{\Omega }),\delta \rho _p=O(\mathrm{\Omega }),\mathrm{and}\delta \mathrm{\Phi }=O(\mathrm{\Omega }).$$
(23)
Examining the zero-frequency subspace of the perturbation in a non-rotating superfluid core, Andersson & Comer (2001) showed that solutions obeying the ordering (23) may be allowed in the equations describing hydrodynamics in the superfluid core. The existence of solutions obeying (23) is similar to that in a barotropic normal fluid star (Lockitch & Friedman 1999). Those solutions can be interpreted as infinitely degenerate $`g`$-modes of zero-frequency if we consider no rotational effect. It is convenient to introduce dimensionless frequency $`\kappa _0`$ in the limit of $`\mathrm{\Omega }0`$, which is defined by

$$\kappa _0=\underset{\mathrm{\Omega }0}{lim}\omega /\mathrm{\Omega }.$$
(24)
Then, basic equations describing the lowest order inertial modes are given by

$$_a(\rho _n\delta V_n^a)=0,_a(\rho _p\delta V_p^a)=0,$$
(25)
$`q_n^a=ϵ^{abc}_b\left[\left({\displaystyle \frac{\rho _{11}}{\rho _n}}g_{cd}2i\kappa _0^1_d\phi _c\right)\delta V_n^d+{\displaystyle \frac{\rho _{12}}{\rho _n}}g_{cd}\delta V_p^d\right]=0,`$ (26)
$`q_p^a=ϵ^{abc}_b\left[\left(\left({\displaystyle \frac{\rho _{22}}{\rho _p}}+{\displaystyle \frac{m_e}{m_p}}\right)g_{cd}2i\kappa _0^1\left(1+{\displaystyle \frac{m_e}{m_p}}\right)_d\phi _c\right)\delta V_p^d+{\displaystyle \frac{\rho _{12}}{\rho _p}}g_{cd}\delta V_n^d\right]=0.`$ (27)
Note that $`q_n^a`$ ($`q_p^a`$) has only two independent components because $`_aq_n^a=0`$ ($`_aq_p^a`$=0) is automatically satisfied.

¿From similar consideration to that for the superfluid core, basic equations for pulsations in the normal fluid envelope are given by

$$_a(\rho \delta v^a)=0,$$
(28)
$`q^a=ϵ^{abc}_b\left[\left(g_{cd}2i\kappa _0^1_d\phi _c\right)\delta v^d\right]=0,`$ (29)
where $`\delta v^a`$ denotes the perturbation of the fluid velocity.

The perturbations in a uniformly rotating star are expanded in terms of spherical harmonic functions with different $`l`$’s for a given $`m`$ (e.g., Lee & Saio 1986). The vector $`\delta V_X^a`$ is therefore given by

$$\delta V_X^r=r\underset{l|m|}{}S_l^X(r)Y_l^m(\theta ,\phi )e^{i\sigma t},$$
(30)
$$\delta V_X^\theta =\underset{l|m|}{}\left[H_l^X(r)\frac{}{\theta }Y_l^m(\theta ,\phi )iT_l^{}^X(r)\frac{1}{\mathrm{sin}\theta }\frac{}{\phi }Y_l^{}^m(\theta ,\phi )\right]e^{i\sigma t},$$
(31)
$$\delta V_X^\phi =\frac{1}{\mathrm{sin}^2\theta }\underset{l|m|}{}\left[H_l^X(r)\frac{}{\phi }Y_l^m(\theta ,\phi )+iT_l^{}^X(r)\mathrm{sin}\theta \frac{}{\theta }Y_l^{}^m(\theta ,\phi )\right]e^{i\sigma t},$$
(32)
where $`l=|m|+2(k1)`$ and $`l^{}=l+1`$ for even modes, and $`l=|m|+2k1`$ and $`l^{}=l1`$ for odd modes, where $`k=1,2,3,\mathrm{}`$. Here, we have employed the constituent index $`X`$ which can be either $`n`$, $`p`$ or $`o`$, where $`o`$ means the velocity field in the normal fluid envelope, and the normal spherical polar coordinate $`(r,\theta ,\phi )`$ is used. In this paper, the even and odd modes are respectively characterized by symmetry and antisymmetry of their eigenfunctions with respect to the equatorial plan. Substituting these expansions into equations (25) through (27), we obtain oscillation equations given as a set of simultaneous linear ordinary differential equations of the expansion coefficients (see Appendix), which is to be integrated in the superfluid core. The oscillation equations solved in the normal fluid envelope are similar to those given in Lockitch & Friedman (1999).

To obtain a unique solution of an oscillation mode, solutions in the superfluid core and in the normal fluid envelope are matched at the interface between the two domains by imposing jump conditions given by

$$\delta V_n^r=\delta V_p^r=\delta V_o^r.$$
(33)
The boundary conditions at the stellar center is the regularity condition of the perturbations $`\delta V_n^a`$ and $`\delta V_p^a`$. The boundary condition at the stellar surface is $`\mathrm{\Delta }p=0`$, where $`\mathrm{\Delta }p`$ is the Lagrangian change in the pressure.

For numerical computation, oscillation equations of a finite dimension are obtained by disregarding the terms with $`l`$ larger than $`l_{max}=|m|+2k_{max}1`$ in the expansions of perturbations such as given by (30) to (32), where $`l_{max}`$ is determined so that the eigenfrequency and the eigenfunctions are well converged as $`l_{max}`$ increases. We solve the oscillation equations of a finite dimension as an eigenvalue problem with the scaled oscillation frequency $`\kappa _0`$ using a Henyey type relaxation method (see, e.g., Unno et al. 1989).

## 3 Dissipations

In the previous section, we consider no dissipation effects. Thus, the conserved energy $`E`$ of oscillations exits and is given by (Mendell 1991b)

$`E=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{V_S}}\left(\rho \delta v_a^{}\delta v^a+{\displaystyle \frac{\stackrel{~}{\rho }^2}{\rho }}\delta w_a^{}\delta w^a+{\displaystyle \underset{i,j}{}}𝒫_{ij}\delta \rho _i^{}\delta \rho _j{\displaystyle \frac{1}{4\pi G}}_a\delta \mathrm{\Phi }^{}^a\delta \mathrm{\Phi }\right)\sqrt{g}d^3x`$ (34)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{V_N}}\left(\rho \delta v_a^{}\delta v^a+\left({\displaystyle \frac{\delta p}{\rho }}+\delta \mathrm{\Phi }\right)\delta \rho ^{}\right)\sqrt{g}d^3x`$
where $`g`$ is the determinant of the metric, and asterisk $`(^{})`$ denotes the complex conjugate of the quantity. Here, $`_{V_S}d^3x`$ and $`_{V_N}d^3x`$ means the integration over the superfluid core and the normal fluid envelope, respectively. In this equation, we have employed new variables, defined by

$`\delta v^a=(\rho _n\delta v_n^a+\rho _p\delta v_p^a)/\rho =(\rho _n\delta V_n^a+\rho _p\delta V_p^a)/\rho ,`$
$`\delta w^a=\delta v_n^a\delta v_p^a={\displaystyle \frac{\rho _n\rho _p}{\stackrel{~}{\rho }^2}}\left(\delta V_n^a\delta V_p^a\right),`$ (35)
$`\delta \rho _1=\delta \rho _n,\delta \rho _2=\delta \rho _p.`$
Note that $`\rho \delta v^a`$ is interpreted as the total momentum current vector of perturbation. In reality, however, there are a lot of dissipation mechanisms in real neutron stars. Thus, the oscillation energy $`E`$ is not conserved in real neutron stars. In other words, the amplitude of an oscillation mode should damp or grow depending on considered dissipation mechanisms. The stability of an oscillation mode of a star is then determined by summing up all contributions from various dissipation mechanisms. In this study, we consider the contributions from gravitational radiation reaction, viscous processes, and mutual friction in the superfluid core. Then, the averaged energy loss (or gain) rate $`dE/dt`$ of a normal mode in a rotating neutron star may be given by

$`{\displaystyle \frac{dE}{dt}}=`$ $`\sigma \omega {\displaystyle \underset{l^{}|m|}{}}N_l\sigma ^{2l}\left|\delta J_{l^{}m}\right|^22{\displaystyle _{V_S+V_N}}\eta _s\delta \sigma ^{ab}\delta \sigma _{ab}\sqrt{g}d^3x`$
$`2\mathrm{\Omega }{\displaystyle _{V_S}}\rho _nB_n\left({\displaystyle \frac{\stackrel{~}{\rho }^2}{\rho _n\rho _p}}\right)^2(g_{ab}z_az_b)\delta w^a\delta w^b\sqrt{g}d^3x`$
$`=`$ $`{\displaystyle \underset{l^{}|m|}{}}\left({\displaystyle \frac{dE}{dt}}\right)_{GJ,l^{}}+\left({\displaystyle \frac{dE}{dt}}\right)_S+\left({\displaystyle \frac{dE}{dt}}\right)_{MF},`$ (36)
where $`z^a`$ is the unit vector that parallels the rotation axis, defined by $`z_a=ϵ_{abc}^b\phi ^c/2`$.

The terms $`(dE/dt)_{GJ,l}`$ on the right hand side of equation (36) denote the energy loss (or gain) rates due to gravitational radiation reaction associated with the mass current multipole moment $`\delta J_{lm}`$, where

$$\delta J_{lm}=\frac{2}{c(l+1)}_{V_S+V_N}r^{l+1}\rho \delta v_aϵ^{abc}_bY_l^m_cr\sqrt{g}d^3x$$
(37)
and

$$N_l=\frac{4\pi G}{c^{2l+1}}\frac{(l+1)(l+2)}{l(l1)[(2l+1)!!]^2},$$
(38)
and $`c`$ is the velocity of light (Thorne 1980, Lindblom et al. 1998).

The term $`(dE/dt)_S`$ is the energy dissipation rate due to shear viscosity, and $`\eta _s`$ is the shear viscosity coefficient, and $`\delta \sigma _{ab}`$ is the traceless rate of strain tensor for the perturbed velocity field, which is defined by

$$\delta \sigma _{ab}=\frac{1}{2}\left(_au_b+_bu_a\frac{2}{3}g_{ab}_cu^c\right),$$
(39)
for the velocity field $`u^a`$. In this paper, we ignore the contribution from the bulk viscosity, which is important only for newly born hot neutron stars without superfluids in the core. The shear viscosity coefficient we use in the superfluid core is

$$\eta _S=6.0\times 10^{18}\left(\frac{\rho }{10^{15}\mathrm{g}/\mathrm{cm}^3}\right)^2\left(\frac{10^9K}{T}\right)^2\mathrm{g}/\mathrm{cm}\mathrm{s}$$
(40)
(Cutler & Lindblom 1987, Flowers & Itoh 1979), and that in the normal fluid envelope is given by

$$\eta _S=1.95\times 10^{18}\left(\frac{\rho }{10^{15}\mathrm{g}/\mathrm{cm}^3}\right)^{9/4}\left(\frac{10^9K}{T}\right)^2\mathrm{g}/\mathrm{cm}\mathrm{s}$$
(41)
(Cutler & Lindblom 1987, Flowers & Itoh 1979). The rate of strain $`\delta \sigma _{ab}`$ is evaluated by using $`\delta v_e^a`$ in the superfluid core, and $`\delta v^a`$ in the normal fluid envelope (Lindblom & Mendell 2000).

The term $`(dE/dt)_{MF}`$ is the energy loss rate due to the mutual friction in the superfluid core, and the dimensionless coefficient $`B_n`$ is given by (Mendell 1991b)

$$B_n=0.011\times \frac{\rho _p}{\rho _n}\left(\frac{\rho _{pp}}{\rho _p}\right)^{1/2}\left(\frac{\rho _{pn}}{\rho _{pp}}\right)^2\left(\frac{\rho _p}{10^{14}\mathrm{g}/\mathrm{cm}^3}\right)^{1/6}.$$
(42)
The mutual friction is a dissipation mechanism inherent to a rotating system of superfluids, and it is caused by scattering of normal fluid particles off the vortices in the superfluids. Since we have assumed perfect charge neutrality between the electrons and protons, we consider scattering between the normal electrons and vortices of the neutron superfluid (Mendell 1991b, see also Alpar, Langer, & Sauls, 1984).

As is indicated by the first term on the right hand side of equation (36), if the normal mode has an oscillation frequency that satisfies $`\sigma \omega >0`$, the oscillation energy $`E`$ in the corotating frame, in the absence of other damping mechanisms, increases as a result of gravitational wave radiation, indicating instability of the mode (Friedman & Schutz 1978).

The damping (or growth) time-scale $`\tau `$ of a normal mode may be given by

$$\frac{1}{\tau }=\frac{1}{2E}\left(\frac{dE}{dt}\right)=\underset{l^{}>|m|}{}\frac{1}{\tau _{GJ,l^{}}}+\frac{1}{\tau _S}+\frac{1}{\tau _{MF}},$$
(43)
where $`\tau _i=2E/(dE/dt)_i`$. For the inertial modes, it is convenient to derive an extrapolation formula of the time-scale $`\tau `$ given as a function of $`\mathrm{\Omega }`$ and the interior temperature $`T`$ (e.g., Lindblom et al. 1998, Lindblom & Mendell 2000):

$$\frac{1}{\tau (\mathrm{\Omega },T)}=\underset{l^{}|m|}{}\frac{1}{\tau _{GJ,l^{}}^0}\left(\frac{\mathrm{\Omega }^2}{\pi G\overline{\rho }}\right)^l^{}+\frac{1}{\tau _S^0}\left(\frac{10^9K}{T}\right)^2+\frac{1}{\tau _{MF}^0}\left(\frac{\mathrm{\Omega }^2}{\pi G\overline{\rho }}\right)^{1/2},$$
(44)
where $`\overline{\rho }=M/(4\pi R^3/3)`$. Here, $`M`$ and $`R`$ are the mass and radius of the star. The quantities $`\tau _{GJ,l^{}}^0`$, $`\tau _S^0`$, and $`\tau _{MF}^0`$ are assumed to be independent on $`\mathrm{\Omega }`$ and $`T`$.

## 4 Numerical Results

Following Lindblom & Mendell (2000), we employ a polytropic model of index $`N=1`$ with the mass $`M=1.4M_{}`$ and the radius $`R=12.57`$km as a background model for modal analysis. The model is divided into a superfluid core and a normal fluid envelope, the interface of which is set at $`\rho =\rho _s=2.8\times 10^{14}`$g/cm<sup>3</sup>. In the normal fluid envelope, we use the polytropic equation of state given by $`p=K\rho ^2`$ both for the equilibrium structure and for the oscillation equations. The core is assumed to be filled with neutron and proton superfluids and a normal fluid of electron, for which an equation of state, labeled A18+$`\delta v`$+UIX (Akmal, Pandharipande, & Ravenhall 1998), is used to give the thermodynamical quantities used for the oscillation equations. For the mass density coefficients $`\rho _{nn}`$, $`\rho _{pp}`$, and $`\rho _{np}`$ in the core, we employ an empirical relation given by (see Lindblom & Mendell 2000)

$$\rho _p/\rho 0.031+8.8\times 10^{17}\rho ,$$
(45)
and a formula given by

$$\rho _{np}=\eta \rho _n,$$
(46)
where $`\eta `$ is a parameter of order of $`0.04`$ (Borumand, Joynt, & Kluźniak 1996). In this paper we call $`\eta `$ the entrainment parameter. In this study, only the modes having $`m=2`$ are examined because they are the most significant for the gravitational radiation driven instability. For the range of the entrainment parameter $`\eta `$, we choose $`0.01\eta 0.04`$. Although negative values of $`\eta `$ might not be suitable for physical models of neutron stars, we include in our calculations negative $`\eta `$ to demonstrate clearly the $`\eta `$ dependence of the inertial modes around $`\eta 0`$.

We find that inertial modes in the superfluid core are split into two families, which we call ordinary fluid inertial modes ($`i^o`$-mode), and superfluid inertial modes ($`i^s`$-mode). This mode splitting is common in other oscillation modes in the superfluid core, because the dynamical degree of freedom in the superfluid core simply doubles (see, e.g., Epstein 1988; Lindblom & Mendell 1994; Lee 1995; Andersson & Comer 2001; Prix & Rieutord 2002). As a subclass of inertial modes, the $`r`$-modes of neutron stars with a superfluid core are also split into ordinary fluid-like $`r`$-modes and superfluid-like $`r`$-modes, which we respectively call $`r^o`$-modes and $`r^s`$-modes as in Lee & Yoshida (2002). Since basic properties of the eigenfunctions $`\delta v^a`$ (or $`\delta w^a`$) of an inertial mode are quite similar to the eigenfunctions $`\delta v^a`$ of the corresponding inertial mode in the normal fluid star without core superfluidity, we employ the same classification scheme as that employed by Lockitch & Friedman (1999) and Yoshida & Lee (2000a) for inertial modes in the normal fluid star. In this scheme, inertial modes are classified in terms of the two quantum numbers $`m`$ and $`l_0|m|`$, and the scaled frequency $`\kappa _0lim_{\mathrm{\Omega }0}\omega /\mathrm{\Omega }`$. For given $`m`$ and $`l_0|m|`$, the number of inertial modes of one kind is $`l_0|m|`$. Even parity inertial modes have even $`l_0|m|`$, and odd parity inertial modes have odd $`l_0|m|`$. Note that the $`r`$-modes are classified as inertial modes with $`l_0|m|=1`$ (e.g., Yoshida & Lee 2000a). It is interesting to note that, although there exist three quantum numbers to classify the normal modes of non-rotating stars, a third quantum number may be missing for inertial modes. Usually, the ratio $`E_S/E`$ of the oscillation energy $`E_S`$ in the core to the total oscillation energy $`E`$ is a good indicator for a distinction between $`i^s`$\- and $`i^o`$-modes, since we find $`E_S/E<0.5`$ for the $`i^o`$-modes and $`E_S/E>0.9`$ for the $`i^o`$-modes for the models composed of a superfluid core and a normal fluid envelope.

Let us first consider the case of $`\eta =0`$, that is, the case of no entrainment. In Table 1, $`\kappa _0`$’s for $`i^o`$\- and $`i^s`$-modes of low $`l_0|m|`$ are tabulated together with $`\kappa _0`$’s for the corresponding inertial modes in the $`N=1`$ polytrope model without superfluidity in the core, which we simply call $`i`$-mode (see, Lockitch & Friedman 1999; Yoshida & Lee 2000a). Our values of $`\kappa _0`$ for the $`i^o`$\- and $`i^s`$-modes are in good agreement with those computed by Lee & Yoshida (2002) with a different method. Note that we have shown five digits for numerical values of $`\kappa _0`$ to indicate that our numerical results include at most a few percent numerical error. We can see that $`\kappa _0`$’s of paired $`i^o`$\- and $`i^s`$-modes have quite similar values to that of the corresponding inertial mode in the $`N=1`$ polytropic normal fluid model. In order to demonstrate model properties of the inertial modes, we exhibit the expansion coefficients $`S_l^v`$, $`T_l^{}^v`$, $`S_l^w`$, and $`T_l^{}^w`$ of the eigenfunctions $`\delta v^a`$ and $`\delta w^a`$, defined by

$$S_l^v=(\rho _nS_l^n+\rho _pS_l^p)/\rho ,T_l^{}^v=(\rho _nT_l^{}^n+\rho _pT_l^{}^p)/\rho ,$$
(47)
$$S_l^w=\rho _n\rho _p(S_l^nS_l^p)/\stackrel{~}{\rho }^2,T_l^{}^w=\rho _n\rho _p(T_l^{}^nT_l^{}^p)/\stackrel{~}{\rho }^2.$$
(48)
Note that $`S_l^w`$ and $`T_l^{}^w`$ vanish in the normal fluid envelope because neutrons and protons co-move. In Figures 1 through 4, the first three expansion coefficients $`S_l^v`$, $`T_l^{}^v`$, $`S_l^w`$, and $`T_l^{}^w`$ are shown for the $`i^o`$\- and $`i^s`$-modes of $`l_0|m|=2`$ which have $`\kappa _0=0.5685`$ and $`\kappa _0=0.5115`$, respectively. Here, the amplitude normalization given by $`T_{m+1}^v=1`$ has been used. It is noted that $`\delta v^a`$ of the $`i^o`$-modes and $`\delta w^a`$ of the $`i^s`$-modes resemble $`\delta v^a`$ of the corresponding inertial modes of the normal fluid star model (compare the figures with Figure 1 of Yoshida & Lee 2000a). It is also found that although the amplitudes of $`\delta v^a`$ and $`\delta w^a`$ are similar to each other for the $`i^o`$-modes, the amplitude of $`\delta w^a`$ dominates that of $`\delta v^a`$ for the $`i^s`$-modes, which indicates that the two superfluids in the core basically counter-move for the $`i^s`$-modes.

In Figures 5 to 8, $`\kappa _0`$’s of inertial modes of low $`l_0|m|`$ are plotted as functions of the entrainment parameter $`\eta `$, where Figure 5 is for even parity modes with positive $`\kappa _0`$, Figure 6 for even parity modes with negative $`\kappa _0`$, Figure 7 for odd parity modes with positive $`\kappa _0`$, and Figure 8 for odd parity modes with negative $`\kappa _0`$. In Figure 7, we have also given $`\kappa _0`$ of the $`r^s`$-mode, the locus of which as a function of $`\eta `$ is in good agreement with the locus calculated by Lee & Yoshida (2002) with a different method. Note that we cannot obtain $`r^o`$-modes by using our numerical method, the reason for which is not clear (see the next section) <sup>1</sup><sup>1</sup>1Recently, we have generalized our formalism to the case of general relativistic superfluid neutron stars. In that study, our numerical code for Newtonian dynamics have been straightforwardly extended into the one for general relativistic dynamics. With this relativistic version of the code, both $`r^o`$\- and $`r^s`$-modes can be obtained even in almost the Newtonian regime with small values of the relativistic parameter $`M/R`$. Fundamental differences in the two formalisms reside only in the master equations. As shown in equations (49)–(52), the basic equations for neutron and proton fluids in the Newtonian case get to degenerate in lowest order of $`\mathrm{\Omega }`$ when $`r^o`$-modes are considered, while the basic equations in the relativistic case don’t even if we are concerned with $`r^o`$-modes (see, Yoshida & Lee 2003). We therefore guess that the reason why $`r^o`$-modes cannot be obtained in the Newtonian case is this possible degeneracy of the two sets of superfluid equations in Newtonian dynamics.. From Figures 5 through 8, we find that, although $`\kappa _0`$ for the $`i^o`$-modes is only weakly dependent on the parameter $`\eta `$, $`|\kappa _0|`$ for the $`i^s`$-modes increases almost linearly as $`\eta `$ increases. Because of these different dependences of $`\kappa _0`$ on the parameter $`\eta `$, frequent avoided crossings are inevitable between $`i^o`$-modes and $`i^s`$-modes. As shown by Figure 7, for example, the $`i^s`$-mode of $`l_0|m|=3`$ which has $`\kappa _0=1.3865`$ at $`\eta =0`$ suffers three avoided crossings with different $`i^o`$-modes as $`\eta `$ is increased from $`0.01`$ to $`0.04`$. Since $`\kappa _0`$’s of paired $`i^o`$\- and $`i^s`$-modes at $`\eta =0`$ are close to each other, they always make an avoided crossing around $`\eta =0`$, which is also shown by Figures 5 through 8.

In Figures 9 through 12, we display the first three coefficients $`S_l^v`$, $`T_l^{}^v`$, $`S_l^w`$, and $`T_l^{}^w`$ for the $`i^o`$\- and $`i^s`$-modes of $`l_0|m|=2`$, which respectively have $`\kappa _0=0.5403`$ and $`\kappa _0=0.6103`$ at $`\eta =0.01`$. Here, the amplitude normalization given by $`T_{m+1}^v=1`$ has been used. Note that these two inertial modes have $`\kappa _0=0.5685`$ and $`\kappa _0=0.5115`$ at $`\eta =0`$. Comparing these figures with Figures 1 to 4, we observe that the basic properties of the eigenfunctions of the inertial modes do not depend on the entrainment parameter $`\eta `$, except at avoided crossings.

In Table 2, we tabulate $`\kappa _0`$ and the scaled damping (or growth) timescales $`\tau _i^0`$ for the $`i^o`$\- and $`i^s`$-modes with positive $`\kappa _0`$ for $`\eta =0`$, $`0.02`$, and $`0.04`$. We find that the growth timescales $`\tau _{GJ,l^{}}^0`$ for all the $`i^o`$\- and $`i^s`$-modes in the table are longer than that of the $`r`$-mode of the normal fluid neutron star, which means that the instability of the $`i^o`$\- and $`i^s`$-modes caused by the current multipole radiation is not strong. This result is the same as that for the inertial modes in a normal fluid star. It is interesting to note that the two $`i^s`$-modes in the table having $`\kappa _0=2.1661`$ and $`\kappa _0=2.4083`$ at $`\eta =0.04`$ are stable against gravitational radiation reaction since they have positive $`\tau _{GJ,l^{}}^0`$. The reason for this can be understood as follows. As mentioned in the last section, the gravitational radiation driven instability sets in when the frequency of the mode satisfies $`\sigma (\sigma m\mathrm{\Omega })<0`$, which leads to the instability criterion given by $`0<\kappa _0<m`$ for positive $`m`$. For inertial modes of normal fluid stars, it is known that all inertial modes having positive $`\kappa _0`$ for $`m>0`$ satisfy the instability criterion (Lindblom & Ipser 1999; Lockitch & Friedman 1999; Yoshida & Lee 2000a). For inertial modes in the superfluid core, on the other hand, since the frequency $`|\kappa _0|`$ of the $`i^s`$-modes monotonically increases as the value of $`\eta `$ is increased, the $`\kappa _0`$ of some positive-frequency $`i^s`$-modes exceeds the limit $`m`$ when $`\eta `$ gets larger than a critical value $`\eta _{crit}`$, the value of which depends on the mode (see, Figures 5 and 7). Thus, some of the positive-frequency $`i^s`$-modes that are unstable satisfying the instability criterion at $`\eta =0`$ become stable invalidating the criterion for $`\eta `$ larger than $`\eta _{crit}`$. As for the effects of mutual friction, we observe that the damping timescales $`\tau _{MF}^0`$ for the inertial modes are much shorter than those due to other dissipation processes unless mutual friction between the superfluids is extremely weak, that is, the parameter $`\eta `$ is extremely small. Note that the amount of dissipation caused by the mutual friction is roughly proportional to $`B_n|\delta w^a|^2`$ where $`B_n\eta ^2`$. Accordingly, if we assume $`\eta 0.04`$ as a typical value for the parameter, the dissipation due to mutual friction in the superfluid core is very effective to damp out the $`i^o`$\- and $`i^s`$-modes.

## 5 Discussion

As mentioned in the previous section, with our method of calculation we fail in obtaining ordinary fluid $`r`$-modes, which are called $`r^o`$-modes in Lee & Yoshida (2002). However, this does not necessarily mean that our basic equations (A2)–(A7) do not contain such $`r`$-mode solutions. To prove this, let us consider purely axial parity solutions, whose angular behavior is given by a purely axial harmonic function having $`l^{}=m`$. In this section, we also assume for simplicity that the densities $`\rho _n`$, $`\rho _p`$, and $`\rho _{12}`$ are constant and do not depend on $`r`$. In virtue of the assumption of purely axial mode, we can reduce our basic equations to the following equations (see, Appendix):

For the superfluid core,

$$\nu m(m+2)J_{m+1}^m\left\{r\frac{d}{dr}T_m^n(m1)T_m^n\right\}=0,$$
(49)
$$\nu m(m+2)J_{m+1}^m\left\{r\frac{d}{dr}T_m^p(m1)T_m^p\right\}=0,$$
(50)
$$\left(1\frac{m\nu }{m(m+1)}\right)T_m^n+\frac{\rho _{12}}{\rho _n}(T_m^pT_m^n)=0,$$
(51)
$$\left(1\frac{m\nu }{m(m+1)}\right)T_m^p+\frac{\rho _{12}}{\rho _p}(T_m^nT_m^p)=0,$$
(52)
For the normal fluid envelope,

$$\nu m(m+2)J_{m+1}^m\left\{r\frac{d}{dr}T_m^o(m1)T_m^o\right\}=0,$$
(53)
$$\left(1\frac{m\nu }{m(m+1)}\right)T_m^o=0,$$
(54)
where $`\nu =2\mathrm{\Omega }/\sigma `$ and $`J_l^m=[(l^2m^2)/(4l^21)]^{1/2}`$. Here, we have ignored the terms due to $`m_e/m_p`$ for simplicity. The condition for the continuity of the total momentum flux at the interface between the superfluid core and the normal fluid envelope is given by

$$T_m^o=(\rho _nT_m^n+\rho _pT_m^p)/\rho \mathrm{at}r=r_s,$$
(55)
where $`r_s`$ is the radius of the superfluid core. Because the factor $`\nu m(m+2)J_{m+1}^m`$ in equations (49), (50) and (53) does not vanish for $`\nu 0`$ and $`m2`$, we can obtain solutions to differential equations (49), (50), and (53) as

$$T_m^n=A_n(r/r_s)^{m1},T_m^p=A_p(r/r_s)^{m1},T_m^o=B(r/r_s)^{m1},$$
(56)
where $`A_n`$, $`A_p`$, and $`B`$ are constants, which are to be determined by boundary condition (55) and constraint equations (51), (52), (54). Substituting the solutions (56) into algebraic equations (51), (52), (54) and (55), we obtain

$$\left(\begin{array}{cc}f\rho _{12}/\rho _n& \rho _{12}/\rho _n\\ \rho _{12}/\rho _p& f\rho _{12}/\rho _p\end{array}\right)\left(\begin{array}{c}A_n\\ A_p\end{array}\right)W\left(\begin{array}{c}A_n\\ A_p\end{array}\right)=0,$$
(57)
$$fB=0,$$
(58)
and

$$B=(\rho _n(r_s)A_n+\rho _p(r_s)A_p)/\rho (r_s),$$
(59)
where $`f=1\nu /(m+1)`$. For non-trivial solutions to exist in the core, we must have $`detW=0`$, which leads to $`f=0`$ or $`f=\rho _{12}(\rho _n+\rho _p)/(\rho _n\rho _p)`$. For the case of $`f=0`$ (i.e. $`\kappa _0=2/(m+1)`$), solutions may be given by

$$A_n=A_p=B.$$
(60)
This solution obviously corresponds to the $`r^o`$-mode since the oscillation frequency $`\kappa _0`$ is exactly the same as that of the $`r^o`$-mode, and the two superfluids flow together in the core. For the case of $`f=\rho _{12}(\rho _n+\rho _p)/(\rho _n\rho _p)`$ (i.e., $`\kappa _02/(m+1)`$), on the other hand, solutions are given by

$$(\rho _n(r_s)A_n+\rho _p(r_s)A_p)/\rho (r_s)=B=0.$$
(61)
This solution corresponds to the $`r^s`$-modes found by Lee & Yoshida (2002), for which the two superfluids counter-move in the core (see, also, Andersson & Comer 2001). It is to be noted that the solution (60) for arbitrary $`\eta `$ and the solution (61) for $`\eta =0`$ can be derived even without the assumption that the densities $`\rho _n`$, $`\rho _p`$, and $`\rho _{12}`$ do not depend on the radial distance $`r`$.

## 6 Conclusion

We have discussed the modal properties of inertial modes of rotating neutron stars with the core filled with neutron and proton superfluids, taking account of entrainment effects between the superfluids. We find that the inertial modes are split into two families, ordinary fluid inertial modes ($`i^o`$-mode) and superfluid inertial modes ($`i^s`$-mode). The two superfluids in the core counter-move for the $`i^s`$-modes. For the $`i^o`$-modes, the dimensionless frequency $`\kappa _0`$ is only weakly dependent on the entrainment parameter $`\eta `$. For the $`i^s`$-modes, on the other hand, $`|\kappa _0|`$ almost linearly increases as the entrainment parameter $`\eta `$ is increased. Because of the different dependences of $`\kappa _0`$ on the parameter $`\eta `$, avoided crossings as functions of $`\eta `$ are common between $`i^o`$\- and $`i^s`$-modes. We also find that some of the $`i^s`$-modes that are unstable against the gravitational radiation reaction at $`\eta =0`$ become stable invalidating the instability criterion when $`\eta `$ becomes larger than some critical $`\eta _{crit}`$, the value of which depends on the mode. The radiation driven instability associated with the current multipole radiation is quite weak for the inertial modes. Since the mutual friction damping in the superfluid core is quite strong for a typical value of $`\eta `$ compared with other damping mechanisms considered in this paper, we conclude that the inertial modes are easily damped by the mutual friction.

Recently, Arras et al. (2002) considered non-linear couplings between an unstable $`r`$-mode and stable inertial modes of rotating normal fluid neutron stars as a mechanism that limits the amplitude growth of the $`r`$-mode, and discussed in the WKB limit the modal properties of the inertial modes, which are important to determine the saturation amplitude of the $`r`$-mode. In their analysis, however, no effects of buoyancy and possible core superfluidity on the inertial modes are correctly included, which substantially affect the modal properties such as the frequency spectra and damping rates (see Yoshida & Lee 2000b). If we consider cold neutron stars in a neutron star binary, inertial modes in the superfluid core of the stars will be excited by the tidal interaction with the companion, and the mutual friction dissipation produced by the excited inertial modes will play an important role in the energy and angular momentum transfer between the orbital motion and the neutron stars, since no $`g`$-modes are propagating in the superfluid core (Lee 1995, Andersson & Comer 2001). In fact, the damping timescale of the inertial modes propagating in the superfluid core may be given by $`\tau \tau _{MF}^0/\overline{\mathrm{\Omega }}`$ sec, and it is for $`\overline{\mathrm{\Omega }}0.1`$ of order of 0.1 to 1 sec, which is extremely short compared with the damping timescales of $`g`$-modes due to viscous processes and emission of gravitational waves (see, e.g., McDermott, Van Horn, & Hansen 1988; Reisenegger & Goldreich 1992). Considering that the orbital period decrease $`d\mathrm{ln}P/dt`$ due to gravitational radiation reactions is proportional to $`D^4`$ (e.g., Shapiro & Teukolsky 1983) but that due to the tidal interaction to $`D^5`$ (e.g., Savonije & Papaloizou 1983) where $`D`$ is the binary separation, the tidal interaction in terms of the inertial modes in the superfluid core becomes important as the binary separation decreases.

Acknowledgments: S.Y. acknowledges financial support from Fundação para a Ciência e a Tecnologia (FCT) through project SAPIENS 36280/99.

## Appendix A Oscillation Equations

For the oscillation equations in the superfluid core, we employ vectors $`𝒔^n`$, $`𝒔^p`$, $`𝒕^n`$, $`𝒕^p`$, $`𝒉^n`$, and $`𝒉^p`$, whose components are given by

$$s_k^n=S_l^n,s_k^p=S_l^p,t_k^n=T_l^{}^n,t_k^p=T_l^{}^p,h_k^n=H_l^n,h_k^p=H_l^p,$$
(62)
where $`l=|m|+2(k1)`$, $`l^{}=l+1`$ for even modes and $`l=|m|+2k1`$, $`l^{}=l1`$ for odd modes, where $`k=1,2,3,\mathrm{}`$. Using these vectors, the perturbed continuity equations (1) and (2) are written as

$`r{\displaystyle \frac{d𝒔^n}{dr}}=\left({\displaystyle \frac{d\mathrm{ln}\rho _n}{d\mathrm{ln}r}}+3\right)𝒔^n+\mathrm{\Lambda }_0𝒉^n,`$ (63)
$`r{\displaystyle \frac{d𝒔^p}{dr}}=\left({\displaystyle \frac{d\mathrm{ln}\rho _n}{d\mathrm{ln}r}}+3\right)𝒔^p+\mathrm{\Lambda }_0𝒉^p.`$ (64)
The independent components of equations (26) and (27) are reduced to

$`L_0r{\displaystyle \frac{d𝒉^n}{dr}}M_1r{\displaystyle \frac{d𝒕^n}{dr}}m\nu \mathrm{\Lambda }_0^1r{\displaystyle \frac{d𝒔^n}{dr}}+\{(2+m\nu )\mathrm{𝟏}2m\nu \mathrm{\Lambda }_0^1\}𝒉^n(\mathrm{𝟏}+2m\nu \mathrm{\Lambda }_0^1)𝒔^p`$
$`+\nu \mathrm{\Lambda }_0^1C(\mathrm{\Lambda }_12\mathrm{𝟏})𝒕^n+{\displaystyle \frac{\rho _{12}}{\rho _n}}\left\{r{\displaystyle \frac{d}{dr}}(𝒉^p𝒉^n)+2(𝒉^p𝒉^n)(𝒔^p𝒔^n)\right\}`$ (65)
$`+r{\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{\rho _{12}}{\rho _n}}\right)(𝒉^p𝒉^n)=0,`$
$`L_0r{\displaystyle \frac{d𝒉^p}{dr}}M_1r{\displaystyle \frac{d𝒕^p}{dr}}m\nu \mathrm{\Lambda }_0^1r{\displaystyle \frac{d𝒔^p}{dr}}+\{(2+m\nu )\mathrm{𝟏}2m\nu \mathrm{\Lambda }_0^1\}𝒉^p(\mathrm{𝟏}+2m\nu \mathrm{\Lambda }_0^1)𝒔^p`$
$`+\nu \mathrm{\Lambda }_0^1C(\mathrm{\Lambda }_12\mathrm{𝟏})𝒕^p+{\displaystyle \frac{\rho _{12}}{\rho _p}}\left(1+{\displaystyle \frac{m_e}{m_p}}\right)^1\left\{r{\displaystyle \frac{d}{dr}}(𝒉^n𝒉^p)+2(𝒉^n𝒉^p)(𝒔^n𝒔^p)\right\}`$ (66)
$`+r{\displaystyle \frac{d}{dr}}\left\{{\displaystyle \frac{\rho _{12}}{\rho _p}}\left(1+{\displaystyle \frac{m_e}{m_p}}\right)^1\right\}(𝒉^n𝒉^p)=0,`$
$$L_1𝒕^nM_0𝒉^n+\nu K𝒔^n+\frac{\rho _{12}}{\rho _n}(𝒕^p𝒕^n)=0,$$
(67)
$`L_1𝒕^pM_0𝒉^p+\nu K𝒔^p+{\displaystyle \frac{\rho _{12}}{\rho _p}}\left(1+{\displaystyle \frac{m_e}{m_p}}\right)^1(𝒕^n𝒕^p)=0,`$ (68)
where $`\nu =2\mathrm{\Omega }/\sigma =2/\kappa _0`$ and $`\mathrm{𝟏}`$ stands for a unit matrix. Here, the quantities, $`L_0`$, $`L_1`$, $`C`$, $`K`$, $`\mathrm{\Lambda }_0`$, $`\mathrm{\Lambda }_1`$, $`M_0`$, and $`M_1`$ are matrices, given as follows:

For even modes,

$$(L_0)_{i,i}=1\frac{m\nu }{l(l+1)},(L_1)_{i,i}=1\frac{m\nu }{(l+1)(l+2)},$$

$$(C)_{i,i}=(l+2)J_{l+1}^m,(C)_{i+1,i}=(l+1)J_{l+2}^m,$$

$$(K)_{i,i}=\frac{J_{l+1}^m}{l+1},(K)_{i,i+1}=\frac{J_{l+2}^m}{l+2},$$

$$(\mathrm{\Lambda }_0)_{i,i}=l(l+1),(\mathrm{\Lambda }_1)_{i,i}=(l+1)(l+2),$$

$$(M_0)_{i,i}=\nu \frac{l}{l+1}J_{l+1}^m,(M_0)_{i,i+1}=\nu \frac{l+3}{l+2}J_{l+2}^m,$$

$$(M_1)_{i,i}=\nu \frac{l+2}{l+1}J_{l+1}^m,(M_1)_{i+1,i}=\nu \frac{l+1}{l+2}J_{l+2}^m,$$

where $`l=|m|+2i2`$ for $`i=1,2,3,\mathrm{}`$, and

$$J_l^m\left[\frac{(l+m)(lm)}{(2l1)(2l+1)}\right]^{1/2}.$$
(69)
For odd modes,

$$(L_0)_{i,i}=1\frac{m\nu }{l(l+1)},(L_1)_{i,i}=1\frac{m\nu }{l(l1)},$$

$$(C)_{i,i}=(l1)J_l^m,(C)_{i,i+1}=(l+2)J_{l+1}^m,$$

$$(K)_{i,i}=\frac{J_l^m}{l},(K)_{i+1,i}=\frac{J_{l+1}^m}{l+1},$$

$$(\mathrm{\Lambda }_0)_{i,i}=l(l+1),(\mathrm{\Lambda }_1)_{i,i}=l(l1),$$

$$(M_0)_{i,i}=\nu \frac{l+1}{l}J_l^m,(M_0)_{i+1,i}=\nu \frac{l}{l+1}J_{l+1}^m,$$

$$(M_1)_{i,i}=\nu \frac{l1}{l}J_l^m,(M_1)_{i,i+1}=\nu \frac{l+2}{l+1}J_{l+1}^m,$$

where $`l=|m|+2i1`$ for $`i=1,2,3,\mathrm{}`$.

The oscillation equations in the superfluid core which we use are given as a set of linear ordinary differential equations for the variables $`𝒔^X`$ and $`𝒕^X`$ with $`X=n,p`$, which are obtained by eliminating the vectors $`𝒉^X`$ in equations (A2)–(A5) using the algebraic equations (A6) and (A7). The oscillation equations in the normal fluid envelope can be obtained by substituting $`𝒔^n=𝒔^p=𝒔^o`$, $`𝒕^n=𝒕^p=𝒕^o`$ and $`\rho _n=\rho _p=\rho `$ into the oscillation equations in the superfluid core, where $`𝒔^o`$ and $`𝒕^o`$ means the fluid velocity fields in the normal fluid envelope.