text:
# On the dynamics of superfluid neutron star cores

## 1 Introduction

In the forty or so years since Migdal’s \[Migdal 1959\] initial work, three important areas of research on superfluidity in neutron stars have emerged: (i) Nuclear physics studies of the pairing gap energies as functions of mass density (for recent reviews, see \[Lombardo 1999, Lombardo & Schulze 2000\]); (ii) observations and subsequent modeling of glitches (see, for instance, \[Sauls 1989, Lyne 1993\] and references therein) and neutron star cooling \[Tsuruta 1994, Heiselberg & Hjorth-Jensen 2000\]; and (iii) studies of the dynamics of superfluid neutron star cores in the Newtonian \[Epstein 1988, Mendell & Lindblom 1991, Mendell 1991, Lindblom & Mendell 1994, Lindblom & Mendell 1995, Lee 1995, Prix 1999, Lindblom & Mendell 2000\] and general relativistic regimes \[Carter 1989, Comer & Langlois 1994, Carter & Langlois 1995a, Carter & Langlois 1995b, Carter & Langlois 1998, Langlois et al 1998, Comer et al 1999, Andersson & Comer 2001\]. It is this third area that is the focus of the present discussion. We want to understand how the pulsation properties of a neutron star changes when it cools below the temperature (a few times $`10^9\mathrm{K}`$) at which the bulk of the interior is expected to become superfluid. One of the underlying motivations for this work is the possibility that future gravitational-wave detectors may be sensitive enough to observe pulsating neutron stars, following (say) a glitch or a starquake \[Andersson & Kokkotas 1996, Andersson & Kokkotas 1998, Kokkotas et al 2001\]. An analysis of such observed data will provide a probe of the neutron star interior and should, at least in principle, allow us to infer the details of the supranuclear equation of state. In particular, one would hope to be able to establish beyond any doubt that the core of a neutron star contains a superfluid and perhaps constrain some of the relevant parameters for (say) entrainment. Obviously we can only hope to achieve this goal if the pulsation properties of a superfluid star differ significantly from those of an “ordinary fluid” neutron star.

The equations that describe ordinary (cold, inviscid etcetera) fluid neutron stars consist of a single mass continuity equation and a single Euler equation that determines the fluid velocity. But even this “simple” system yields an amazing diversity of modes of pulsation (see McDermott et al \[McDermott et al 1988\] for a fine discussion of many classes of modes), which includes the spheroidal f-, p-, and g-modes and the toroidal r-modes. In reality cold neutron stars are significantly more complicated than what is implied by the ordinary fluid model. The current understanding is that the outer regions of a mature neutron star consist of crust nuclei and an electron gas, which (beyond neutron drip) are everywhere permeated by superfluid neutrons. In the core, nuclei have dissolved leaving behind neutrons in a superfluid state, superconducting protons, and an ultra-relativistic gas of degenerate electrons. In addition, more exotic particles (e.g. hyperons) may be present. Deep in the core various condensates (of kaons, pions etcetera) may form, and deconfined quarks may also play a crucial role. A priori, one might expect the pulsation spectrum of such a “real” neutron star to be very complex, consisting of the various modes that exist for the ordinary fluid case, plus additional modes that arise because of new fluid degrees of freedom due to different particle species moving (more or less) independently of each other.

The most striking ways that (pure) superfluids differ from ordinary fluids is that superfluids have zero viscosity and are locally irrotational \[Tilley & Tilley 1986\]. The latter property, however, is compensated by the superfluid being threaded by an array of quantized vortices so that it can on macroscopic scales mimic the rotational behaviour of an ordinary fluid. The neutron superfluid in a rotating neutron star is believed to contain such threading. The vortices also lead to an effective viscosity through the so-called entrainment effect and a consequence of it known as mutual friction. In a mixture of the two superfluids Helium three and Helium four, it is known that a momentum induced in one of the constituents will cause some of the mass of the other to be carried along, or entrained \[Andreev & Bashkin 1975, Tilley & Tilley 1986\]. The analog in neutron stars is an entrainment of some of the protons, say, by the neutrons \[Alpar et al 1984, Borumand et al 1996\]. Because of entrainment, the flow of neutrons around the neutron fluid vortices will also induce a flow in a fraction of the protons, leading to magnetic fields being formed around the vortices. But since the electrons are coupled to the protons on very short timescales \[Alpar et al 1984\], some will track closely the entrained protons. Mutual friction is the dissipative scattering of these electrons off of the magnetic fields associated with the vortices.

The comparatively simple model of neutron star superfluidity that will be used in this paper considers just two fluids, and the entrainment effect that acts between them. One fluid consists of the superfluid neutrons that exist in the inner crust and core, whereas the other is a conglomerate of all the charged particles (i.e. the nuclei and electrons in the crust, and the superconducting protons and electrons in the core) that will be loosely referred to as “protons.” In comparison to the ordinary fluid case, the superfluid system of equations will consist of two mass continuity equations, and two Euler equations that determine the neutron and proton fluid velocities. This model has been used in previous numerical studies of the linearized pulsations of superfluid neutron stars \[Lindblom & Mendell 1994, Lee 1995, Comer et al 1999\]. However, it is important to point out that Mendell \[Mendell 1991\] has discussed the more general case where electromagnetic effects are explicitly included, and where the electrons (and even muons) are free to move independently of the protons, a net result of which is to delimit the dynamical timescales for which the simplified two-fluid model is appropriate.

It is also relevant to mention that even before the equations were analyzed in detail, Epstein \[Epstein 1988\] (and then later Mendell \[Mendell 1991\]) argued, using a simple counting of the fluid degrees of freedom in a superfluid neutron star, that there should be a new class of modes that were later dubbed superfluid modes. They arise primarily because the neutrons and protons in the core are superfluid and no longer locked together in nuclei, thus leading to an increase in the number of fluid degrees of freedom. It is also impressive that—using as an analogy a system of coupled pendulums—Mendell \[Mendell 1991\] argued that the distinguishing characteristic of the new modes should be the counter-motion of the neutrons with respect to the protons. That is, in terms of a projection along the radial direction, the superfluid modes should have the protons moving oppositely to the neutrons, unlike an ordinary fluid mode that would have the neutrons and protons moving more or less in “lock-step.” This picture has been confirmed by numerical (and simplified analytical) studies in both the Newtonian \[Lindblom & Mendell 1994, Lee 1995\] and general relativistic regimes \[Comer et al 1999\]. Another feature, shown by Lindblom and Mendell \[Lindblom & Mendell 1994\], is that the superfluid modes are driven by deviations from chemical equilibrium between the neutrons and protons.

At first sight, the existence of the superfluid modes appears to confirm one’s intuition that the new fluid degrees of freedom would double the number of pulsation mode-families. However, there is an open question in the literature regarding the g-modes in the superfluid case. In particular, Lee’s \[Lee 1995\] numerical results did not reveal any independent set of pulsating g-modes, much less a doubling. This is perhaps strange, and certainly not consistent with the simple counting argument given above, since Reisenegger and Goldreich \[Reisenegger & Goldreich 1992\] have shown convincingly that a composition gradient (for instance, a stable stratification in terms of the proton fraction) leads to the existence of propagating g-modes in non-rotating cold ordinary fluid neutron stars. Even more puzzling is the fact that Lee showed that the ordinary fluid definition of the so-called Brunt-Väisälä frequency leads to real values, which is the usual indication that pulsating g-modes should exist. His final conclusion was that the superfluid modes are to be found among the ordinary fluid f- and p-modes, and that there are no propagating g-modes in the cores of superfluid neutron stars, even if there exists a composition gradient. Here, we will use a “local” analysis of the linearized superfluid equations to show that superfluid modes are predominately acoustic and to confirm Lee’s numerical results that superfluidity prevents g-modes from pulsating.

Of course, if there are no distinct g-modes then one is still left with the question of where they have gone. To find the answer, we will analyze the zero-frequency subspace of solutions to the superfluid equations of motion, following recent studies by Lockitch and Friedman \[Lockitch & Friedman 1999\] and Lockitch et al \[Lockitch et al 2001\] in the ordinary fluid case. This subspace is composed of the time-independent solutions to the linearized equations. We find that the zero-frequency subspace appropriate to the superfluid equations is spanned by solutions that separate into two distinct classes: (i) Non-trivial local scalar matter (e.g. number densities, pressure, etc) perturbations accompanied by zero velocity perturbations, and (ii) degenerate velocity perturbations with vanishing local matter perturbations. The first class of solutions simply takes static and spherically symmetric configurations and deforms them into other, nearby static and spherically symmetric configurations. The second class separates into independent sets of spheroidal and toroidal velocity perturbations. This is exacly analogous to the results for an ordinary fluid, but there is one important difference. We find two sets of decoupled perturbations for each of the spheroidal/toroidal velocity perturbations.

We interpret these results in the following way: Because the second class of solutions have vanishing local matter variations, it is natural to take the spheroidal solutions to be the missing g-modes. The existence of two independent sets of such solutions would then be in agreement with the intuitive notion that the fluid degrees of freedom are doubled in the superfluid case. The presence of spheroidal modes in the zero-frequency subspace is analogous to the ordinary fluid case when there is no stratification in the background star \[Lockitch & Friedman 1999, Lockitch et al 2001\]. The doubling of modes in the class of toroidal solutions is also unique to the superfluid. Given these results, an analysis of rotating configurations becomes of prime importance. A comparison to the ordinary fluid case \[Lockitch & Friedman 1999, Lockitch et al 2001\]) suggests that rotation will break the degeneracy of the various zero-frequency and lead to inertial modes with true dynamics. Our results for the zero-frequency subspace suggests that, just like in a barotropic ordinary fluid star, the typical inertial mode will be a hybrid mixture of spheroidal and toroidal velocity components in the non-rotating limit.

We take the first steps towards the study of inertial modes in a superfluid star by considering the so-called r-modes (which correspond to purely toroidal velocities in the non-rotating limit). These modes are of particular relevance since it was recently discovered \[Andersson 1998, Friedman & Morsink 1998\] that they are generically unstable to gravitational radiation via the so-called CFS (Chandrasekhar-Friedman-Schutz) mechanism \[Chandrasekhar 1970, Friedman & Schutz 1978, Friedman 1978\]. In principle, the CFS mechanism can apply to any mode in a rotating neutron star, but viscous damping of various kinds tends to kill the instability before substantial amounts of gravitational radiation are emitted, except, apparently, for the r-modes \[Lindblom et al 1998, Andersson et al 1999, Owen et al 1998, Andersson & Kokkotas 2001\]. Remarkably, even though mutual friction in superfluid neutron stars is extremely effective at killing the CFS instability for f- and p-modes, Lindblom and Mendell \[Lindblom & Mendell 2000\] have shown that it does not damp out the r-modes (except for a very small subset of the entrainment models employed). But the main point for the present discussion is whether or not superfluidity leads to a new set of r-modes, thus making the problem richer than the normal fluid one. Unfortunately, this is hard to discern from the discussion of Lindblom and Mendell \[Lindblom & Mendell 2000\]. While they argue for the existence of some type of superfluid r-modes from their numerical results, in the sense that they see the counter-motion of the neutrons with respect to the protons, they find that these modes differ little from the ordinary fluid r-modes at both the lowest- and second-order in the angular velocity of the background configuration. Also, they argue that any analogy with the spheroidal type of superfluid modes discussed earlier is only superficial, since the fluid motion is not dominated by a deviation from chemical equilibrium. We will argue below that a distinct set of superfluid r-modes does, in fact, exist. Moreover, we will show that the frequencies of these modes, even at the lowest-order in angular velocity, differ from those of the ordinary fluid r-modes, and depend in an essential way on the strength of the entrainment effect. These results shed light on some of the (so far unexplained) peculiarities on the results regarding the effect of mutual friction on the unstable r-modes.

Carter, Langlois and collaborators \[Carter 1989, Carter & Langlois 1995a, Carter & Langlois 1995b, Carter & Langlois 1998, Langlois et al 1998, Comer & Langlois 1994\] have been developing a system of equations that can describe superfluids in general relativistic neutron stars. They are the relativistic analogs of the equations used by Lindblom and Mendell and Lee for the Newtonian regime. Until recently, the relativistic equations had not been used for modeling various scenarios but this situation is rapidly changing. Comer, Langlois and Lin \[Comer et al 1999\] have used the relativistic equations to model the quasinormal modes of a simple neutron star model, and were successful at extracting the general relativistic version of the superfluid modes. Andersson and Comer \[Andersson & Comer 2001\] have built models of “slowly rotating” general relativistic superfluid neutron stars. A unique feature of their formalism is that the neutrons are not a priori forced to corotate with the protons. We will use the Newtonian limit of these equations for the main analysis presented below. One reason for this is to show how to connect the general relativistic formalism with the Newtonian one. Another underlying motivation is to show how to incorporate the current microscopic models of entrainment in the general relativistic case.

The layout of the paper is as follows: In Section II, we review the standard reasoning that leads to a demonstration of p- and g-modes in cold (single fluid) neutron stars. This will provide a convenient background for the analysis of the following three sections (III, IV, and V) of the perturbations of superfluid neutron stars. Aside from the main points discussed above, a practical outcome of these calculations will be a simple, approximate formula for the superfluid mode frequencies that has not been given previously in the literature. Largely just to simplify the analysis we use the so-called Cowling approximation, i.e. the variation in the gravitational field will be ignored. In Section VI we unambiguously demonstrate the existence of two distinct classes of r-modes in a superfluid star. A brief summary is given in Section VII. Finally, in the Appendix we give the steps used to derive the Newtonian limit of the general relativistic superfluid field equations, and comment on the relation between our formulae and those used by Lindblom and Mendell.

## 2 Stellar pulsation primer

### 2.1 The standard picture: p/g-modes

Before turning to the superfluid case, it is useful to recall the “standard” theory of stellar pulsation \[Cox 1980, Unno et al 1989\]. For an ordinary fluid star, we need to consider the linearized version of the perturbed Euler equation

$$\frac{^2\stackrel{}{\xi }}{t^2}=\frac{\delta \rho }{\rho ^2}P\frac{1}{\rho }\delta P\delta \mathrm{\Phi }$$
(1)
(where $`\stackrel{}{\xi }`$ is the fluid displacement vector, $`P`$ denotes the pressure, $`\rho `$ is the density and $`\mathrm{\Phi }`$ is the gravitational potential) as well as the (integrated form of) the continuity equation

$$\delta \rho +(\rho \stackrel{}{\xi })=0.$$
(2)
We are using $`\delta `$ to denote Eulerian perturbations, while $`\mathrm{\Delta }`$ will be used (in this section only) to indicate Lagrangian variations.

It is customary to introduce the adiabatic index of the perturbations $`\mathrm{\Gamma }_1`$ via

$$\frac{\mathrm{\Delta }P}{P}=\mathrm{\Gamma }_1\frac{\mathrm{\Delta }\rho }{\rho }$$
(3)
or, in terms of the Eulerian variations,

$`\delta P`$ $`=`$ $`{\displaystyle \frac{P\mathrm{\Gamma }_1}{\rho }}\delta \rho +\stackrel{}{\xi }\left[{\displaystyle \frac{P\mathrm{\Gamma }_1}{\rho }}\rho P\right]`$ (4)
$``$ $`{\displaystyle \frac{P\mathrm{\Gamma }_1}{\rho }}\delta \rho +P\mathrm{\Gamma }_1(\stackrel{}{\xi }\stackrel{}{A})`$ (6)
which defines the Schwarzschild discriminant $`\stackrel{}{A}`$.

For spherical stars we can now readily rewrite the Euler equation as

$$\frac{^2\stackrel{}{\xi }}{t^2}=\left(\frac{\delta P}{\rho }\right)+\frac{P\mathrm{\Gamma }_1}{\rho }\stackrel{}{A}(\stackrel{}{\xi })\delta \mathrm{\Phi }.$$
(7)
Once the equation is written in this form we can deduce that the fluid motion is affected by (neglecting $`\delta \mathrm{\Phi }`$) two restoring forces: the pressure variation and the buoyancy associated with $`\stackrel{}{A}`$. The latter is relevant whenever the star is stratified, either by entropy or compositional variations. Since we are considering neutron stars we can to good approximation assume that the temperature is zero. This means that we can neglect any internal entropy gradients. Still, as was pointed out by Reisenegger and Goldreich \[Reisenegger & Goldreich 1992\], we cannot assume that $`\stackrel{}{A}=0`$ since any variation of the internal composition will lead to an effective buoyancy force acting on a fluid element. For neutron stars, the contribution due to the varying proton fraction is likely to be the most important effect.

We want to study (7) and try to infer the nature of the various modes of pulsation that the star may have. In doing this we assume that the Cowling approximation is used, i.e. we neglect the variation $`\delta \mathrm{\Phi }`$ in the gravitational potential. If we further assume that the fluid element remains in pressure equilibrium with its surroundings (in such a way that $`\delta P=0`$) we have

$$\mathrm{\Delta }P\stackrel{}{\xi }P\rho \stackrel{}{\xi }\stackrel{}{g},$$
(8)
where we have used the standard definition of the local gravitational acceleration $`\stackrel{}{g}`$. If we also use the continuity equation,

$$\mathrm{\Delta }\rho \delta \rho +\stackrel{}{\xi }\rho =\rho \stackrel{}{\xi },$$
(9)
we find that the radial component of the Euler equation becomes

$$\frac{^2\xi ^r}{t^2}=g\xi ^rAN^2\xi ^r$$
(10)
which defines the so-called Brunt-Väisälä frequency $`N`$. In other words, we have oscillatory motion whenever $`N^2>0`$. The resultant modes of pulsation are known as the g-modes, as they are essentially governed by gravity and internal stratification in the star. Whenever $`N^2<0`$ the perturbation is unstable.

We can study the nature of the star’s modes in more detail in the following way (cf. Reisenegger & Goldreich \[Reisenegger & Goldreich 1992\]). Since $`\stackrel{}{g}`$ is purely radial for a spherical star the horizontal component of the Euler equation (7) leads to

$$\stackrel{}{\xi }_{}=\frac{1}{\omega ^2\rho }_{}\delta P,$$
(11)
where we have assumed that the perturbation has time-dependence $`e^{i\omega t}`$. Combining this with (9) we find that

$`\mathrm{\Delta }\rho `$ $`=`$ $`{\displaystyle \frac{\rho }{r^2}}{\displaystyle \frac{}{r}}(r^2\xi ^r)\rho _{}\stackrel{}{\xi }`$ (12)
$`=`$ $`{\displaystyle \frac{\rho }{r^2}}{\displaystyle \frac{}{r}}(r^2\xi ^r)+{\displaystyle \frac{l(l+1)}{\omega ^2r^2}}\delta P,`$ (14)
where we have assumed that the angular dependence of $`\delta P`$ can be represented by a single spherical harmonic $`Y_{lm}(\theta ,\phi )`$ (which is natural since the pressure variation is a scalar) such that

$$_{}_{}\delta P=\frac{l(l+1)}{r^2}\delta P.$$

We also have the radial component of (7):

$$\omega ^2\xi ^r=\frac{}{r}\left(\frac{\delta P}{\rho }\right)\frac{P\mathrm{\Gamma }_1A\mathrm{\Delta }\rho }{\rho ^2}.$$
(15)
In these two equations we can replace $`\mathrm{\Delta }\rho `$ by noticing that

$$\mathrm{\Delta }\rho =\frac{\rho }{P\mathrm{\Gamma }_1}(\delta P\rho g\xi ^r)$$
(16)
(having used the fact that $`\stackrel{}{g}=g\widehat{e}_r`$), and we get

$$\frac{1}{r^2}\frac{}{r}(r^2\xi ^r)\frac{\rho g}{P\mathrm{\Gamma }_1}\xi ^r=\left[\frac{l(l+1)}{\omega ^2r^2}\frac{\rho }{P\mathrm{\Gamma }_1}\right]\frac{\delta P}{\rho }$$
(17)
and

$$\frac{1}{\rho }\frac{}{r}\delta P+\frac{g}{P\mathrm{\Gamma }_1}\delta P=(\omega ^2+gA)\xi ^r.$$
(18)
Let us now introduce new variables

$`\widehat{\xi }^r`$ $`=`$ $`{\displaystyle \frac{r^2\xi ^r}{\varphi }},`$ (19)
$`\delta \widehat{P}`$ $`=`$ $`\varphi \delta P,`$ (20)
where

$$\varphi =\mathrm{exp}\left[g/c_s^2𝑑r\right]$$
(21)
and the sound speed is

$$c_s^2=\frac{\mathrm{\Delta }P}{\mathrm{\Delta }\rho }=\frac{P\mathrm{\Gamma }_1}{\rho }.$$
(22)
With these definitions our two equations can be written

$$\frac{\widehat{\xi }^r}{r}=\left[L_l^2\omega ^2\right]\frac{r^2\delta \widehat{P}}{\rho \omega ^2c_s^2\varphi ^2},$$
(23)
where we have introduced the so-called Lamb frequency

$$L_l^2=\frac{l(l+1)c_s^2}{r^2},$$
(24)
and

$$\frac{\delta \widehat{P}}{r}=\left[\omega ^2+gA\right]\frac{\rho \widehat{\xi }^r\varphi ^2}{r^2}.$$
(25)
Given this we can reduce the problem to the following ordinary differential equation for $`\widehat{\xi }^r`$:

$`{\displaystyle \frac{d}{dr}}\left\{{\displaystyle \frac{\rho \omega ^2c_s^2\varphi ^2}{r^2}}\left[L_l^2\omega ^2\right]^1{\displaystyle \frac{d\widehat{\xi }^r}{dr}}\right\}`$ (26)
(27)
$`[\omega ^2N^2]{\displaystyle \frac{\rho \widehat{\xi }^r\varphi ^2}{r^2}}=0.`$ (28)
From this we can draw some very important conclusions. We see that the problem reduces to the Sturm-Liouville form both for high and low frequencies. For large $`\omega ^2`$ we get

$$\frac{d}{dr}\left\{\frac{\rho c_s^2\varphi ^2}{r^2}\frac{d\widehat{\xi }^r}{dr}\right\}+[\omega ^2N^2]\frac{\rho \widehat{\xi }^r\varphi ^2}{r^2}=0.$$
(29)
Then standard Sturm-Liouville theory tells us that there will be an infinite set of modes which can be labelled by the number of radial nodes ($`n`$) of the various eigenfunctions, and for which $`\omega _n\mathrm{}`$ as $`n\mathrm{}`$. In the opposite limit, when $`\omega ^2`$ is small, the problem becomes

$$\frac{d}{dr}\left\{\frac{\rho \varphi ^2}{l(l+1)}\frac{d\widehat{\xi }^r}{dr}\right\}+\left[\frac{N^2}{\omega ^2}1\right]\frac{\rho \widehat{\xi }_r\varphi ^2}{r^2}=0$$
(30)
and we deduce that there will be another set of modes, with eigenfrequencies such that $`\omega _n0`$ as $`n\mathrm{}`$.

These sets of modes are the p- and g-modes, respectively, and we can estimate their frequencies in the following way: Assume that the perturbations have a characteristic wavelength $`k^1`$ (such that the various functions are proportional to $`\mathrm{exp}(ikr)`$). Then we can readily deduce the simple dispersion relation

$$k^2=\frac{1}{c_s^2\omega ^2}(N^2\omega ^2)(L_l^2\omega ^2).$$
(31)
Here we must have $`\omega ^2>0`$ for stability, and we see that we will have pulsating modes ($`k^2>0`$) in two different cases. Either $`\omega ^2`$ must be smaller than both $`N^2`$ and $`L_l^2`$, or it must be larger than both these quantities. For $`l>>kr`$ we estimate the mode-frequencies as $`\omega ^2L_l^2`$ for the p-modes and $`\omega ^2N^2`$ for the g-modes.

So far, we have discussed familiar textbook results \[Cox 1980, Unno et al 1989\]. We did this in order to be able to compare and constrast these results with the corresponding ones in the superfluid case. This is the aim of the rest of the paper.

### 2.2 The first step towards two fluids

As a first step towards considering a superfluid star we approach the above pulsation problem in a somewhat indirect way. Let us suppose that the star is composed of two distinct species of particles, which we will think of as the superfluid neutrons in the core and the “protons” as defined in the Introduction. We denote their respective number densities by $`n_n`$ and $`n_p`$. Then standard thermodynamical considerations tell us that the pressure can be determined from

$$dP=n_nd\mu _n+n_pd\mu _p,$$
(32)
where $`\mu _i`$ ($`i=n,p`$) are the two chemical potentials (which do not include the rest-masses). In general, there should also be a term proportional to $`d|\stackrel{}{v}_n\stackrel{}{v}_p|^2`$, where $`\stackrel{}{v}_n`$ and $`\stackrel{}{v}_p`$ are the neutron and proton velocities (cf. the Appendix), respectively, but since the background velocities are zero (or at least the same for both fluids) in the cases we will consider, such a term will not contribute in what follows. By introducing

$$\stackrel{~}{\mu }_n=\frac{\mu _n}{m_n},\stackrel{~}{\mu }_p=\frac{\mu _p}{m_p},$$
(33)
we get the following relation between the Eulerian pressure variation and the variations in the two chemical potentials:

$$\delta P=\rho _n\delta \stackrel{~}{\mu }_n+\rho _p\delta \stackrel{~}{\mu }_p,$$
(34)
where $`\rho _n=m_nn_n`$ and $`\rho _p=m_pn_p`$.

Chemical equilibrium for our system corresponds to $`\mu _n=\mu _p`$. It is natural to introduce a variable that describes the deviation from equilibrium introduced by the fluid motion. Thus we define $`\delta \beta `$ by

$$\delta \stackrel{~}{\mu }_p=\delta \stackrel{~}{\mu }_n+\delta \beta .$$
(35)
For simplicity, we have assumed that the two particle masses are equal: $`m_n=m_p=m`$. To define $`\delta \beta `$ in this particular way (in terms of the tilde variables) may seem peculiar, but it is useful since this variable then has exactly the same meaning as $`\delta \beta `$ in the series of papers by Lindblom and Mendell \[Lindblom & Mendell 1994, Lindblom & Mendell 1995\]. Anyway, the above relations enable us to write

$$\delta \stackrel{~}{\mu }_n=\frac{\delta P}{\rho }\frac{\rho _p}{\rho }\delta \beta .$$
(36)
The corresponding thermodynamic condition is that

$$d\stackrel{~}{\mu }_n=\frac{1}{\rho }dP\frac{\rho _p}{\rho }d\beta .$$
(37)
From this we can immediately read off that

$$\frac{1}{\rho }=\left(\frac{\stackrel{~}{\mu }_n}{P}\right)_\beta ,\frac{\rho _p}{\rho }=\left(\frac{\stackrel{~}{\mu }_n}{\beta }\right)_P,$$
(38)
and using the equality of mixed partial derivatives, we can deduce the useful identity:

$$\rho ^2\frac{}{P}\left(\frac{\rho _p}{\rho }\right)_\beta =\left(\frac{\rho }{\beta }\right)_P.$$
(39)
Let us now assume that the motion of each species of particles is determined by the variations in the chemical and gravitational potentials according to (see the Appendix for a justification)

$$\frac{\delta \stackrel{}{v}_n}{t}+(\delta \stackrel{~}{\mu }_n+\delta \mathrm{\Phi })=0$$
(40)
and

$$\frac{\delta \stackrel{}{v}_p}{t}+(\delta \stackrel{~}{\mu }_p+\delta \mathrm{\Phi })=0.$$
(41)
In view of the earlier relations, the second of the two equations can be written

$$\frac{\delta \stackrel{}{v}_p}{t}+(\delta \stackrel{~}{\mu }_n+\delta \beta +\delta \mathrm{\Phi })=0.$$
(42)
Given (40) and (42) we can make an important observation: The neutrons and protons will only move together (in the sense that $`\delta \stackrel{}{v}_n=\delta \stackrel{}{v}_p=_t\stackrel{}{\xi }`$) if $`\delta \beta =0`$. Intuitively this is obvious, since it simply says that the perturbation keeps the fluid in chemical equilibrium if the two species move together. Yet, it will prove a useful observation later. We note that if the two particle species move together we retain the standard Euler equation (7) in the particular case $`\stackrel{}{A}=0`$. This is the first hint of a result that has repercussions on the pulsation properties of superfluid stars. As we will now show, the case of barotropic perturbations ($`\stackrel{}{A}=0`$) in the standard description, corresponds to $`\delta \beta =0`$ in the two-fluid picture.

We elucidate this correspondence in the following way: Introducing a new variable corresponding to the average velocity

$$\frac{\stackrel{}{\xi }_+}{t}=\frac{\rho _n}{\rho }\delta \stackrel{}{v}_n+\frac{\rho _p}{\rho }\delta \stackrel{}{v}_p,$$
(43)
(i.e. $`\rho _t\stackrel{}{\xi }_+`$ is the total mass-density current) and combining (40) and (42), we readily get

$$\frac{^2\stackrel{}{\xi }_+}{t^2}+\delta \stackrel{~}{\mu }_n+\delta \mathrm{\Phi }+\frac{\rho _p}{\rho }\delta \beta =0$$
(44)
or

$$\frac{^2\stackrel{}{\xi }_+}{t^2}+\left(\frac{\delta P}{\rho }\right)+\delta \mathrm{\Phi }\frac{}{P}\left(\frac{\rho _p}{\rho }\right)\delta \beta P=0.$$
(45)
Comparing (45) to the standard Euler equation (7), and identifying $`\stackrel{}{\xi }_+=\stackrel{}{\xi }`$, we see that we must have

$$\frac{}{P}\left(\frac{\rho _p}{\rho }\right)_\beta \delta \beta P=\frac{P\mathrm{\Gamma }_1}{\rho }\stackrel{}{A}(\stackrel{}{\xi }_+).$$
(46)
In other words, the one-fluid Schwarzschild discriminant is intimately linked to the variation $`\delta \beta `$ in the two-fluid picture, and indeed $`\delta \beta =0`$ implies $`\stackrel{}{A}=0`$. This is quite natural since (in the case of a chemical composition gradient) the Schwarzschild discriminant describes the extent to which a given perturbation drives a fluid element away from chemical equilibrium. This observation implies that the pulsation properties of a two-fluid model with $`\delta \beta =0`$ should be analogous to the case of barotropic perturbations for which there are no non-trivial g-modes. In other words, in this case one expects only to find a set of p-modes.

## 3 The equations governing a perturbed superfluid

We now consider the full two-fluid model for superfluid neutron stars, i.e. consider the superfluid neutrons as being coupled to a conglomerate of charged particles (protons, electrons, nuclei etcetera). Adopting the notation of the discussion in the previous section we first rewrite the linearized form of the equations (211)-(222) from the Appendix in terms of a set of variables that are intimately linked to the physics of the system. The first two of these are the average velocity $`_t\stackrel{}{\xi }_+`$ and the deviation from beta-equilibrium $`\delta \beta `$ that we introduced in Section 2. Together with these we use the pressure variation $`\delta P`$ and a variable reflecting the relative motion of the neutrons and the protons. In order to retain complete correspondence with the equations used by Lindblom and Mendell<sup>1</sup><sup>1</sup>1This is obviously not necessary but we want to avoid unnecessary confusion among readers who are well acquainted with the relevant literature. we define

$$\frac{\stackrel{}{\xi }_{}}{t}=\frac{\rho _n\rho _p}{det\rho }(\delta \stackrel{}{v}_p\delta \stackrel{}{v}_n),$$
(47)
where we note that the case $`\stackrel{}{\xi }_{}=0`$ corresponds to the two species of particles moving together. We have introduced

$`det\rho `$ $`=`$ $`\rho _{nn}\rho _{pp}\rho _{np}^2,`$ (48)
$`\rho _n`$ $`=`$ $`\rho _{nn}+\rho _{np},`$ (50)
$`\rho _p`$ $`=`$ $`\rho _{pp}+\rho _{np}.`$ (52)
The mass-density matrix element $`\rho _{np}`$ is what allows for the entrainment effect \[Andreev & Bashkin 1975\], and it can be rewritten as

$$\rho _{np}=ϵ\rho _n,$$
(53)
where $`ϵ`$ is the so-called entrainment parameter. We also have

$$\rho _{nn}=\rho _n(1+ϵ),\rho _{pp}=\rho _p+ϵ\rho _n.$$
(54)
In the model used by Lindblom and Mendell \[Lindblom & Mendell 2000\] one has

$$ϵ=\frac{\rho _p}{\rho _n}\left(\frac{m_p}{m_p^{}}1\right),$$
(55)
where $`m_p^{}`$ is the effective proton mass (which enters because the protons form a Fermi liquid and it is thus their associated quasiparticle features that are paramount, see Sjöberg ).

The first continuity equation (i.e. the sum of the linearized versions of (211) and (212) of the Appendix) takes the form

$$\delta \rho +(\rho \stackrel{}{\xi }_+)=0$$
(56)
and we also arrive (by adding the linearized versions of (217) and (222) of the Appendix) at the same equation of motion as in Section II:

$$\frac{^2\stackrel{}{\xi }_+}{t^2}+\left(\frac{\delta P}{\rho }\right)\frac{1}{\rho ^2}\left(\frac{\rho }{\beta }\right)_P\delta \beta P=0.$$
(57)
We have also assumed (for reasons of clarity) that the Cowling approximation is made, i.e. $`\delta \mathrm{\Phi }=0`$.

In addition (by subtracting the linearized versions of (217) and (222) of the Appendix) we have a second equation of motion relating $`\stackrel{}{\xi }_{}`$ and $`\delta \beta `$:

$$\frac{^2\stackrel{}{\xi }_{}}{t^2}+\delta \beta =0$$
(58)
and (after some manipulations<sup>2</sup><sup>2</sup>2The derivation of this equation is somewhat involved. The first step is to “invert” the velocities to find    $`\delta \stackrel{}{v}_n`$ $`=`$ $`{\displaystyle \frac{}{t}}\left(\stackrel{}{\xi }_+{\displaystyle \frac{det\rho }{\rho \rho _n}}\stackrel{}{\xi }_{}\right),`$    $`\delta \stackrel{}{v}_p`$ $`=`$ $`{\displaystyle \frac{}{t}}\left(\stackrel{}{\xi }_++{\displaystyle \frac{det\rho }{\rho \rho _p}}\stackrel{}{\xi }_{}\right).`$   Next one inserts the relationships obtained from (242) into the linearized form of the proton mass-continuity equation, using also the inverted velocities given above. Finally, one uses the linearized form of the neutron mass-continuity equation to obtain an equation for $`_i\xi _+^i`$.) the final continuity equation can be written

$`\rho ^2{\displaystyle \frac{}{P}}\left({\displaystyle \frac{\rho _p}{\rho }}\right)_\beta {\displaystyle \frac{\delta P}{\rho }}+{\displaystyle \frac{\rho _n^2}{\rho }}{\displaystyle \frac{}{\beta }}\left({\displaystyle \frac{\rho _p}{\rho _n}}\right)_P\delta \beta +`$ (60)
(61)
$`\rho \stackrel{}{\xi }_+\left({\displaystyle \frac{\rho _p}{\rho }}\right)+\left({\displaystyle \frac{det\rho }{\rho }}\stackrel{}{\xi }_{}\right)=0.`$ (62)
These four equations are identical to those derived and used by Lindblom and Mendell \[Lindblom & Mendell 1994, Lindblom & Mendell 1995\].

In these equations it is, however, difficult to discern the meaning, and thus importance, of the various thermodynamic derivatives that appear. Fortunately, the situation can be clarified somewhat by introducing the two local sound speeds that our two-fluid system of equations admit (just like their mathematical twins, the superfluid equations constructed for superfluid helium by Landau \[Putterman 1974\]). Using the standard technique of expanding the fluid variables as plane waves, and considering the background fluid to be homogeneous and at rest, we find that these sound speeds are obtained as solutions for $`u^2`$ from the following quadratic:

$`0=[u^2c_n^2][u^2{\displaystyle \frac{m_p}{m_p^{}}}c_p^2]+{\displaystyle \frac{\rho _p}{\rho _n}}[({\displaystyle \frac{m_p}{m_p^{}}}1)\times `$ (63)
(64)
$`([2c_{np}^2c_n^2]u^2+c_n^2c_p^2[1+{\displaystyle \frac{\rho _p}{\rho _n}}]c_{np}^4)c_{np}^4],`$ (65)
where

$`c_n^2`$ $``$ $`\rho _n{\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _n}},`$ (66)
$`c_p^2`$ $``$ $`\rho _p{\displaystyle \frac{\stackrel{~}{\mu }_p}{\rho _p}},`$ (68)
$`c_{np}^2`$ $``$ $`\rho _n{\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _p}}=\rho _n{\displaystyle \frac{\stackrel{~}{\mu }_p}{\rho _n}}.`$ (70)
When the ratio $`\rho _p/\rho _n`$ is small we see that the sound speeds are essentially $`c_n^2`$ and $`(m_p/m_p^{})c_p^2`$. In terms of these definitions we can now write the various thermodynamic derivatives in the (still exact) forms

$`{\displaystyle \frac{1}{c_{\mathrm{eq}}^2}}`$ $``$ $`\left({\displaystyle \frac{\rho }{P}}\right)_\beta ={\displaystyle \frac{x_n}{c_n^2}}[1{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_{np}}{c_n}}\right)^2\left({\displaystyle \frac{c_{np}}{c_p}}\right)^2]^1\times `$ (73)
$`\left[1+{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_n}{c_p}}\right)^22{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_{np}}{c_p}}\right)^2\right],`$
$`\left({\displaystyle \frac{\rho }{\beta }}\right)_P`$ $`=`$ $`\rho _n{\displaystyle \frac{x_p}{c_p^2}}[1{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_{np}}{c_n}}\right)^2\left({\displaystyle \frac{c_{np}}{c_p}}\right)^2]^1\times `$ (76)
$`\left[1\left({\displaystyle \frac{c_p}{c_n}}\right)^2\left(1{\displaystyle \frac{x_p}{x_n}}\right)\left({\displaystyle \frac{c_{np}}{c_n}}\right)^2\right],`$
$`{\displaystyle \frac{}{\beta }}\left({\displaystyle \frac{\rho _p}{\rho _n}}\right)_P`$ $`=`$ $`{\displaystyle \frac{x_p}{c_p^2}}[1{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_{np}}{c_n}}\right)^2\left({\displaystyle \frac{c_{np}}{c_p}}\right)^2]^1\times `$ (79)
$`\left[1+{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_p}{c_n}}\right)^2+2{\displaystyle \frac{x_p}{x_n}}\left({\displaystyle \frac{c_{np}}{c_n}}\right)^2\right],`$
and

$$\frac{1}{\mathrm{\Delta }}=\frac{c_n^2c_p^2}{\rho ^2x_nx_p}\left[1\frac{x_p}{x_n}\left(\frac{c_{np}}{c_n}\right)^2\left(\frac{c_{np}}{c_p}\right)^2\right],$$
(80)
where $`x_{n,p}=\rho _{n,p}/\rho `$ and we have used (233)–(240) given in the Appendix. (Also, one should not confuse the use of the $`\mathrm{\Delta }`$ symbol here with its earlier role as a Lagrangian variation.) The utility of these expressions is that they will illuminate the role of the proton fraction in determining the order of magnitude contributions from individual terms in the dispersion relations that will be written below.

An identity that follows from (73)–(80) above, and which will be used later, is

$$c_{\mathrm{eq}}^2\left(\frac{\rho }{\beta }\right)_P^2\rho _n^2\frac{}{\beta }\left(\frac{\rho _p}{\rho _n}\right)_P=c_{\mathrm{eq}}^2\mathrm{\Delta }.$$
(81)
In anticipation of later results, and in order to facilitate comparisons between the ordinary fluid and superfluid cases, we define for the superfluid system of equations two frequencies:

$$_l^2\frac{l(l+1)c_{\mathrm{eq}}^2}{r^2}$$
(82)
and

$$𝒩^2\frac{g^2}{c_{\mathrm{eq}}^2\mathrm{\Delta }}\left(\frac{\rho }{\beta }\right)_P^2.$$
(83)
The first of these corresponds to the standard Lamb frequency, while the second has (as will be discussed later) similar character to the Brunt-Väisälä frequency.

## 4 Oscillations of nonrotating superfluid stars

We now want to analyze the superfluid perturbation equations in a way similar to that used in Section II for the standard pulsation problem. There are several reasons for doing this. The most obvious one being that we want to contrast the pulsation properties of a superfluid star with the standard ordinary fluid results in order to see whether future observations may be able to distinguish between the two cases. From the general nature of the equations one might expect that the character of the various modes of oscillation may be rather different in the two cases. After all, in the superfluid case we are allowing the two fluids to move more or less independently and so would seem to have brought in additional fluid dynamical degrees of freedom. On the other hand, the g-modes in the standard case depend crucially on the stable stratification mainly associated with the varying chemical composition \[Reisenegger & Goldreich 1992\]. It is thus not at all clear what will happen to these modes if we allow the neutrons and the protons to move relative to one another.

### 4.1 Two simple limiting cases

Just as in the ordinary fluid case it is interesting, and potentially instructive, to consider what happens if we freeze various degrees of freedom. For example, the equations governing $`\stackrel{}{\xi }_+`$ and $`\delta P`$ reduce to the standard one-fluid equations (cf. Section II) and it seems reasonable to think that these variables should therefore reflect the ordinary fluid properties. At the same time, we have seen that $`\delta \beta `$ was to a certain extent accounted for in the standard picture via the buoyancy and the Schwarzschild discriminant $`A`$. Finally, the variable $`\stackrel{}{\xi }_{}`$ is clearly unique to the two-fluid system and could therefore be expected to bring some new features to the pulsation problem.

Let us first consider the case when the neutrons and protons move in such a way that they remain in chemical equilibrium. This would correspond to $`\delta \beta =0`$, and we immediately see that this requires $`\stackrel{}{\xi }_{}=0`$ as well. In other words, the neutrons and the protons must move together (quite intuitively). From the analysis in Section II (cf. the discussion that leads to (46)) we already know that the two equations (56) and (57) reduce to the standard equations for barotropic (nonstratified) stars ($`A=0`$). Hence, the only non-trivial modes we expect to find are the p-modes. In this case one can easily show that they will have frequency

$$\omega _o^2_l^2$$
(84)
noting, however, that $`c_s^2`$ is being replaced by $`c_{\mathrm{eq}}^2`$, i.e. the sound speed defined in (73). This is, of course, exactly the ordinary fluid result of Section II in the limit $`A0`$.

As a side remark it is interesting to note that one would not normally expect the above assumptions to lead to the two remaining equations, (58) and (62), also being satisfied. However, in the present case (58) is trivial and, for particular models, one can also satisfy (62). This requires

$$\frac{}{P}\left(\frac{\rho _p}{\rho }\right)_\beta =0.$$
(85)
In other words, we need to have $`\rho _p\rho `$ or $`\rho _p=0`$. The latter does, of course, correspond to a star in which there is only one species of particles so the equations of motion must reduce to the standard ones for a single fluid.

In Section II we saw that the nature of the g-modes could be deduced by requiring that the motion be such that the fluid remains in pressure equilibrium, $`\delta P=0`$. Let us consider the equations from the previous section in this case. Conveniently “forgetting” (for the moment) the conservation of mass equation (56) we need to consider

$$\frac{d\delta \beta }{dr}\frac{\omega ^2\rho }{r^2det\rho }Z=0,$$
(86)
and

$`{\displaystyle \frac{dZ}{dr}}`$ $`=`$ $`[({\displaystyle \frac{det\rho }{c_{eq}^4\mathrm{\Delta }}}_l^2+𝒩^2){\displaystyle \frac{1}{\omega ^2}}(1+\left[{\displaystyle \frac{c_{\mathrm{eq}}}{g}}𝒩\right]^2)]\times `$ (89)
$`{\displaystyle \frac{r^2c_{\mathrm{eq}}^2\mathrm{\Delta }}{\rho }}\delta \beta ,`$
which follows from (57)–(62), once we separate the radial and horizontal components of the equations (as in Section II) and introduce the new variable ($`\xi _{}^r`$ is the radial component of $`\stackrel{}{\xi }_{}`$)

$$Z=\frac{det\rho }{\rho }r^2\xi _{}^r,$$
(90)
as well as use the various definitions and relations given at the end of Section III. We now combine the two equations to find

$`{\displaystyle \frac{d}{dr}}\{{\displaystyle \frac{\rho }{r^2c_{\mathrm{eq}}^2\mathrm{\Delta }}}[(1+\left[{\displaystyle \frac{c_{\mathrm{eq}}}{g}}𝒩\right]^2)`$ (91)
(92)
$`({\displaystyle \frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}}_l^2+𝒩^2){\displaystyle \frac{1}{\omega ^2}}]^1{\displaystyle \frac{dZ}{dr}}\}+{\displaystyle \frac{\omega ^2\rho }{r^2det\rho }}Z=0`$ (93)
and it is easy to see that, just as in Section II, we have a Sturm-Liouville problem for large $`\omega ^2`$. Hence, we expect there to exist a set of modes for which $`\omega _n^2\mathrm{}`$ for the large overtones (with index $`n\mathrm{}`$).

To elucidate the nature of these modes, let us derive the appropriate dispersion relation. Assuming that the perturbation variables behave as $`\mathrm{exp}(ikr)`$ we readily arrive at

$$k^2=\frac{c_{\mathrm{eq}}^2\mathrm{\Delta }}{det\rho }\left[\left(1+\left[\frac{c_{\mathrm{eq}}}{g}𝒩\right]^2\right)\omega ^2\left(\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}_l^2+𝒩^2\right)\right]$$
(94)
and in the limit when $`l>>kr`$ we find

$$\omega _s^2\left(1+\left[\frac{c_{\mathrm{eq}}}{g}𝒩\right]^2\right)^1\left(\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}_l^2+𝒩^2\right).$$
(95)
Assuming that the sound speeds are well-behaved if $`x_p`$ is small, then we can infer from (73) that $`c_{\mathrm{eq}}^2`$, and hence $`_l^2`$, are roughly independent of the proton fraction to leading order. We also infer from (52) and (80) that the combination $`det\rho /\mathrm{\Delta }`$ is independent of the proton fraction to leading order. Finally, we see from (80) and (76) that

$$\frac{1}{\mathrm{\Delta }}\left(\frac{\rho }{\beta }\right)_P^2x_p,$$
(96)
which implies that

$$𝒩^2x_p.$$
(97)
Thus, to linear order in the proton fraction we see

$$\omega _s^2_l^2\left(1\left[\frac{c_{\mathrm{eq}}}{g}𝒩\right]^2\right)\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}+𝒩^2.$$
(98)
To the lowest-order in the proton fraction, the dominant contribution to $`\omega _s`$ is

$$\omega _s^2_l^2\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}\frac{m_p}{m_p^{}}\frac{l(l+1)}{r^2}c_p^2.$$
(99)
Hence we deduce that these “superfluid” pulsation modes are essentially governed by the proton sound speed $`c_p^2`$, and also that the entrainment effect plays an important role in determining their frequencies. We also notice that contributions due to $`𝒩`$ are negligible at this order. This approximate result provides a natural explanation for the result of Comer et al \[Comer et al 1999\] that the ordinary and superfluid frequencies are virtually the same when the adiabatic index of the neutrons (treated as a relativistic polytrope) is set equal to the index for the protons (also treated as a relativistic polytrope). In that case, the sound speeds of the two fluids are virtually the same. Furthermore, in the models considered in Comer et al \[Comer et al 1999\] there is no entrainment since $`m_p^{}`$ is set equal to $`m_p`$ .

### 4.2 The fully coupled case

Having established the existence of two distinct classes of pulsation modes in a superfluid star we now return to the full problem, and consider the modes of pulsation for the system of equations (56)-(62).

Proceeding in the now familiar way, we separate the radial and horizontal parts of (57) and (58), and then combine the results with the two conservation laws (56) and (62). This way we arrive at a set of four equations

$`{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{d}{dr}}\left(r^2\xi _+^r\right){\displaystyle \frac{g}{c_{\mathrm{eq}}^2}}\xi _+^r+\left[{\displaystyle \frac{1}{c_{\mathrm{eq}}^2}}{\displaystyle \frac{l(l+1)}{\omega ^2r^2}}\right]{\displaystyle \frac{\delta P}{\rho }}+`$ (100)
(101)
$`{\displaystyle \frac{1}{\rho }}\left({\displaystyle \frac{\rho }{\beta }}\right)_P\delta \beta =0,`$ (102)
(103)
$`{\displaystyle \frac{d}{dr}}\delta P+{\displaystyle \frac{g}{c_{\mathrm{eq}}^2}}\delta P\rho \omega ^2\xi _+^r+g\left({\displaystyle \frac{\rho }{\beta }}\right)_P\delta \beta =0,`$ (104)
(105)
$`Z{\displaystyle \frac{r^2}{\omega ^2}}{\displaystyle \frac{det\rho }{\rho }}{\displaystyle \frac{d}{dr}}\delta \beta =0,`$ (106)
(107)
$`{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{dZ}{dr}}+\left[{\displaystyle \frac{\rho _n^2}{\rho }}{\displaystyle \frac{}{\beta }}\left({\displaystyle \frac{\rho _p}{\rho _n}}\right){\displaystyle \frac{l(l+1)}{\omega ^2r^2}}{\displaystyle \frac{det\rho }{\rho }}\right]\delta \beta +`$ (108)
(109)
$`\left({\displaystyle \frac{\rho }{\beta }}\right)_P\left[{\displaystyle \frac{\delta P}{\rho }}g\xi _+^r\right]=0,`$ (110)
where $`\xi _+^r`$ is the radial component of $`\stackrel{}{\xi }_+`$ and the variable $`Z`$ was defined in Section IVA.

As in Section II the first two equations simplify considerably if we introduce the integrating factor $`\varphi `$ as defined by (21) (replacing $`c_s`$ with $`c_{\mathrm{eq}}`$, of course), and then work with the variables $`\delta \widehat{P}=\varphi \delta P`$ and $`\widehat{Y}=r^2\xi _+^r/\varphi `$. If we also introduce the characteristic wavelength of the pulsation through $`\mathrm{exp}(ikr)`$, our first two equations take the form

$$ik\widehat{Y}+\left[\frac{r^2}{c_{\mathrm{eq}}^2}\frac{l(l+1)}{\omega ^2}\right]\frac{\delta \widehat{P}}{\rho \varphi ^2}=\frac{r^2}{\rho }\left(\frac{\rho }{\beta }\right)_P\frac{\delta \beta }{\varphi }$$
(111)
and

$$ik\delta \widehat{P}\frac{\rho \omega ^2}{r^2}\varphi ^2\widehat{Y}=g\left(\frac{\rho }{\beta }\right)_P\varphi \delta \beta .$$
(112)
These two relations can easily be combined to provide a relation between $`\delta \widehat{P}`$ and $`\delta \beta `$:

$`\left[{\displaystyle \frac{l(l+1)+k^2r^2}{r^2}}c_{\mathrm{eq}}^2\omega ^2\right]\delta \widehat{P}`$ $`=`$ $`[ikg+\omega ^2]\times `$ (115)
$`c_{\mathrm{eq}}^2\varphi \left({\displaystyle \frac{\rho }{\beta }}\right)_P\delta \beta .`$
After some similar manipulations, one can show that the two remaining perturbation equations lead to

$`[\omega ^2ikg]\left({\displaystyle \frac{\rho }{\beta }}\right)_P\delta \widehat{P}=[{\displaystyle \frac{l(l+1)+k^2r^2}{r^2}}det\rho `$ (116)
(117)
$`\omega ^2\rho _n^2{\displaystyle \frac{}{\beta }}\left({\displaystyle \frac{\rho _p}{\rho _n}}\right)+g^2\left({\displaystyle \frac{\rho }{\beta }}\right)_P^2]\varphi \delta \beta .`$ (118)
These two equations can be combined to provide the dispersion relation for waves propagating in a superfluid star

$`[\omega ^4+k^2g^2]c_{\mathrm{eq}}^2\left({\displaystyle \frac{\rho }{\beta }}\right)_P^2=[{\displaystyle \frac{l(l+1)+k^2r^2}{r^2}}c_{\mathrm{eq}}^2\omega ^2]\times `$ (119)
(120)
$`[{\displaystyle \frac{l(l+1)+k^2r^2}{r^2}}det\rho \omega ^2\rho _n^2{\displaystyle \frac{}{\beta }}\left({\displaystyle \frac{\rho _p}{\rho _n}}\right)_P+`$ (121)
(122)
$`g^2\left({\displaystyle \frac{\rho }{\beta }}\right)_P^2].`$ (123)
In principle, this equation contains information equivalent to the dispersion relation of Lindblom and Mendell \[Lindblom & Mendell 1994\] (cf. their equation (82)).

We can do one final rewriting of (123) as a quadratic in $`\omega ^2`$, which we could then solve. Using the various definitions and relations listed at the end of Section III, the fully coupled dispersion relation takes the final form

$`0=\omega ^4\{[1+{\displaystyle \frac{k^2r^2}{l(l+1)}}][1+{\displaystyle \frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}}]_l^2+`$ (124)
(125)
$`[1+(1+{\displaystyle \frac{k^2r^2}{l(l+1)}})\left({\displaystyle \frac{c_{\mathrm{eq}}}{g}}_l\right)^2]𝒩^2\}\omega ^2+`$ (126)
$`_l^2\left\{\left[1+{\displaystyle \frac{k^2r^2}{l(l+1)}}\right]^2{\displaystyle \frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}}_l^2+𝒩^2\right\}.`$ (127)
It is, however, not particularly instructive to write down the formal solution to this equation. It is much better to first simplify it somewhat. This provides a better insight into the relevant physics of the solution.

### 4.3 Results for a small proton fraction

Our aim now is to use the results of the previous section to infer how the presence of a superfluid affects the pulsation modes of a nonrotating neutron star core. In principle, we have drawn the main conclusions already. There will exist two more or less distinct classes of modes. One of these comprise the standard p-modes, while the other has unique properties (although see comments below) due to the presence of the superfluid. Both these families of modes are such that the eigenfrequencies $`\omega _n^2\mathrm{}`$ as the mode-number $`n\mathrm{}`$. In other words, the two sets of modes will be interlaced in the pulsation spectrum of the star. These are, however, only qualitative results and it would clearly be useful to make more quantitative estimates.

The key point to simplifying the fully coupled dispersion relation is that $`𝒩^2x_p`$ which means that the frequencies can easily be determined to linear order in the proton fraction. Doing this (and considering also the limit where $`l>>kr`$) yields the two solutions

$$\omega _o^2_l^2\left(1+\left[1\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}\right]^1\left[\frac{c_{\mathrm{eq}}}{g}𝒩\right]^2\right)$$
(128)
and

$$\omega _s^2\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}_l^2\left(1\left[1\frac{det\rho }{c_{\mathrm{eq}}^4\mathrm{\Delta }}\right]^1\left[\frac{c_{\mathrm{eq}}}{g}𝒩\right]^2\right)+𝒩^2.$$
(129)
To lowest-order in the proton fraction, we recover the ordinary, i.e. $`\omega _o^2`$, and superfluid, i.e. $`\omega _s^2`$, solutions respectively.

Let us now address one of the main questions that motivated the present investigation: What happens to the g-modes when the neutron star core becomes superfluid? We recall the discussion by Reisenegger and Goldreich \[Reisenegger & Goldreich 1992\] that showed that a non-superfluid neutron star will have a distinct set of g-modes whose nature is determined by the chemical composition gradient associated with the varying proton fraction. From Section II we know that these modes will have eigenfrequencies such that $`\omega _n^20`$ as $`n\mathrm{}`$. We have already shown that there will not be a family of modes with this character in the superfluid case. Thus, even though we clearly have a varying proton fraction, a superfluid neutron star core will have no pulsating g-modes, in accordance with the numerical results of Lee \[Lee 1995\]. This is part of the answer to the question posed above, and in order to provide the rest of the answer, we will need to consider the space of time-independent solutions to the linearized equations, i.e. the zero-frequency subspace. But before doing that, we want to emphasize some consequences of having a varying proton fraction.

The effect of a composition gradient will be analyzed as follows: In the ordinary fluid case one can estimate the Schwarzschild discriminant due to a variation in the proton fraction as (see Unno et al \[Unno et al 1989\] and the discussion in Lee \[Lee 1995\] that surrounds Lee’s equation (47))

$$A=\frac{1}{\rho }\left(\frac{\rho }{\rho _p}\right)_P\widehat{e}_r\rho _p$$
(130)
(assuming the star has zero temperature). This then leads to g-mode frequencies

$$\omega _g^2N^2=gA=\frac{g}{\rho }\left(\frac{\rho }{\rho _p}\right)_P\widehat{e}_r\rho _p,$$
(131)
or

$$\omega _g^2N^2g^2\left(\frac{\rho }{\rho _p}\right)_P\left(\frac{\rho _p}{P}\right)_\beta x_p\left(\frac{g}{c_p}\right)^2.$$
(132)
Now consider the dominant contribution to the superfluid mode frequency (129) when $`x_p`$ and $`c_p^2`$ are both small; it is simply

$$\omega _s^2𝒩^2x_p\left(\frac{g}{c_p}\right)^2.$$
(133)
From this we see that $`\omega _g^2N^2𝒩^2\omega _s^2`$ when the proton fraction and $`c_p^2`$ are small. Thus we conclude that, even though a superfluid neutron star core does not support a distinct class of pulsating g-modes, the buoyancy associated with internal composition gradients can be relevant. Of course, the effect depends crucially on the proton fraction and cannot be distinguished unless $`c_p^2`$ is small.

## 5 The zero-frequency subspace

We will now focus our attention on the perturbations of the superfluid equations that belong to the zero-frequency subspace. That is, after perturbing (and linearizing) the equations, we assume that the velocity perturbations, as well as $`\delta n_n`$ and $`\delta n_p`$, are time-independent (cf. \[Lockitch & Friedman 1999, Lockitch et al 2001\]) . We then get from (211) and (212):

$$0=_i\left(\rho _n\delta v_n^i\right),0=_i\left(\rho _p\delta v_p^i\right).$$
(134)
Meanwhile, the two linearized Euler equations and the equation for the linearized Newtonian potential lead to

$$0=_i(\delta \mathrm{\Phi }+\delta \stackrel{~}{\mu }_n),0=_i(\delta \mathrm{\Phi }+\delta \stackrel{~}{\mu }_p),$$
(135)
and

$$^2\delta \mathrm{\Phi }=4\pi G\left(\delta \rho _n+\delta \rho _p\right)=4\pi G\delta \rho .$$
(136)
From these equations we immediately deduce that we can split the perturbations into two independent sets. The first set has the perturbed velocities vanishing but the other perturbations non-zero. The second set has $`\delta v_n^i`$ and $`\delta v_p^i`$ nonvanishing while the local matter perturbations $`\delta \rho _n=\delta \rho _p=\delta \stackrel{~}{\mu }_n=\delta \stackrel{~}{\mu }_p=\delta \mathrm{\Phi }=0`$. The general solution will be a superposition of solutions from each set.

### 5.1 Neighbouring spherical equilibria

We begin by interpreting the meaning of the first set of perturbations distinguished above. By taking the difference between the two equations in (135) it is easy to show that we must have

$$_i\delta \beta =0.$$
(137)
Given this, we can also show that the sum of the two equations leads to

$$_i\delta P=\delta \rho _i\mathrm{\Phi }\rho _i\delta \mathrm{\Phi }.$$
(138)
These two equations follow simply from the perturbed versions of the equations that determine the unperturbed configuration, i.e.

$$_i(\mathrm{\Phi }+\stackrel{~}{\mu }_n)=0,_i(\mathrm{\Phi }+\stackrel{~}{\mu }_p)=0,$$
(139)
using also the earlier equation that relates pressure variations to chemical potential variations. And just as for the equilibrium equations we can show that the solutions must have spherical symmetry. Thus, the interpretation of this set of perturbations becomes clear: They simply correspond to a neighbouring spherical equilibrium model. Provided that $`\delta \beta =0`$ the perturbed model is in chemical equilibrium.

### 5.2 Non-trivial velocity perturbations

We now turn to the other set of perturbations, for which the three-velocities are non-zero and represent convective currents. We first have to account for the fact that perturbations of a spherical star can be decoupled into two different classes: spheroidal and toroidal.

Spheroidal velocity perturbations take the form

$$\delta \stackrel{}{v}_j=\frac{1}{r}(W_j(r),V_j(r)_\theta ,\frac{V_j(r)}{\mathrm{sin}\theta }_\phi )Y_{lm},j=n,p$$
(140)
where we have used a standard decomposition in terms of spherical harmonics. With these definitions we readily get from (134)

$`0`$ $`=`$ $`l(l+1)n_nV_n{\displaystyle \frac{d}{dr}}\left(rn_nW_n\right),`$ (141)
$`0`$ $`=`$ $`l(l+1)n_pV_p{\displaystyle \frac{d}{dr}}\left(rn_pW_p\right).`$ (143)
Foregoing here a detailed discussion of boundary conditions, we can establish that non-trivial solutions exist. In particular, we see that there are two functions ($`W_n`$ and $`W_p`$, say) that can be freely specified. Because they are spheroidal, and have vanishing pressure and chemical potential perturbations, we identify these solutions as the g-modes “missing” from the time-dependent perturbations. In this respect, we see that superfluid g-modes behave like those in the barotropic ordinary fluid case, i.e. when the equation of state of the perturbations is the same as that of the background \[Lockitch & Friedman 1999, Lockitch et al 2001\]. However, the doubling of the modes (because there are two free functions) is unique to the superfluid.

Toroidal three-velocity perturbations can be written as

$$\delta \stackrel{}{v}_j=\frac{1}{r}(0,\frac{U_j(r)}{\mathrm{sin}\theta }_\phi ,U_j(r)_\theta )Y_{lm},j=n,p$$
(144)
For velocity perturbations of this form Eq. (134) is automatically satisfied for arbitrary $`U_n`$ and $`U_p`$. Since there are two functions that can be freely specified we deduce that there will be two sets of toroidal modes in the zero-frequency subspace. When the star is rotating this should lead to the presence of two classes of r-modes. We will verify this in the next section.

## 6 Rotating superfluid stars: Qualitative insights into the r-modes

In the last two years various aspects of rotating neutron stars have been under intense scrutiny following the discovery \[Andersson 1998, Friedman & Morsink 1998\] that the emission of gravitational waves drives the so-called r-modes unstable. Initial studies of the problem indicated that the r-mode instability might cause a newly born neutron star to spin down to a rotation rate comparable to that inferred from observation for the Crab pulsar \[Lindblom et al 1998, Andersson et al 1999\], in the process radiating gravitational waves that may well prove detectable with the generation of large-scale interferometers due to come online in the next few years \[Owen et al 1998\]. One issue of utmost importance for the r-mode problem concerns the role of superfluidity \[Andersson & Kokkotas 2001\]. So for example was it originally thought that dissipation due to the superfluid mutual friction would counteract the instability in a significant way. The only available calculation of this effect, due to Lindblom and Mendell \[Lindblom & Mendell 2000\], suggests a rather different picture, however. Their results indicate that the mutual friction will typically not be able to suppress the r-mode instability. But the results also suggest that the effect becomes dominant for certain values of the entrainment parameter $`ϵ`$. Clearly, this problem is far from well understood at the present time.

Our aim in this section is to study the r-modes at a qualitative level (comparable to our study of spheroidal p- and g-modes in the previous sections). This would seem a natural starting point for a discussion of the r-modes of superfluid stars, and as we will see it provides insights that may well explain various features observed in the numerical work by Lindblom and Mendell \[Lindblom & Mendell 2000\].

To analyze the r-mode problem we will extend Saio’s \[Saio 1982\] vorticity argument to the two-fluid problem. Thus we focus on the two Euler equations. When expressed in terms of the variables $`\stackrel{}{\xi }_+`$ and $`\stackrel{}{\xi }_{}`$, these equations can be written

$`_t^2\stackrel{}{\xi }_++(\stackrel{}{u})_t\stackrel{}{\xi }_++_t\stackrel{}{\xi }_+\stackrel{}{u}+\delta \mathrm{\Phi }+`$ (145)
(146)
$`\delta \stackrel{~}{\mu }_n+x_p\delta \beta =0`$ (147)
and

$`_t^2\stackrel{}{\xi }_{}+(\stackrel{}{u})_t\stackrel{}{\xi }_{}_t\stackrel{}{\xi }_{}\stackrel{}{u}+\delta \beta +`$ (148)
(149)
$`2{\displaystyle \frac{det\rho }{\rho _n\rho _p}}(_t\stackrel{}{\xi }_{})\stackrel{}{u}=0.`$ (150)
Here we have assumed that the protons and neutrons corotate in the unperturbed case. The relevant rotation velocity is denoted by $`\stackrel{}{u}`$ in the above equations. We note that it may be desirable to relax the assumption of corotation of the two background fluids in order to model a real neutron star (cf. the discussion in Andersson and Comer \[Andersson & Comer 2001\]). The resultant problem is, however, much more complicated than the present one and we will return to it in future investigations.

We first translate these equations into the rotating frame. Then we get

$$_t^2\stackrel{}{\xi }_++2\stackrel{}{\mathrm{\Omega }}\times _t\stackrel{}{\xi }_++\delta \mathrm{\Phi }+\delta \stackrel{~}{\mu }_n+x_p\delta \beta =0$$
(151)
and

$`_t^2\stackrel{}{\xi }_{}+2\stackrel{}{\mathrm{\Omega }}\times _t\stackrel{}{\xi }_{}+\delta \beta 2[1{\displaystyle \frac{det\rho }{\rho _n\rho _p}}]\times `$ (152)
(153)
$`(_t\stackrel{}{\xi }_{})\stackrel{}{u}=0.`$ (154)
The next step involves assuming that the mode is horizontal to leading order, taking the curl of the above two equations and focusing on the radial component of the resultant relations. Consider first (151), which readily yields

$$_t\left[\times _t\stackrel{}{\xi }_++2\times \stackrel{}{\mathrm{\Omega }}\times \stackrel{}{\xi }_+\right]_r=0$$
(155)
since

$$\times \left(x_p\delta \beta \right)=x_p\times \delta \beta $$
(156)
has a vanishing radial component for a slowly rotating (still spherical) star. We now use

$`\left[\times \stackrel{}{\mathrm{\Omega }}\times \stackrel{}{\xi }_+\right]_r`$ $`=`$ $`\{\stackrel{}{\mathrm{\Omega }}(\stackrel{}{\xi }_+)\stackrel{}{\xi }_+(\stackrel{}{\mathrm{\Omega }})+`$ (159)
$`(\stackrel{}{\xi }_+)\stackrel{}{\mathrm{\Omega }}(\stackrel{}{\mathrm{\Omega }})\stackrel{}{\xi }_+\}_r`$
$``$ $`\left[(\stackrel{}{\xi }_+)\stackrel{}{\mathrm{\Omega }}\right]_r.`$ (161)
After these manipulations we have arrived at

$$_t\left[_t(\times \stackrel{}{\xi }_+)+2(\stackrel{}{\xi }_+)\stackrel{}{\mathrm{\Omega }}\right]_r0.$$
(162)
Finally, we insert in this equation the assumption that the mode we are interested in is purely toroidal, i.e. is of the form .

$$\stackrel{}{\xi }_+=r(0,\frac{T}{\mathrm{sin}\theta }_\phi ,T_\theta )e^{i\omega _rt}Y_{lm},$$
(163)
Given this assumption we immediately find that these modes must have frequency

$$\omega _r=\frac{2m\mathrm{\Omega }}{l(l+1)},$$
(164)
in the rotating frame. This is, of course, the standard r-mode result \[Andersson & Kokkotas 2001\].

Let us now go through the same exercise for the second Euler equation, (154), that appears in the superfluid case. The radial component of the curl of this equation can be written (in the rotating frame)

$`_t\{(\times _t\stackrel{}{\xi }_{})+2\times (\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{\xi }_{})`$ (165)
(166)
$`2\times \left[(1{\displaystyle \frac{det\rho }{\rho _n\rho _p}})(\stackrel{}{\xi }_{})\stackrel{}{u}\right]\}_r=0.`$ (167)
Here the first two terms are identical to those we encountered in the analysis of (151), but the last term is new. This term can be written

$$2\times \left(\gamma \stackrel{}{\xi }_{}\stackrel{}{u}\right)=2\left(\gamma \right)\times (\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{\xi }_{})+2\gamma \times (\stackrel{}{\mathrm{\Omega }}\times \stackrel{}{\xi }_{}),$$
(168)
where we have introduced

$$\gamma =1\frac{det\rho }{\rho _n\rho _p}.$$
(169)
It is clear that the first term of the right-hand side of (168) does not contribute to the radial component. The second term, however, does have a non-vanishing radial component. We readily find that

$`[2\times (\gamma \stackrel{}{\xi }_{}\stackrel{}{u})]_r={\displaystyle \frac{2\gamma \mathrm{\Omega }}{r\mathrm{sin}\theta }}\times `$ (170)
(171)
$`\left[_\theta (\mathrm{sin}\theta \mathrm{cos}\theta \xi _{}^\theta )+_\phi (\mathrm{cos}\theta \xi _{}^\phi )\right]`$ (172)
and if we assume that the vector $`\stackrel{}{\xi }_{}`$ is toroidal, i.e. takes the form (163), we have

$`[2\times \left(\gamma \stackrel{}{\xi }_{}\stackrel{}{u}\right)]_r`$ $`=`$ $`{\displaystyle \frac{2\gamma \mathrm{\Omega }T}{\mathrm{sin}\theta }}\times `$ (175)
$`\left[_\theta (\mathrm{cos}\theta _\phi Y_{lm})_\phi (\mathrm{cos}\theta _\theta Y_{lm})\right]`$
$`=`$ $`2im\mathrm{\Omega }\gamma TY_{lm}.`$ (177)
We can now combine this result with the final result obtained from (151) to deduce that we will have modes with frequency

$$\omega _r=\frac{2m\mathrm{\Omega }}{l(l+1)}\frac{det\rho }{\rho _n\rho _p}.$$
(178)
From the definition of $`det\rho `$ it follows that this is identical to

$$\omega _r=\frac{2m\mathrm{\Omega }}{l(l+1)}\left[1+ϵ\left(1+\frac{\rho _n}{\rho _p}\right)\right],$$
(179)
where $`ϵ`$ is the entrainment coefficient. This is a very interesting result since it demonstrates the existence of a distinct class of superfluid r-modes, the properties of which are to a large extent determined by the entrainment parameter.

Based on the above analysis we can now discuss the general nature of the toroidal modes of rotating superfluid neutron star cores (even though we should advice some caution since (179) is $`r`$-dependent and therefore strictly speaking only describes an r-mode in a thin shell). We have seen that, just like in the non-rotating case, there will be two distinct classes of modes. The first corresponds to fluid motion such that the neutrons and the protons flow together, and the resultant modes are analogous to the standard perfect fluid r-modes. These are the modes that the calculations of Lindblom and Mendell \[Lindblom & Mendell 2000\] were aimed at studying. For the second class, the neutrons and the protons are counter-moving. This class of r-modes has not previously been discussed in the literature (although see the discussion below). It is interesting to note that the two classes of modes are (to leading order) degenerate in the absence of entrainment. This is, of course, rather intuitive since the two fluids are then effectively uncoupled and the Coriolis force affects each fluid separately.

We note that the existence of two distinct classes of r-modes in superfluid stars provides a likely explanation for the somewhat peculiar results obtained by Lindblom and Mendell \[Lindblom & Mendell 2000\] regarding the effect of mutual friction on the r-mode instability. In their study Lindblom and Mendell found that the mutual friction dissipation timescale was sufficiently long that this effect would not suppress the unstable r-modes for most values of $`ϵ`$. But they also discovered that mutual friction would become dominant for certain values of $`ϵ`$, cf. their figure 6. Furthermore, they noted that while the variable $`\delta \beta `$ was small in most cases, it became large for the particular values of $`ϵ`$ at which mutual friction was found to be important. Given our present conclusions regarding the r-modes in superfluid stars, we can interpret the Lindblom-Mendell result in the following way: There are two distinct classes of r-modes in a superfluid star. One of these classes, the one for which the neutrons and protons flow together, is typically not rapidly damped by mutual friction. The other class of modes, however, is such that the two fluids are counter-moving. Since mutual friction tends to damp relative motion between the two fluids these modes will be strongly affected by mutual friction. This then makes the Lindblom-Mendell results quite natural. For certain values of the entrainment parameter they have simply found the second class of modes rather than the first.

The situation may, however, be yet more complicated. It is likely relevant to make the following observation: In many situations where the mode-spectrum is studied for a family of stellar models for which a single parameter is varied (such as the crustal shear modulus) one finds that the modes undergo what are known as “avoided crossings.” So for example have recent studies of neutron stars with an elastic crust shown that the r-modes associated with the core fluid show avoided crossings with the toroidal shear modes in the crust for certain values of the rotation rate. Avoided crossings between pulsation modes are, in fact, a regular occurence in studies of more complicated stellar models \[Cox 1980, Unno et al 1989\]. It is interesting to note that, well away from the crossing point, the different modes have a distinct nature, but as one approaches the avoided crossings the two modes become similar in nature and when one reaches “the other side” they have completely exchanged properties (the core r-mode has become a crust shear mode and vice versa). In view of our qualitative analysis and the numerical results of Lindblom and Mendell, it would seem natural to predict that a more detailed investigation of the superfluid r-mode problem will unveil analogous avoided crossings between the co- and counter-moving r-modes as the entrainment coefficient is varied.

## 7 Conclusions

In this paper we have compared and contrasted the pulsation properties of a superfluid neutron star core (represented by a relatively simple two-fluid model) with the familiar (textbook) results for a normal fluid star. This study provides a qualitative understanding of the nature of the oscillation modes of a superfluid neutron star, and thus provides a theoretical fundament for various numerical results in this area \[Lindblom & Mendell 1994, Lee 1995\]. We have shown that a non-rotating superfluid neutron star core exhibits two distinct families of pulsation modes. Our approximate results provide a deep insight into the physics of these modes, with one set of modes corresponding to the two fluids moving in “lock-step” while the second family (the modes that are unique to a superfluid) correspond to the two fluids being counter-moving. We have shown that both sets of modes are crucially governed by the acoustic properties of the two fluids. This is in clear contrast to the normal fluid case, where the two families of modes, the p- and g-modes, are acoustic and governed by buoyancy (mainly due to chemical composition gradients), respectively. We have investigated what happens to the g-modes as the star becomes superfluid and have confirmed Lee’s numerical results that there are no propagating g-modes in a superfluid core. In addition, we have shown that the “missing” g-modes can be found in the zero-frequency subspace. Finally, we have shown how the various modes are affected by the parameters of entrainment.

Basically, our results are important for two reasons. Firstly, they provide a clear and concise description of the problem and the general nature of pulsating superfluid stars. Thus they fill what we perceive as a gap in the existing literature. In particular, we believe we have resolved some open questions regarding g-modes in the superfluid case. Secondly, the results indicate that future observations of neutron star mode oscillations, eg. by a highly sensitive generation of gravitational wave detectors, could potentially help constrain the parameters of the large scale superfluidity that is believed to exist in the core of mature neutron stars. We will discuss this exciting possibility elsewhere.

In addition to studying non-rotating stellar cores, we have presented a qualitative analysis of the r-modes in a slowly rotating superfluid core. We have demonstrated that there will be two distinct families of such modes, and argued that the interplay between these modes (eg. via so-called avoided crossings as the entrainment parameter is varied) may provide an understanding of puzzling numerical results obtained by Lindblom and Mendell \[Lindblom & Mendell 2000\]. Our current analysis clearly shows that the r-mode problem for superfluid stars is likely to be much richer than has previously been appreciated. In particular since, in addition to the r-modes, we anticipate that there will be a large class of “hybrid” inertial modes \[Lockitch & Friedman 1999, Lockitch et al 2001\]. No attempts to study such modes in a superfluid star have yet been made. A better understanding of these various issues is clearly needed in view of the fact that the r-modes have been shown to be unstable due to the emission of gravitational waves, and that this instability may govern the spin-evolution of neutron stars during various phases of their lives. We therefore plan to investigate this problem in greater detail in the near future.

## Acknowledgements

We wish to thank J. Friedman, K. D. Kokkotas, D. Langlois, and R. Prix for fruitful discussions. GLC gratefully acknowledges partial support from a Saint Louis University SLU2000 Faculty Research Leave Award as well as EPSRC in the UK via grant number GR/R52169/01, and the warm hospitality of the Center for Gravitation and Cosmology of the University of Wisconsin-Milwaukee and the University of Southampton where part of this research was carried out. NA acknowledges generous support from PPARC via grant numberPPA/G/S/1998/00606, the EU Network “Sources of Gravitational Radiation,” as well as the EU programme “Improving the Human Research Potential and the Socio-Economic Knowledge Base” (research training network contract HPRN-CT-2000-00137).

## Appendix A Newtonian limit of the general relativistic superfluid equations

The purpose of this Appendix is to derive the Newtonian version of the general relativistic superfluid field equations. There are several motivations for doing this. For conceptual reasons we want to be able to compare our relativistic calculations \[Comer et al 1999, Andersson et al 2001\] to previous results in this research area (in particular those of Lindblom and Mendell \[Lindblom & Mendell 1994, Lindblom & Mendell 1995, Lindblom & Mendell 2000\]) most of which were obtained in the Newtonian regime. We also feel that it is useful to address various problems on a qualitative level before working out the detailed answers. Given that Newtonian calculations are often simpler, it is natural to address the Newtonian problem first. In addition to this, we hope that a study of the Newtonian limit will provide insight into ways of extending various Newtonian entrainment models to the relativistic regime.

### A.1 The general relativistic formalism and Newtonian limit

We begin by recounting the formalism that has been used to model general relativistic superfluid neutron stars \[Carter & Langlois 1995a, Carter & Langlois 1995b, Langlois et al 1998, Comer et al 1999, Andersson & Comer 2001\]. The fluid dynamical degrees of freedom are described by $`n_n^\mu `$, the conserved neutron number density current, and $`n_p^\mu `$, the conserved proton number density current. The fundamental scalar in the formalism is the so-called “master” function $`\mathrm{\Lambda }`$, which is a function of the three scalars $`n_n^2=n_\rho ^nn_n^\rho `$, $`n_p^2=n_\rho ^pn_p^\rho `$, and $`x^2=n_\rho ^pn_n^\rho `$, where $`n_\mu ^n=g_{\mu \nu }n_n^\nu `$ and $`n_\mu ^p=g_{\mu \nu }n_p^\nu `$. The quantity $`\mathrm{\Lambda }`$ corresponds to the total thermodynamic energy density of the entire fluid.

There are two momenta $`\mu _\mu `$ and $`\chi _\mu `$ which are dynamically, and thermodynamically, conjugate to $`n_n^\mu `$ and $`n_p^\mu `$. They are defined by a general variation (that keeps the spacetime metric fixed) of $`\mathrm{\Lambda }(n_n^2,n_p^2,x^2)`$, i.e.

$$\delta \mathrm{\Lambda }=\mu _\rho \delta n_n^\rho +\chi _\rho \delta n_p^\rho ,$$
(180)
where

$$\mu _\mu =𝒦^{nn}n_\mu ^n+𝒦^{np}n_\mu ^p,\chi _\mu =𝒦^{pp}n_\mu ^p+𝒦^{np}n_\mu ^n,$$
(181)
and

$$𝒦^{np}=\frac{\mathrm{\Lambda }}{x^2},𝒦^{nn}=2\frac{\mathrm{\Lambda }}{n_n^2},𝒦^{pp}=2\frac{\mathrm{\Lambda }}{n_p^2}.$$
(182)
In this general relativistic context the entrainment effect is seen in that the momentum $`\mu _\mu `$, say, is a linear combination of $`n_n^\mu `$ and $`n_p^\mu `$. The generalized pressure $`\mathrm{\Psi }`$ is given by

$$\mathrm{\Psi }=\mathrm{\Lambda }n_n^\rho \mu _\rho n_p^\rho \chi _\rho .$$
(183)
The equations of motion consist of two conservation equations,

$$_\mu n_n^\mu =0,_\mu n_p^\mu =0,$$
(184)
and two Euler type equations, which can be conveniently written in the compact form

$$n_n^\mu _{[\mu }\mu _{\nu ]}=0,n_p^\mu _{[\mu }\chi _{\nu ]}=0,$$
(185)
where the square braces ‘\[ \]’ indicate antisymmetrization on the indices.

In order to derive the Newtonian limit, the general relativistic field equations will be written to order $`c^0`$ where $`c`$ is the speed of light. Formally, the Newtonian equations will then be obtained in the limit that the speed of light $`c`$ becomes infinite. The gravitational potential, denoted $`\mathrm{\Phi }`$, is assumed to be small in the sense that

$$1<<\frac{\mathrm{\Phi }}{c^2}0.$$
(186)
To order $`c^0`$ the metric can be written as

$$ds^2=c^2\left(1+\frac{2\mathrm{\Phi }}{c^2}\right)dt^2+\delta _{ij}dx^idx^j,$$
(187)
where the $`x^i`$ ($`i=1,2,3`$) are Cartesian-like coordinates.

In the superfluid field equations there are two four-velocities that must be considered. This means that we must take into account the fact that the fluids define two different proper times: one for the neutrons, to be denoted $`\tau _n`$, and the second for the protons, to be denoted $`\tau _p`$. The two different fluid trajectories are then obtained from the functions

$`x_n^\mu (\tau _n)`$ $`=`$ $`(t(\tau _n),x_n^i(\tau _n)),`$ (188)
$`x_p^\mu (\tau _p)`$ $`=`$ $`(t(\tau _p),x_p^i(\tau _p)).`$ (190)
The respective four-velocities of the neutron and proton fluids are thus given by

$$u_n^\mu =\frac{dx_n^\mu }{d\tau _n},u_p^\mu =\frac{dx_p^\mu }{d\tau _p}.$$
(191)
Because of the choice of coordinates, the four-velocities satisfy $`u_\mu ^nu_n^\mu =c^2`$ and $`u_\mu ^pu_p^\mu =c^2`$, where $`u_\mu ^n=g_{\mu \nu }u_n^\nu `$ and $`u_\mu ^p=g_{\mu \nu }u_p^\nu `$.

We will turn this around and use a global time $`t(\tau _n)=t(\tau _p)t`$ as the parameter for both curves. In this case, the two proper times will be given by

$`ds_n^2`$ $`=`$ $`c^2d\tau _n^2=c^2\left(1+{\displaystyle \frac{2\mathrm{\Phi }}{c^2}}{\displaystyle \frac{\delta _{ij}v_n^iv_n^j}{c^2}}\right)dt^2,`$ (192)
$`ds_p^2`$ $`=`$ $`c^2d\tau _p^2=c^2\left(1+{\displaystyle \frac{2\mathrm{\Phi }}{c^2}}{\displaystyle \frac{\delta _{ij}v_p^iv_p^j}{c^2}}\right)dt^2,`$ (194)
where $`v_n^i=dx_n^i/dt`$ and $`v_p^i=dx_p^i/dt`$ are the Newtonian three-velocities of the neutron and proton fluids, respectively. Each three-velocity is considered to be small in the sense that

$$\frac{\left|v_n^i\right|}{c}<<1,\frac{\left|v_p^i\right|}{c}<<1.$$
(195)
Hence, to the correct order the four-velocity components are given by

$$u_n^t=1\frac{\mathrm{\Phi }}{c^2}+\frac{v_n^2}{2c^2},u_n^i=v_n^i,$$
(196)
and

$$u_p^t=1\frac{\mathrm{\Phi }}{c^2}+\frac{v_p^2}{2c^2},u_p^i=v_p^i,$$
(197)
where $`v_n^2=\delta _{ij}v_n^iv_n^j`$ and $`v_p^2=\delta _{ij}v_p^iv_p^j`$.

Note that the two particle number currents are now written as

$$n_n^\mu =n_n(u_n^\mu /c),n_p^\mu =n_p(u_p^\mu /c).$$
(198)
To the correct order, one finds for $`x^2=n_\mu ^nn_p^\mu `$ that

$$x^2=n_nn_p\left(1+\frac{w^2}{2c^2}\right),$$
(199)
where

$$w^2=\delta _{ij}\left(v_n^iv_p^i\right)\left(v_n^jv_p^j\right).$$
(200)
In order to write the Euler equations with the terms to the required order, it is necessary to explicitly break up the “master” function into its mass part and internal energy part $`E`$, i.e., to write $`\mathrm{\Lambda }`$ as

$$\mathrm{\Lambda }=\left(m_nn_n+m_pn_p\right)c^2E(n_n^2,n_p^2,x^2).$$
(201)
In this context, $`E`$ is small in the sense that, for instance,

$$0\frac{E}{m_nn_nc^2}<<1.$$
(202)
It is also convenient to use a different choice for the independent variables that more closely agrees with what is used for Newtonian superfluids, which is the triplet of variables $`(n_n^2,n_p^2,w^2)`$. Thus, from now on we assume that $`E=E(n_n^2,n_p^2,w^2)`$.

With this choice, a variation of $`\mathrm{\Lambda }`$ that leaves the metric fixed yields the following:

$$d\mathrm{\Lambda }=(m_nc^2+\mu _n)dn_n(m_pc^2+\mu _p)dn_p\alpha dw^2,$$
(203)
where

$$\mu _n=\frac{E}{n_n},\mu _p=\frac{E}{n_p},\alpha =\frac{E}{w^2}.$$
(204)
Now the $`𝒦^{np}`$, $`𝒦^{nn}`$, and $`𝒦^{pp}`$ coefficients of the general relativistic formalism are related to the coefficients defined above via

$`𝒦^{np}`$ $`=`$ $`{\displaystyle \frac{2c^2}{n_nn_p}}\alpha ,`$ (205)
$`𝒦^{nn}`$ $`=`$ $`{\displaystyle \frac{m_nc^2+\mu _n}{n_n}}{\displaystyle \frac{2c^2}{n_n^2}}\left(1+{\displaystyle \frac{w^2}{2c^2}}\right)\alpha ,`$ (207)
$`𝒦^{pp}`$ $`=`$ $`{\displaystyle \frac{m_pc^2+\mu _p}{n_p}}{\displaystyle \frac{2c^2}{n_p^2}}\left(1+{\displaystyle \frac{w^2}{2c^2}}\right)\alpha .`$ (209)
In terms of these variables, the pressure $`P`$ is seen to be

$$P=E+\mu _nn_n+\mu _pn_p.$$
(210)
The Newtonian limit of the general relativistic superfluid field equations reduce to the following set of 8 equations:

$$0=\frac{n_n}{t}+_i\left(n_nv_n^i\right),$$
(211)
$$0=\frac{n_p}{t}+_i\left(n_pv_p^i\right),$$
(212)
$`0`$ $`=`$ $`{\displaystyle \frac{}{t}}\left(v_n^i+{\displaystyle \frac{2\alpha }{m_nn_n}}\left[v_p^iv_n^i\right]\right)+`$ (217)
$`v_n^j_j\left(v_n^i+{\displaystyle \frac{2\alpha }{m_nn_n}}\left[v_p^iv_n^i\right]\right)+\delta ^{ij}_j\left(\mathrm{\Phi }+{\displaystyle \frac{\mu _n}{m_n}}\right)+`$
$`{\displaystyle \frac{2\alpha }{m_nn_n}}\delta ^{ij}\delta _{kl}\left(v_p^lv_n^l\right)_jv_n^k,`$
and

$`0`$ $`=`$ $`{\displaystyle \frac{}{t}}\left(v_p^i+{\displaystyle \frac{2\alpha }{m_pn_p}}\left[v_n^iv_p^i\right]\right)+`$ (222)
$`v_p^j_j\left(v_p^i+{\displaystyle \frac{2\alpha }{m_pn_p}}\left[v_n^iv_p^i\right]\right)+\delta ^{ij}_j\left(\mathrm{\Phi }+{\displaystyle \frac{\mu _p}{m_p}}\right)+`$
$`{\displaystyle \frac{2\alpha }{m_pn_p}}\delta ^{ij}\delta _{kl}\left(v_n^lv_p^l\right)_jv_p^k.`$
Finally, the equation for the gravitational potential $`\mathrm{\Phi }`$ is

$$_i^i\mathrm{\Phi }=4\pi G\left(m_nn_n+m_pn_p\right).$$
(223)
These equations have been derived independently by Prix \[Prix 2001\], using a Newtonian Lagrangian-based variational principle.

### A.2 Comparison with the Formalism of Lindblom and Mendell

In their studies of oscillating superfluid stars Lindblom and Mendell \[Mendell & Lindblom 1991, Mendell 1991, Lindblom & Mendell 1994, Lindblom & Mendell 1995, Lindblom & Mendell 2000\] use the mass densities, $`\rho _n`$ and $`\rho _p`$, and the fluid three-“velocities” $`\stackrel{}{V}_n`$ and $`\stackrel{}{V}_p`$ as their main variables. The latter are not the same as the three-velocities used in the main text of this paper, rather these “velocities” are the “macroscopically averaged” gradients of the phases of the mesoscopic wave functions that describe the superfluid neutrons and the superconducting protons and they are actually proportional to the spatial components of our neutron and proton momenta. By comparing our mass-conservation equations with those of Lindblom and Mendell, which are given by

$`0`$ $`=`$ $`{\displaystyle \frac{\rho _n}{t}}+\stackrel{}{}\left(\rho _{nn}\stackrel{}{V}_n+\rho _{np}\stackrel{}{V}_p\right),`$ (224)
$`0`$ $`=`$ $`{\displaystyle \frac{\rho _p}{t}}+\stackrel{}{}\left(\rho _{pn}\stackrel{}{V}_n+\rho _{pp}\stackrel{}{V}_p\right),`$ (226)
where $`\rho _{nn}`$, $`\rho _{pp}`$ and $`\rho _{np}=\rho _{pn}`$ were defined earlier in the main text, we find that

$`\rho _n\stackrel{}{v}_n`$ $`=`$ $`\rho _{nn}\stackrel{}{V}_n+\rho _{np}\stackrel{}{V}_p,`$ (227)
$`\rho _p\stackrel{}{v}_p`$ $`=`$ $`\rho _{pn}\stackrel{}{V}_n+\rho _{pp}\stackrel{}{V}_p.`$ (229)
A final comparison to our fluid momenta yields

$$\alpha =\frac{1}{2}\frac{\rho _n\rho _p}{\rho _{np}^2\rho _{nn}\rho _{pp}}\rho _{np}.$$
(230)
This is a very useful result since it enables us to make contact with the model that Lindblom and Mendell use to describe the entrainment effect.

In Section III of the main text, we made reference to equations (69)-(72) of Lindblom and Mendell \[Lindblom & Mendell 1994\]. Neglecting terms of order $`(m_e/m_n)^2`$, where $`m_e`$ is the electron mass, then these equations are (recalling that $`\stackrel{~}{\mu }_n=\mu _n/m_n`$ etcetera)

$`\left({\displaystyle \frac{\rho }{P}}\right)_\beta `$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{\rho }}[\left({\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _n}}\right)_{\rho _p}2\left({\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _p}}\right)_{\rho _n}+`$ (233)
$`\left({\displaystyle \frac{\stackrel{~}{\mu }_p}{\rho _p}}\right)_{\rho _n}],`$
$`\left({\displaystyle \frac{\rho }{\beta }}\right)_P`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{\rho }}[\rho _n\left({\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _n}}\right)_{\rho _p}\rho _p\left({\displaystyle \frac{\stackrel{~}{\mu }_p}{\rho _p}}\right)_{\rho _n}+`$ (236)
$`(\rho _p\rho _n)\left({\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _p}}\right)_{\rho _n}],`$
$`{\displaystyle \frac{}{\beta }}\left({\displaystyle \frac{\rho _p}{\rho _n}}\right)_P`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{\rho }}[\left({\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _n}}\right)_{\rho _p}+2{\displaystyle \frac{\rho _p}{\rho _n}}\left({\displaystyle \frac{\stackrel{~}{\mu }_n}{\rho _p}}\right)_{\rho _n}+`$ (239)
$`{\displaystyle \frac{\rho _p^2}{\rho _n^2}}\left({\displaystyle \frac{\stackrel{~}{\mu }_p}{\rho _p}}\right)_{\rho _n}],`$
and

$$\frac{1}{\mathrm{\Delta }}=\left(\frac{\stackrel{~}{\mu }_n}{\rho _n}\right)_{\rho _p}\left(\frac{\stackrel{~}{\mu }_p}{\rho _p}\right)_{\rho _n}\left(\frac{\stackrel{~}{\mu }_n}{\rho _p}\right)_{\rho _n}^2.$$
(240)
Note that these follow by considering $`\mu _n`$ and $`\mu _p`$ to be functions of $`\rho _n`$ and $`\rho _p`$, and thus re-writing (34) and (35) in the forms

$$\delta P=(\mathrm{})\delta \rho _n+(\mathrm{})\delta \rho _p,\delta \beta =(\mathrm{})\delta \rho _n+(\mathrm{})\delta \rho _p,$$
(241)
and then inverting to find

$$\delta \rho _n=(\mathrm{})\delta P+(\mathrm{})\delta \beta ,\delta \rho _p=(\mathrm{})\delta P+(\mathrm{})\delta \beta .$$
(242)
Finally, by introducing the various sounds speeds we arrive at Eqns (73) - (80).