text:
# Adiabatic oscillations of non–rotating superfluid neutron stars

## 1 Introduction

The study of stellar oscillations has proved very fruitful in improving our understanding of the inner structure and dynamics of stars (the terms helio– and astro–seismology have been coined), for which the oscillation modes can often be observed rather directly. The best developed example of this probing of the internal structure of an astrophysical body via its oscillations is probably the Earth. In the case of neutron stars, the observation of oscillations is unfortunately not possible in such a direct way, and has not yet been achieved. In practically all cases we can only observe the regular radio–pulses of neutron stars, which are virtually unaffected by its oscillations and give information mostly about their rotation rate. Nevertheless this field bears great potential interest: on one hand the better understanding of neutron star oscillations could eventually help to elucidate the phenomenon of glitches, which is probably the most striking and puzzling aspect of observed neutron star dynamics. This phenomenon still represents somewhat of a mystery, even though the crucial role of superfluidity seems well established (see Link et al. (2000); Carter et al. (2000) for recent discussions). On the other hand, several highly sensitive gravitational wave detectors are expected to reach their full sensitivity within the next few years, and neutron star oscillations represent one of the potentially most interesting sources of gravitational waves. Gravitational wave detection could open a new and complementary observational window onto neutron stars, which would allow us to learn much more about their inner structure and dynamics than it is currently possible with the purely electro-magnetic observations.

Currently most studies of neutron star oscillations are based on simple perfect fluid models, which neglects the crucial importance of superfluidity in neutron stars. The presence of substantial amounts of superfluid matter in neutron stars is backed by a number of theoretical calculations of the state of matter at these extreme densities (e.g. see Baldo et al. (1992); Sjöberg (1976)), and by the qualitative success of superfluid models to accommodate observed features of glitches and their relaxation.

The first study to point out the importance of superfluidity for the oscillation properties of neutron stars was by Epstein (1988), who has argued in a local (sound wave) analysis that superfluidity should lead to new modes and modify the previously known modes. Lindblom & Mendell (1994) have argued further for the existence of these modes, but failed to find them numerically. Lee (1995) presented the first numerical results indicating the presence of new modes that did not exist in perfect fluid models, and the absence of g–modes which would have been present in the non–superfluid case. A local analysis by Andersson & Comer (2001a) has given further analytic evidence for the absence of g–modes in simple superfluid models. The relativistic numerical analysis by Comer et al. (1999) has shown an effective doubling of acoustic modes in superfluid models with respect to the normal fluid case. Recently this work has been extended by Andersson et al. (2002) to include entrainment, and they have shown that avoided mode crossings occur when one varies the entrainment parameter. The relevance of superfluid oscillations for gravitational wave detection has been discussed by Andersson & Comer (2001b). Some studies have also started to look at oscillations of rotating superfluid neutron stars (Lindblom & Mendell 2000) and Sedrakian & Wasserman (2000).

Despite the number of studies on oscillations of non–rotating superfluid neutrons stars, we think that this problem still deserves attention and that several points needed to be clarified. In particular it is worth emphasising the importance of stratification for the nature of superfluid oscillations, a point that has not yet been fully appreciated. We demonstrate here that only in non–stratified models can the eigenmode spectrum be separated into two families of modes, one of which is *identical* to the case of a normal–fluid star, while the other is characterised by counter–motion of the two fluids and vanishing gravitational perturbation. These two distinct families are usually referred to as “ordinary” and “superfluid” modes. Stratification of the background star, however, couples these distinct mode-families and renders them non–separable. As a consequence every mode shares qualitative properties of both families to some extent, and the resulting mode spectrum consists of modes that bear no direct connection to the normal fluid case. The two–fluid model used here to describe superfluid neutron stars is practically equivalent to those use in previous studies, and we refer the reader to Andersson & Comer (2001a) and Prix et al. (2002) for a more extensive discussion about its physical motivations and justification.

The plan of this paper is as follows: in Sect. 2 we introduce the basic equations for the general two–fluid neutron star model, and in Sect. 3 we develop its linear perturbation equations and show how to recover the special case of a single perfect fluid. In Sect. 4 we specialise to the simpler case of adiabatic oscillations of free, cold fluids, and we derive the necessary boundary conditions. In this section we also show that the separation into two distinct mode families is possible only for non–stratified models. Sect. 5 presents the numerical results concerning the background models, the eigenmode spectrum and its dependence on entrainment (resulting in avoided crossings), as well as the one–fluid results where we recover the expected composition g–modes. We present our conclusions and a discussion of necessary future work in Sect. 6.

## 2 The general two–fluid neutron star model

### 2.1 The general two–fluid equations

The equations and notation for the Newtonian two–fluid neutron star model used here are based on a more general formalism described in Prix (2002), which is the Newtonian analogue of a generally relativistic framework developed by Carter (1989). In this section we will briefly summarise the general model and equations relevant for the present work, and we refer the reader to Prix (2002) for the derivation and more detailed discussion of this model.

We describe a neutron star as a system consisting of two fluids: a superfluid of neutrons, and a normal fluid of comoving constituents, which include protons, electrons and entropy (and generally further particles like muons etc). We denote the particle number densities for neutrons, protons and electrons as $`n_\mathrm{n}`$, $`n_\mathrm{p}`$ and $`n_\mathrm{e}`$ respectively, and we use $`s`$ for the entropy density. The velocities of the two fluids are $`v_\mathrm{n}`$ for the neutron fluid, and $`v_\mathrm{c}`$ for the fluid of comoving constituents, the relative velocity $`\mathrm{\Delta }`$ between the two fluids is therefore

$$\mathrm{\Delta }v_\mathrm{c}v_\mathrm{n}.$$
(1)
On the local (“mesoscopic”) level, a superfluid is constrained to be in a state of irrotational flow, and its angular momentum will be accommodated by a lattice of “microscopic” vortices. However, for the global description of a neutron star we are more interested in a “macroscopic” description of the superfluid, obtained by averaging over volume elements containing a large number of vortices, but which are small compared to the dimensions of the neutron star. In this macroscopic description, the superfluid has a continuous vorticity and nearly behaves like an ordinary fluid (apart from small anisotropies due to the vortex-tension, which we will neglect). The absence of (local) viscosity still allows the superfluid to move relative to the normal fluid, but the presence of the vortex lattice now allows for a direct force between the two fluids. In the case of this force being zero, the vortex lattice moves with the superfluid and we refer to this situation as free fluids. If the vortices are “locked” to the normal fluid (e.g. as can happen in the crust), the mutual force will be non-zero but strictly non–dissipative. This is known as “vortex pinning”. Only in the intermediate cases, where a friction force causes the vortices to have a different velocity from both the normal fluid and the superfluid, energy is dissipated. In this case the mutual force is usually referred to as “mutual friction”.

An essential simplification of the present two–fluid model is that we neglect all electrodynamic effects, as we assume the charge densities of protons and electrons to be strictly balanced, i.e.

$$n_\mathrm{e}=n_\mathrm{p}.$$
(2)
Therefore we can effectively eliminate the electrons from our description, as their density and velocity is entirely specified by the protons.

We note that in “transfusive” models (as first set up in Langlois et al. (1998)), i.e. models which allow for $`\beta `$–reactions ($`\mathrm{n}\mathrm{p}+\mathrm{e}+\overline{\nu }`$) between the two fluids, the total mass in the reaction has to be conserved in a consistent Newtonian description. Therefore we set

$$m_\mathrm{b}m_\mathrm{n}m_\mathrm{p}+m_\mathrm{e}.$$
(3)
The respective mass densities of the two fluids can now be written as

$$\rho _\mathrm{n}=m_\mathrm{b}n_\mathrm{n},\text{and}\rho _\mathrm{c}=m_\mathrm{b}n_\mathrm{p},$$
(4)
and the total mass density $`\rho `$ is simply $`\rho =\rho _\mathrm{n}+\rho _\mathrm{c}`$.

The (local) *kinematics* of the system is completely described (up to arbitrary rotations and boosts) in terms of $`\rho _\mathrm{n}`$, $`\rho _\mathrm{c}`$, $`s`$ and $`\mathrm{\Delta }^2`$. The *dynamics* is determined by the internal energy density function $``$ or equation of state, which is a function of the form $`=(s,\rho _\mathrm{n},\rho _\mathrm{c},\mathrm{\Delta }^2)`$. This energy function defines the first law of thermodynamics for this system by its total differential, namely

$$d=Tds+\stackrel{~}{\mu }^\mathrm{c}d\rho _\mathrm{c}+\stackrel{~}{\mu }^\mathrm{n}d\rho _\mathrm{n}+\alpha d\mathrm{\Delta }^2.$$
(5)
This differential defines the dynamic quantities, namely the temperature $`T`$, the specific chemical potentials $`\stackrel{~}{\mu }^\mathrm{c}`$, $`\stackrel{~}{\mu }^\mathrm{n}`$, and the entrainment $`\alpha `$ as the thermodynamic conjugates of the kinematic quantities $`s`$, $`\rho _\mathrm{c}`$, $`\rho _\mathrm{n}`$ and $`\mathrm{\Delta }`$. The specific chemical potentials $`\stackrel{~}{\mu }`$ are related to the more common definition of the chemical potentials $`\mu `$ via $`\stackrel{~}{\mu }\mu /m_\mathrm{b}`$. We note that $`\mu ^\mathrm{c}`$ is not simply the proton chemical potential, because adding a proton in this model implies adding an electron as well (due to (2)), and therefore one can see that $`\mu ^\mathrm{c}=\mu ^\mathrm{p}+\mu ^\mathrm{e}`$, where $`\mu ^\mathrm{p}`$ and $`\mu ^\mathrm{e}`$ are the respective proton and electron chemical potentials.

The function $`\alpha `$ defined in (5) reflects the dependence of the internal energy on the relative velocity $`\mathrm{\Delta }`$ between the two fluids, which characterises the so–called entrainment effect. This entrainment function $`\alpha `$ has dimensions of a mass density, and it will be useful in the following to define the two dimensionless entrainment functions $`\epsilon _\mathrm{n}`$ and $`\epsilon _\mathrm{c}`$ as

$$\epsilon _\mathrm{c}\frac{2\alpha }{\rho _\mathrm{c}},\text{and}\epsilon _\mathrm{n}\frac{2\alpha }{\rho _\mathrm{n}}=\frac{\rho _\mathrm{c}}{\rho _\mathrm{n}}\epsilon _\mathrm{c}.$$
(6)
The “generalised” pressure $`P`$ is introduced as the Legendre-conjugate of the energy density, namely by the usual thermodynamic relation

$$+P=\rho _\mathrm{n}\stackrel{~}{\mu }^\mathrm{n}+\rho _\mathrm{c}\stackrel{~}{\mu }^\mathrm{c}+sT,$$
(7)
which results in the total differential of the pressure function $`P(\stackrel{~}{\mu }^\mathrm{c},\stackrel{~}{\mu }^\mathrm{n},T,\mathrm{\Delta }^2)`$:

$$dP=\rho _\mathrm{n}d\stackrel{~}{\mu }^\mathrm{n}+\rho _\mathrm{n}d\stackrel{~}{\mu }^\mathrm{c}+sdT\alpha d\mathrm{\Delta }^2.$$
(8)
This generalised pressure $`P`$ can be seen to reduce to the usual definition of the pressure of a perfect fluid in the case of $`\mathrm{\Delta }=0`$.

The equations of motion of this two–fluid system are derived from a “convective” variational principle in Prix (2002), and here we only present the resulting equations. The conservation of energy results in the following equation:

$$T\Gamma _\mathrm{s}=f\mathrm{\Delta }+\beta \Gamma _\mathrm{n},$$
(9)
where $`\Gamma _\mathrm{s}`$ and $`\Gamma _\mathrm{n}`$ are the creation rates of entropy and neutrons respectively, i.e.

$$\Gamma _\mathrm{s}_ts+(sv_\mathrm{c}),\text{and}\Gamma _\mathrm{n}_t\rho _\mathrm{n}+(\rho _\mathrm{n}v_\mathrm{n}),$$
(10)
while the proton creation rate $`\Gamma _\mathrm{c}`$ has to satisfy $`\Gamma _\mathrm{c}=\Gamma _\mathrm{n}`$ for baryon conservation. The quantity $`\beta `$ in (9) characterises the deviation from chemical equilibrium and its explicit expression if found as

$$\beta \stackrel{~}{\mu }^\mathrm{c}\stackrel{~}{\mu }^\mathrm{n}\frac{1}{2}\left(12\epsilon _\mathrm{n}\right)\mathrm{\Delta }^2.$$
(11)
This is the Newtonian analogue of a result first found in the relativistic transfusive model by Langlois et al. (1998). We note that there is an additional kinetic term with respect to the naive $`\stackrel{~}{\mu }^\mathrm{c}\stackrel{~}{\mu }^\mathrm{n}`$ in the case of relative motion. The mutual force density $`f`$ (the sign convention is such that this force acts *on* the neutron fluid) is a direct interaction force between the two fluids, which in the case of a superfluid stems from vortex interactions like pinning or mutual friction. In order to ensure explicitly that the second law of thermodynamics, i.e. $`\Gamma _\mathrm{s}0`$, is satisfied by (9), we can write the neutron creation rate $`\Gamma _\mathrm{n}`$ and the mutual force $`f`$ in the form

$`\Gamma _\mathrm{n}`$ $`=`$ $`\mathrm{\Xi }\beta ,\text{with}\mathrm{\Xi }0,`$ (12)
$`f`$ $`=`$ $`\eta \mathrm{\Delta }+\kappa \times \mathrm{\Delta },\text{with}\eta 0,`$ (13)
where the non–negative functions $`\mathrm{\Xi }`$ and $`\eta `$ govern the beta–reaction rate and the friction force, while the vector $`\kappa `$ allows for a non–dissipative Magnus–type force (i.e. orthogonal to the relative motion). A non–transfusive model, i.e. one that does not allow for beta reactions $`\mathrm{n}\mathrm{c}`$, has $`\mathrm{\Xi }=0`$, and free vortices correspond to $`\eta =0`$ and $`\kappa =0`$, while pinned vortices correspond to $`\eta =0`$ and $`\kappa 0`$.

The momentum equation for the superfluid neutrons is given by

$$(_t+v_\mathrm{n})(v_\mathrm{n}+\epsilon _\mathrm{n}\mathrm{\Delta })+(\stackrel{~}{\mu }^\mathrm{n}+\mathrm{\Phi })+\epsilon _\mathrm{n}\mathrm{\Delta }_av_\mathrm{n}^a=\frac{1}{\rho _\mathrm{n}}f,$$
(14)
while the equation for the normal fluid reads as

$`(_t+v_\mathrm{c})(v_\mathrm{c}\epsilon _\mathrm{c}\mathrm{\Delta })+(\stackrel{~}{\mu }^\mathrm{c}+\mathrm{\Phi })\epsilon _\mathrm{c}\mathrm{\Delta }_av_\mathrm{c}^a`$
$`+{\displaystyle \frac{s}{\rho _\mathrm{c}}}T={\displaystyle \frac{1}{\rho _\mathrm{c}}}f+\left(1\epsilon _\mathrm{c}\epsilon _\mathrm{n}\right){\displaystyle \frac{\Gamma _\mathrm{n}}{\rho _\mathrm{c}}}\mathrm{\Delta }.`$ (15)
The gravitational potential $`\mathrm{\Phi }`$ is related to the mass densities $`\rho _\mathrm{n}`$ and $`\rho _\mathrm{c}`$ via the Poisson equation

$$^2\mathrm{\Phi }=4\pi G(\rho _\mathrm{n}+\rho _\mathrm{c}).$$
(16)
### 2.2 The static equilibrium background

We consider a static background star, so we set

$$v_\mathrm{n}=v_\mathrm{p}=\mathrm{\Delta }=0,\text{and}\Gamma _\mathrm{n}=\Gamma _\mathrm{c}=\Gamma _\mathrm{s}=0,$$
(17)
which by (13) implies the vanishing of the mutual force, i.e. $`f=0`$. Because chemical reactions would not be negligible on long timescales, we also assume the background star to be in chemical equilibrium, i.e.

$$\beta =\stackrel{~}{\mu }^\mathrm{c}\stackrel{~}{\mu }^\mathrm{n}=0.$$
(18)
These equilibrium conditions reduce the equation of motion (14) to

$$\stackrel{~}{\mu }^\mathrm{n}=\stackrel{~}{\mu }^\mathrm{c}=\mathrm{\Phi },$$
(19)
and with (15) this also implies that the background star is in thermal equilibrium, i.e. $`T=0`$.

The static background has to be spherically symmetric, and therefore (19) and (16) lead to the following equation for the background:

$$\stackrel{~}{\mu }^{\prime \prime }(r)+\frac{2}{r}\stackrel{~}{\mu }^{}(r)=4\pi G\rho (r),$$
(20)
where we have introduced the equilibrium chemical potential $`\stackrel{~}{\mu }\stackrel{~}{\mu }^\mathrm{n}=\stackrel{~}{\mu }^\mathrm{c}`$, and the prime () denotes the radial derivative $`d/dr`$.

The equation of state allows one to relate the equilibrium chemical potential $`\stackrel{~}{\mu }`$ directly to the total mass density $`\rho `$ at constant temperature, and therefore the background is fully determined. The numerical method for solving this equation will be discussed in sect. 5.

In the following it will be convenient to use the radius $`R`$ and central density $`\rho _0`$ of the static background as basis units for length and mass density, so the corresponding “natural unit” for frequencies is $`\sqrt{4\pi G\rho _0}`$. All equations in the following are expressed in these natural units except otherwise stated.

### 2.3 Entrainment and effective masses

For small relative velocities $`\mathrm{\Delta }`$, we can separate the “bulk” equation of state from the entrainment by expanding $`(s,\rho _\mathrm{n},\rho _\mathrm{c};\mathrm{\Delta }^2)`$ in terms of $`\mathrm{\Delta }`$, i.e. by writing

$$=^{(0)}(s,\rho _\mathrm{n},\rho _c)+\alpha ^{(0)}(s,\rho _\mathrm{n},\rho _\mathrm{c})\mathrm{\Delta }^2+𝒪(\mathrm{\Delta }^4),$$
(21)
where the quantities $`^{(0)}`$ and $`\alpha ^{(0)}`$ are evaluated at zero relative velocity $`\mathrm{\Delta }`$. The background equation of state $`^{(0)}`$ is therefore decoupled from the entrainment function $`\alpha ^{(0)}`$, and we can specify these two functions independently. The link between the entrainment function $`\alpha `$ and the equivalent description in terms of effective masses $`m^{}`$ (Andreev & Bashkin 1975) has been discussed in previous work (Prix et al. 2002), and it can be shown that one can express $`\alpha `$ in terms of the proton effective mass $`m_\mathrm{p}^{}`$ (which is generally a function of the densities), in the form

$$2\alpha =\rho _\mathrm{c}\left(1\frac{m_\mathrm{p}^{}}{m_\mathrm{b}}\right).$$
(22)
The dimensionless entrainment functions111 We note that Lindblom & Mendell (2000) and more recently Andersson et al. (2002) have used a slightly different dimensionless function $`ϵ`$ to characterise the entrainment effect. The relation between $`ϵ`$ and $`\epsilon _c`$ is given by $`ϵ=\epsilon _c\rho _\mathrm{c}/(\rho _\mathrm{n}\epsilon _c\rho )`$. $`\epsilon _\mathrm{n}`$ and $`\epsilon _\mathrm{c}`$ can then be expressed according to (6) as

$$\epsilon _\mathrm{c}=1\frac{m_\mathrm{p}^{}}{m_\mathrm{b}},\text{and}\epsilon _\mathrm{n}=\frac{\rho _\mathrm{c}}{\rho _\mathrm{n}}\epsilon _\mathrm{c}.$$
(23)
## 3 Linearised perturbation equations

### 3.1 Oscillations of superfluid neutron stars

We consider small perturbations with respect to the static equilibrium background described in Sect. 2.2. Linearising the equations of motion (14) and (15) yields

$`\rho _\mathrm{n}_t(\delta v_\mathrm{n}+\epsilon _\mathrm{n}\delta \mathrm{\Delta })+\rho _\mathrm{n}(\delta \stackrel{~}{\mu }^\mathrm{n}+\delta \mathrm{\Phi })=\delta f,`$ (24)
$`\rho _\mathrm{c}_t(\delta v_\mathrm{c}\epsilon _\mathrm{c}\delta \mathrm{\Delta })+\rho _\mathrm{c}(\delta \stackrel{~}{\mu }^\mathrm{c}+\delta \mathrm{\Phi })=\delta fs\delta T,`$ (25)
where $`\delta Q`$ denotes the *Eulerian* perturbation of the quantity $`Q`$. The perturbation of the mutual force (13) is given by

$$\delta f=\eta \delta \mathrm{\Delta }+\kappa \times \delta \mathrm{\Delta },$$
(26)
and the linearised energy conservation (9) and (12) together with the condition of baryon conservation lead to

$`_t\delta s+(s\delta v_\mathrm{c})`$ $`=`$ $`0,`$ (27)
$`_t\delta \rho _\mathrm{n}+(\rho _\mathrm{n}\delta v_\mathrm{n})`$ $`=`$ $`\mathrm{\Xi }(\delta \stackrel{~}{\mu }^\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{n}),`$ (28)
$`_t\delta \rho _\mathrm{c}+(\rho _\mathrm{c}\delta v_\mathrm{c})`$ $`=`$ $`\mathrm{\Xi }(\delta \stackrel{~}{\mu }^\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{n}).`$ (29)
The perturbed Poisson equation (16) reads (in natural units) as

$$^2\delta \mathrm{\Phi }=\delta \rho _\mathrm{n}+\delta \rho _\mathrm{c}.$$
(30)
The system is closed by specification of the mutual force functions $`\eta `$ and $`\kappa `$, the transfusion function $`\mathrm{\Xi }`$, and an equation of state which allows to express the dynamical quantities $`\delta T`$, $`\delta \stackrel{~}{\mu }^\mathrm{n}`$ and $`\delta \stackrel{~}{\mu }^\mathrm{c}`$ in terms of the kinematic variables $`\delta s`$, $`\delta \rho _\mathrm{n}`$ and $`\delta \rho _\mathrm{c}`$, thereby reducing the number of unknown perturbation quantities to $`13`$, which corresponds exactly to the number of equations.

### 3.2 The special case of normal–fluid neutron stars

It is interesting to compare the superfluid neutron star case with the normal fluid case, where the two constituents $`n`$ and $`\mathrm{c}`$ are moving together and form a single perfect fluid. This case is obviously just a subclass of the two–fluid case discussed so far, namely subject to the additional constraint $`v_\mathrm{n}=v_\mathrm{c}`$, and therefore222Lindblom & Mendell (1994, 1995) have imposed $`\delta \beta =0`$ to recover the perfect fluid case. However, adiabatic oscillations of a perfect fluid only satisfy this condition in *non–stratified* stars (cf. sect. 5.4), the constraint $`\delta \beta =0`$ is therefore generally not met.

$$\delta v_\mathrm{n}=\delta v_\mathrm{c}=\delta v,\delta \mathrm{\Delta }=0.$$
(31)
By linking the two constituents together, the degrees of freedom have been reduced by three, and instead of the individual momentum equations (24) and (25), now only the sum of momenta can be required to be conserved, i.e.

$$\rho \left(_t\delta v+\delta \mathrm{\Phi }\right)+\rho _\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}+\rho _\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c}+s\delta T=0.$$
(32)
We introduce the notation

$$\rho _\mathrm{c}=x_\mathrm{c}\rho ,\rho _\mathrm{n}=x_\mathrm{n}\rho ,\text{and}s=\stackrel{~}{s}\rho ,$$
(33)
for the proton and neutron fractions $`x_\mathrm{c}`$ and $`x_\mathrm{n}`$, and the specific entropy $`\stackrel{~}{s}`$, which allows us to rewrite the one–fluid equation of motion (32) in the slightly more familiar form

$$_t\delta v+\left(\frac{\delta P}{\rho }+\delta \mathrm{\Phi }\right)\left[\delta \beta x_\mathrm{c}+\delta T\stackrel{~}{s}\right]=0,$$
(34)
where we used the fact that the total pressure differential $`dP`$ of (8) in this perfect fluid case reduces to

$$dP=\rho _\mathrm{n}d\stackrel{~}{\mu }^\mathrm{n}+\rho _\mathrm{n}d\stackrel{~}{\mu }^\mathrm{c}+sdT.$$
(35)
We can now compare (32) to the standard expression for the Euler equation of stellar oscillations of non–rotating stars (Cox 1980; Unno et al. 1989), which is usually written as

$$_t\delta v+\left(\frac{\delta P}{\rho }+\delta \mathrm{\Phi }\right)+𝒞(r)(\delta v)A=0,$$
(36)
where $`𝒞(r)`$ is a function of the background, and $`A`$ is the so–called Schwarzschild discriminant that is responsible for presence of g–modes. By comparing eqs. (32) and (34), we see that $`A`$ will be non–zero (indicating the presence of g–modes) whenever $`\stackrel{~}{s}0`$ or $`x_\mathrm{c}0`$. This reflects the well–known fact that any type of stratification, either in specific entropy $`\stackrel{~}{s}`$ or in the chemical composition $`x_\mathrm{c}`$, will result in g–modes, as pointed out by Reisenegger & Goldreich (1992).

## 4 Adiabatic oscillations of ”free”, cold fluids

### 4.1 Reduction to a 1D eigenvalue problem

In order to close the system of perturbation equations (24)–(30), we need specific models for the mutual force $`f`$, the transfusion $`\mathrm{\Xi }`$, in addition to an equation of state of the form $`(\rho _n,\rho _\mathrm{c},\mathrm{s})`$, all of which are highly dependent on microphysical models and are rather poorly known at the present stage. For this reason, we will postpone the inclusion of these effects to future work, and focus on the case of purely *adiabatic* oscillations (i.e. $`\mathrm{\Xi }=0`$) of free fluids (meaning $`f=0`$). We will further neglect temperature effects (which is generally a very good approximation except for very young neutron stars), so we set $`s=0`$ and $`T=0`$. The resulting simplified system of equations is

$`_t\delta \rho _\mathrm{n}+\left(\rho _\mathrm{n}\delta v_\mathrm{n}\right)`$ $`=`$ $`0,`$ (37)
$`_t\delta \rho _\mathrm{c}+\left(\rho _\mathrm{c}\delta v_\mathrm{c}\right)`$ $`=`$ $`0,`$ (38)
$`_t\left[(1\epsilon _\mathrm{n})\delta v_\mathrm{n}+\epsilon _\mathrm{n}\delta v_\mathrm{c}\right]`$ $`=`$ $`\left(\delta \stackrel{~}{\mu }^\mathrm{n}+\delta \mathrm{\Phi }\right),`$ (39)
$`_t\left[(1\epsilon _\mathrm{c})\delta v_\mathrm{c}+\epsilon _\mathrm{c}\delta v_\mathrm{n}\right]`$ $`=`$ $`\left(\delta \stackrel{~}{\mu }^\mathrm{c}+\delta \mathrm{\Phi }\right),`$ (40)
$`^2\delta \mathrm{\Phi }`$ $`=`$ $`\delta \rho _\mathrm{n}+\delta \rho _\mathrm{c},.`$ (41)
We point out that this system of equations is identical to the one used in Andersson & Comer (2001a), which was obtained from a Newtonian limit of the relativistic equations. It is also related to the equations of Lindblom & Mendell (1994), which are expressed in the “orthodox” formulation of superfluids (Landau & Lifshitz 1959), while the present description is based on the “canonical” approach introduced by Carter (1989).

Using the equation of state we can link the density perturbations $`\delta \rho _\mathrm{X}`$ to $`\delta \stackrel{~}{\mu }^\mathrm{X}`$ (with the constituent index notation $`\mathrm{X}=\mathrm{n},\mathrm{c}`$) to linear order, namely

$$\begin{array}{cc}\hfill \delta \rho _\mathrm{n}=& 𝒮_{\mathrm{nn}}\delta \stackrel{~}{\mu }^\mathrm{n}+𝒮_{\mathrm{nc}}\delta \stackrel{~}{\mu }^\mathrm{c},\hfill \\ \hfill \delta \rho _\mathrm{c}=& 𝒮_{\mathrm{cn}}\delta \stackrel{~}{\mu }^\mathrm{n}+𝒮_{\mathrm{cc}}\delta \stackrel{~}{\mu }^\mathrm{c},\hfill \end{array}$$
(42)
where the symmetric “structure matrix” $`𝒮_{\mathrm{XY}}`$ is defined as

$$\left(𝒮^1\right)_{\mathrm{XY}}\frac{\stackrel{~}{\mu }^\mathrm{X}}{\rho _\mathrm{Y}}=\frac{^2}{\rho _\mathrm{X}\rho _\mathrm{Y}},\mathrm{X},\mathrm{Y}\{\mathrm{n},\mathrm{c}\}.$$
(43)
Due to the dual role of the pressure $`P`$ in (7) and (8), we can equivalently express $`𝒮_{\mathrm{XY}}`$ as

$$𝒮_{\mathrm{XY}}=\frac{\rho _\mathrm{X}}{\stackrel{~}{\mu }^\mathrm{Y}}=\frac{^2P}{\stackrel{~}{\mu }^\mathrm{X}\stackrel{~}{\mu }^\mathrm{Y}}.$$
(44)
We note that although we have assumed free fluids, i.e. there is no direct force acting between them, the fluids are nevertheless locally coupled by the equation of state; we can distinguish two sources of this coupling, one is due to the non–diagonal term $`𝒮_{\mathrm{nc}}`$ in (42), while the second is due to the entrainment terms $`\epsilon _\mathrm{X}`$.

The background quantities $`\stackrel{~}{\mu }^\mathrm{X}`$ can be seen in (19) to behave like the gravitational potential $`\mathrm{\Phi }`$; this means in particular that their gradient is always finite, even at the surface. Therefore $`\delta \stackrel{~}{\mu }^\mathrm{X}`$ is finite everywhere, contrary to $`\delta \rho _\mathrm{X}`$ which can diverge at the surface when $`\rho ^{}\mathrm{}`$. This is seen from the relation $`\mathrm{\Delta }\rho =\delta \rho +\xi ^r\rho ^{}(r)`$ between the Lagrangian perturbation $`\mathrm{\Delta }\rho `$ and the Eulerian $`\delta \rho `$, for a radial displacement $`\xi ^r`$. On physical grounds $`\mathrm{\Delta }\rho `$ must be bounded everywhere (as it reflects the physical property of a fluid element), while the first–order Eulerian quantity $`\delta \rho `$ diverges at the surface whenever $`\rho ^{}\mathrm{}`$ and $`\xi ^r0`$ at $`r=R`$. This might seem problematic for the validity of the equations, but it only reflects the fact that in this case even an infinitesimal displacement of the surface will lead to a *finite* (as opposed to infinitesimal) Eulerian density change there. By considering Lagrangian instead of Eulerian variables, it can be shown that the physical solution is still well behaved even if $`\delta \rho \mathrm{}`$ at the surface. In this case the first–order quantity $`\delta \rho `$ no longer approximates the physical Eulerian density change, but the divergence is such that the *Lagrangian* first–order quantity $`\mathrm{\Delta }\rho `$ is still perfectly regular. If one wanted to impose that $`\delta \rho `$ should be bounded everywhere (as did Lindblom & Mendell (1994)), then this situation would be inverted and the Lagrangian quantity $`\mathrm{\Delta }\rho `$ would diverge, which is unphysical indeed.

From a numerical point of view it seemed better to solve directly for the well-behaved $`\delta \stackrel{~}{\mu }^\mathrm{X}`$ instead of the potentially diverging $`\delta \rho _\mathrm{X}`$, by using (42) to substitute for $`\delta \rho _\mathrm{X}`$. We note that the coefficients $`𝒮_{\mathrm{XY}}`$ in this expression will generally diverge (or vanish) at the surface, depending on the equation of state and reflecting the behaviour of $`\delta \rho _\mathrm{X}`$. The system of equations (37-41) for eigenmode solutions of the form $`\delta Q(x,t)=\delta Q(x)e^{i\omega t}`$ now yields

$`(\rho _\mathrm{n}\delta v_\mathrm{n})`$ $`=`$ $`i\omega \left[𝒮_{\mathrm{nn}}\delta \stackrel{~}{\mu }^\mathrm{n}+𝒮_{\mathrm{nc}}\delta \stackrel{~}{\mu }^\mathrm{c}\right],`$ (45)
$`(\rho _\mathrm{c}\delta v_\mathrm{c})`$ $`=`$ $`i\omega \left[𝒮_{\mathrm{cn}}\delta \stackrel{~}{\mu }^\mathrm{n}+𝒮_{\mathrm{cc}}\delta \stackrel{~}{\mu }^\mathrm{c}\right],`$ (46)
$`(\delta \stackrel{~}{\mu }^\mathrm{n}+\delta \mathrm{\Phi })`$ $`=`$ $`i\omega \left[(1\epsilon _\mathrm{n})\delta v_\mathrm{n}+\epsilon _\mathrm{c}\delta v_\mathrm{c}\right],`$ (47)
$`(\delta \stackrel{~}{\mu }^\mathrm{c}+\delta \mathrm{\Phi })`$ $`=`$ $`i\omega \left[(1\epsilon _\mathrm{c})\delta v_\mathrm{c}+\epsilon _\mathrm{n}\delta v_\mathrm{n}\right],`$ (48)
$`^2\delta \mathrm{\Phi }`$ $`=`$ $`k_\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}+k_\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c},`$ (49)
where we have introduced the convenient “structure vector” $`k_\mathrm{X}`$, which is defined as

$$k_\mathrm{X}\underset{\mathrm{Y}=\mathrm{n},\mathrm{c}}{}S_{\mathrm{XY}},\mathrm{X}\{\mathrm{n},\mathrm{c}\}.$$
(50)
For a spherically symmetric background, we can separate the radial and angular dependence and obtain solutions with definite quantum numbers $`l`$ and $`|m|l`$ using the ansatz

$`\delta \mathrm{\Phi }(r,\theta ,\phi )`$ $`=`$ $`\delta \mathrm{\Phi }(r)Y_l^m(\theta ,\phi ),`$ (51)
$`\delta \stackrel{~}{\mu }^\mathrm{X}(r,\theta ,\phi )`$ $`=`$ $`\delta \stackrel{~}{\mu }^\mathrm{X}(r)Y_l^m(\theta ,\phi ),`$ (52)
$`\delta v_\mathrm{X}(r,\theta ,\phi )`$ $`=`$ $`{\displaystyle \frac{W_\mathrm{X}(r)}{r}}R_l^m+V_\mathrm{X}(r)S_l^miU_\mathrm{X}(r)T_l^m,`$ (53)
where the spherical harmonics $`Y_l^m(\theta ,\phi )`$ are the eigenfunctions of $`r^2^2Y_l^m=l(l+1)Y_l^m`$, and $`R`$, $`S`$ and $`T`$ form the orthogonal harmonic basis, defined as

$$R_l^mY_l^mr,S_l^mY_l^m,T_l^m\times R,$$
(54)
see Rieutord (1987) for details. The three–dimensional eigenvalue problem (45)-(49) has now been reduced to the following one–dimensional problem:

$`\left(r\rho _\mathrm{n}W_\mathrm{n}\right)^{}l(l+1)\rho _\mathrm{n}V_\mathrm{n}`$ $`=`$ $`i\omega r^2\left[𝒮_{\mathrm{nn}}\delta \stackrel{~}{\mu }^\mathrm{n}+𝒮_{\mathrm{nc}}\delta \stackrel{~}{\mu }^\mathrm{c}\right],`$ (55)
$`\left(r\rho _\mathrm{c}W_\mathrm{c}\right)^{}l(l+1)\rho _\mathrm{c}V_\mathrm{c}`$ $`=`$ $`i\omega r^2\left[𝒮_{\mathrm{nc}}\delta \stackrel{~}{\mu }^\mathrm{n}+𝒮_{\mathrm{cc}}\delta \stackrel{~}{\mu }^\mathrm{c}\right],`$ (56)
$`r\left(\delta \stackrel{~}{\mu }_{}^{\mathrm{n}}{}_{}{}^{}+\delta \mathrm{\Phi }^{}\right)`$ $`=`$ $`i\omega [(1\epsilon _\mathrm{n})W_\mathrm{n}+\epsilon _\mathrm{n}W_\mathrm{c})],`$ (57)
$`r\left(\delta \stackrel{~}{\mu }_{}^{\mathrm{c}}{}_{}{}^{}+\delta \mathrm{\Phi }^{}\right)`$ $`=`$ $`i\omega [(1\epsilon _\mathrm{c})W_\mathrm{c}+\epsilon _\mathrm{c}W_\mathrm{n})],`$ (58)
$`\delta \stackrel{~}{\mu }^\mathrm{n}\delta \mathrm{\Phi }`$ $`=`$ $`i\omega [(1\epsilon _\mathrm{n})V_\mathrm{n}+\epsilon _\mathrm{n}V_\mathrm{c})],`$ (59)
$`\delta \stackrel{~}{\mu }^\mathrm{c}\delta \mathrm{\Phi }`$ $`=`$ $`i\omega [(1\epsilon _\mathrm{c})V_\mathrm{c}+\epsilon _\mathrm{c}V_\mathrm{n})],`$ (60)
$`\left(r^2\delta \mathrm{\Phi }^{}\right)^{}l(l+1)\delta \mathrm{\Phi }`$ $`=`$ $`r^2(k_\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}+k_\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c}).`$ (61)
The axial velocity component $`U_\mathrm{X}`$ is decoupled and corresponds to zero frequency $`\omega =0`$, therefore all non–zero frequency eigenmodes are purely polar. We note that the horizontal velocity equations (59) and (60) only hold for $`l>0`$, because $`S_0^0=0`$.

### 4.2 Boundary and regularity conditions

#### 4.2.1 At the center

It can be shown that the representation (51)-(53) of a regular physical quantity requires the following asymptotic behaviour of the radial functions as $`r0`$,

$$\delta \stackrel{~}{\mu }^\mathrm{X}\delta \mathrm{\Phi }W_\mathrm{X}V_\mathrm{X}𝒪(r^l),\mathrm{X}\{\mathrm{n},\mathrm{c}\},$$
(62)
where $`𝒪(r^l)`$ means of order $`r^l`$ *or higher*. Another requirement is the regularity of the solution at singular points of the equations. From eqs. (57) and (58) we see that $`W_\mathrm{X}`$ must vanish as $`𝒪(r)`$ at the origin, while from eqs. (59) and (60) we can derive the further regularity requirement

$$\left(1+r\frac{\rho _\mathrm{X}^{}}{\rho _\mathrm{X}}\right)W_\mathrm{X}+rW_\mathrm{X}^{}l(l+1)V_\mathrm{X}𝒪(r^2).$$
(63)
Regularity of solutions of Poisson’s equation (61) requires

$$\delta \mathrm{\Phi }𝒪(r^2)\text{for}l>0,$$
(64)
while for radial oscillations ($`l=0`$) only $`\delta \mathrm{\Phi }^{}𝒪(r)`$ is required. These constraints are automatically satisfied by (62) for $`l2`$, but they are stronger requirements than (62) in the cases $`l=1`$ and $`l=0`$.

#### 4.2.2 At the surface

At the outer surface ($`r=R`$) we need to ensure the continuity of the gravitational potential $`\mathrm{\Phi }`$, which results (e.g. see Ledoux & Walraven (1958)) in the boundary condition

$$\delta \mathrm{\Phi }^{}(R)+\frac{l+1}{R}\delta \mathrm{\Phi }(R)=\xi ^r\rho (R),$$
(65)
where $`\xi ^r`$ is the radial displacement of the surface, and $`\rho (R)`$ is the *inner* limit of $`\rho (r)`$, i.e. $`\rho (R)=lim_{rR}\rho (r)`$. In the present work we will only consider stars with vanishing density at the surface333We note that Lindblom & Mendell (1994, 1995) have set the right–hand side of (65) to zero in their homogeneous model, which is inconsistent if one allows for surface displacements $`\xi ^r`$. i.e. $`\rho (r)0`$. The conservation equations (55) and (56) contain a (regular) singularity at the surface $`r=R`$, because $`\rho ^{}/\rho `$ diverges as $`\rho 0`$. Therefore we need an additional regularity condition that can be obtained by a Frobenius–type expansion around the surface444We can conveniently use $`\stackrel{~}{\mu }`$ as an effective radial coordinate at the surface, because it is well-behaved and monotonic (i.e. $`\stackrel{~}{\mu }(R)=0`$ and $`\stackrel{~}{\mu }^{}(R)<0`$).,i.e. $`\stackrel{~}{\mu }=0`$. We first rewrite these equations as

$$(rW_\mathrm{X})^{}l(l+1)V_\mathrm{X}=\left[r\frac{\rho _\mathrm{X}^{}}{\rho _\mathrm{X}}W_\mathrm{X}+i\omega r^2\frac{𝒮_{\mathrm{XY}}\delta \stackrel{~}{\mu }^\mathrm{Y}}{\rho _\mathrm{X}}\right],$$
(66)
and using the definition (43) of the matrix $`𝒮_{\mathrm{XY}}`$ together with the equilibrium condition $`\stackrel{~}{\mu }^\mathrm{n}{}_{}{}^{}=\stackrel{~}{\mu }^\mathrm{c}{}_{}{}^{}=\stackrel{~}{\mu }^{}`$, we can write

$$\rho _\mathrm{X}^{}=k_\mathrm{X}\stackrel{~}{\mu }^{},$$
(67)
where $`k_\mathrm{X}`$ has been defined in (50). In order to simplify the analysis, we will assume that the equation of state becomes decoupled in the limit of vanishing density, i.e.

$$P(\stackrel{~}{\mu }^\mathrm{n},\stackrel{~}{\mu }^\mathrm{c})P_\mathrm{n}(\stackrel{~}{\mu }^\mathrm{n})+P_\mathrm{c}(\stackrel{~}{\mu }^\mathrm{c}),$$
(68)
such that $`\rho _\mathrm{n}\rho _\mathrm{n}(\stackrel{~}{\mu }^\mathrm{n})`$ and $`\rho _\mathrm{c}\rho _\mathrm{c}(\stackrel{~}{\mu }^\mathrm{c})`$, and $`k_\mathrm{X}𝒮_{\mathrm{XX}}`$ in this limit. In this case we can write the asymptotic form of (66) as

$$(rW_\mathrm{X})^{}l(l+1)V_\mathrm{X}=\frac{𝒮_{\mathrm{XX}}}{\rho _\mathrm{X}}\left[r\stackrel{~}{\mu }^{}W_\mathrm{X}+i\omega r^2\delta \stackrel{~}{\mu }^\mathrm{X}\right],$$
(69)
where the expression in brackets is regular for $`\stackrel{~}{\mu }0`$. The behaviour of $`𝒮_{\mathrm{XX}}/\rho _\mathrm{X}`$ in this limit can be analysed by writing $`\rho _\mathrm{X}`$ in the asymptotic form $`\rho _\mathrm{X}(\stackrel{~}{\mu }^\mathrm{X})^\alpha `$ (with positive $`\alpha `$), which is possible due to the assumed decoupling. Therefore we find $`𝒮_{\mathrm{XX}}=\rho _\mathrm{X}/\stackrel{~}{\mu }^\mathrm{X}\stackrel{~}{\mu }^{\alpha 1}`$, and so we have identified the singularity in (69) as

$$\frac{𝒮_{\mathrm{XX}}}{\rho _\mathrm{X}}\stackrel{~}{\mu }^1.$$
(70)
Regularity of the solution therefore requires the following asymptotic behaviour at the surface:

$$r\stackrel{~}{\mu }^{}W_\mathrm{X}+i\omega r^2\delta \stackrel{~}{\mu }^\mathrm{X}𝒪(\stackrel{~}{\mu }),\text{for}\mathrm{X}\{\mathrm{n},\mathrm{c}\}.$$
(71)
### 4.3 Decoupling “ordinary” and “superfluid” modes?

In this section we discuss a change of variables that has been used in several previous studies of oscillations of superfluid neutron stars (Lindblom & Mendell 1994, 1995, 2000; Sedrakian & Wasserman 2000; Andersson & Comer 2001a), namely

$$\begin{array}{cccc}\hfill \delta \mathrm{\Delta }& \delta v_\mathrm{c}\delta v_\mathrm{n},\hfill & \hfill \delta \beta & \delta \stackrel{~}{\mu }^\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{n},\hfill \\ \hfill \delta v& x_\mathrm{c}\delta v_\mathrm{c}+x_\mathrm{n}\delta v_\mathrm{n},\hfill & \hfill \delta \stackrel{~}{\mu }& x_\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c}+x_\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}.\hfill \end{array}$$
(72)
This choice of variables is motivated by the intuitive idea that the additional degrees of freedom of a second fluid should allow for two different types of motion, characterised roughly by the two fluids being either “co-moving” or “counter-moving”, and which are sometimes referred to as “ordinary” and “superfluid” modes. By choosing such “adapted” coordinates (72), one might hope to separate, or at least simplify the system of equations, but we will see that this is generally not the case.

Using the relations (42) and the definition (67) of $`k_\mathrm{n}`$ and $`k_\mathrm{c}`$, and defining $`kk_\mathrm{n}+k_\mathrm{c}`$, we can express the density perturbations $`\delta \rho _\mathrm{X}`$ as

$$\begin{array}{cc}\hfill \delta \rho _\mathrm{n}=& k_\mathrm{n}\delta \stackrel{~}{\mu }+\left(𝒮_{\mathrm{nc}}x_\mathrm{c}k_\mathrm{n}\right)\delta \beta ,\hfill \\ \hfill \delta \rho _\mathrm{c}=& k_\mathrm{c}\delta \stackrel{~}{\mu }+\left(𝒮_{\mathrm{nc}}+x_\mathrm{n}k_\mathrm{c}\right)\delta \beta ,\hfill \end{array}$$
(73)
and therefore

$$\delta \rho =k\delta \stackrel{~}{\mu }(x_\mathrm{c}k_\mathrm{n}x_\mathrm{n}k_\mathrm{c})\delta \beta .$$
(74)
In the case of a spherically symmetric background considered here, we can use (67) to express

$$x_\mathrm{c}k_\mathrm{n}x_\mathrm{n}k_\mathrm{c}=k\frac{\rho }{\rho ^{}}x_\mathrm{c}^{},$$
(75)
and using these variables and relations, we now rewrite the eigenmode equations (37)–(41) in the form

$`k_t\delta \stackrel{~}{\mu }+(\rho \delta v)`$ $`=`$ $`x_\mathrm{c}^{}\left[k{\displaystyle \frac{\rho }{\rho ^{}}}_t\delta \beta \right],`$ (76)
$`_t\delta v+(\delta \stackrel{~}{\mu }+\delta \mathrm{\Phi })`$ $`=`$ $`x_\mathrm{c}^{}\left[\delta \beta e^r\right],`$ (77)
$`^2\delta \mathrm{\Phi }k\delta \stackrel{~}{\mu }`$ $`=`$ $`x_\mathrm{c}^{}\left[k{\displaystyle \frac{\rho }{\rho ^{}}}\delta \beta \right],`$ (78)
$`𝒦_t\delta \beta +\left(\rho x_\mathrm{n}x_\mathrm{c}\delta \mathrm{\Delta }\right)`$ $`=`$ $`x_\mathrm{c}^{}[k{\displaystyle \frac{\rho }{\rho ^{}}}_t\delta \stackrel{~}{\mu }+\rho \delta v^r],`$ (79)
$`\left(1{\displaystyle \frac{2\alpha \rho }{\rho _\mathrm{n}\rho _c}}\right)_t\delta \mathrm{\Delta }+\delta \beta `$ $`=`$ $`\mathrm{\hspace{0.33em}\hspace{0.33em}0},`$ (80)
where $`e^r`$ is the radial basis vector, and we introduced the abbreviation $`𝒦x_\mathrm{c}^2k_\mathrm{n}+x_\mathrm{n}^2k_\mathrm{c}𝒮_{\mathrm{nc}}`$.

We see that the “ordinary”-type of motion ($`\delta \stackrel{~}{\mu }`$, $`\delta v`$, $`\delta \mathrm{\Phi }`$) does not decouple from the “superfluid”-type variables ($`\delta \mathrm{\Delta }`$, $`\delta \beta `$) whenever there is *stratification*, i.e. when $`x_\mathrm{c}^{}0`$ ! This can be understood as follows: while a non–zero relative velocity $`\delta \mathrm{\Delta }`$ can be regarded as a characteristic of a superfluid mode (as opposed to modes in a single fluid), the chemical equilibrium deviation $`\delta \beta `$ is generally non–zero even for a single (but non–barotropic) fluid. In a stratified fluid, any general adiabatic motion will drive a fluid element out of equilibrium, i.e. $`\delta \beta `$ nonzero is not characteristic for either “superfluid” or “ordinary” modes (contrary to claims in Lindblom & Mendell (1994, 1995)), it is a general feature of modes in stratified fluids, and therefore the choice of variables (72) does not lead to a decoupling of the system in this case.

However, it is interesting to consider for a moment this special case of a non–stratified background (which probably never applies in real neutron stars). Setting the proton fraction $`x_\mathrm{c}`$ to a constant, we can separate the equations into two decoupled sets. One system describes “ordinary” modes, namely

$`k_t\delta \stackrel{~}{\mu }+(\rho \delta v)`$ $`=`$ $`0,`$ (81)
$`_t\delta v+(\delta \stackrel{~}{\mu }+\delta \mathrm{\Phi })`$ $`=`$ $`0,`$ (82)
$`^2\delta \mathrm{\Phi }k\delta \stackrel{~}{\mu }`$ $`=`$ $`0,`$ (83)
which are seen to be independent of the entrainment $`\epsilon `$ as well as of the coupling through the “bulk” equation of state, i.e. $`𝒮_{\mathrm{nc}}`$. The second system of equations governs the “superfluid” modes and reads as

$`\left(k{\displaystyle \frac{𝒮_{\mathrm{nc}}}{x_\mathrm{n}x_\mathrm{c}}}\right)_t\delta \beta +(\rho \delta \mathrm{\Delta })`$ $`=`$ $`0.`$ (84)
$`\left(1{\displaystyle \frac{2\alpha \rho }{\rho _\mathrm{n}\rho _\mathrm{c}}}\right)_t\delta \mathrm{\Delta }+\delta \beta `$ $`=`$ $`0,`$ (85)
We see that contrary to the ordinary modes, the superfluid modes do depend on the coupling through entrainment $`\alpha `$ and the equation of state, i.e. $`𝒮_{\mathrm{nc}}`$, but they are completely decoupled from the gravitational perturbation $`\delta \mathrm{\Phi }`$, as they leave the total density unchanged, i.e. $`\delta \rho =0`$.

We have therefore shown that in the non–stratified case there exist two separate families, namely “ordinary” modes $`(\delta v,\delta \stackrel{~}{\mu },\delta \mathrm{\Phi },0,0)`$ and “superfluid” modes $`(0,0,0,\delta \mathrm{\Delta },\delta \beta )`$. One of our numerical models (see next section) has a constant proton fraction $`x_\mathrm{c}`$, and we will see the present analysis confirmed by the numerical results for this model. In the general stratified case, however, these two mode families are coupled and such a clearcut separation is not possible.

## 5 Numerical results

### 5.1 Equation of state: two–constituent polytropes

We use a simple class of two–constituent equations of state which is very convenient to explore the properties of a two–fluid system, namely the following “generalised polytrope”, defined as555This equation of state has been used previously by Comer et al. (1999) and Andersson et al. (2002) to study two–fluid oscillations in general relativity.

$$^{(0)}(\rho _\mathrm{n},\rho _\mathrm{c})=\kappa _\mathrm{n}\rho _\mathrm{n}^{\gamma _\mathrm{n}}+\kappa _\mathrm{c}\rho _\mathrm{c}^{\gamma _\mathrm{c}}$$
(86)
which simply consists of the sum of two ordinary polytropes. For regularity of the chemical potentials (5) in the limit $`\rho _\mathrm{n}0`$ and $`\rho _\mathrm{c}0`$, the polytropic indices must satisfy $`\gamma _\mathrm{n}1`$ and $`\gamma _\mathrm{c}1`$. This equation of state allows the explicit inversion

$$\rho _\mathrm{n}(\stackrel{~}{\mu }^\mathrm{n})=\left(\frac{\stackrel{~}{\mu }^\mathrm{n}}{\kappa _\mathrm{n}\gamma _\mathrm{n}}\right)^{N_\mathrm{n}},\text{and}\rho _\mathrm{c}(\stackrel{~}{\mu }^\mathrm{c})=\left(\frac{\stackrel{~}{\mu }^\mathrm{c}}{\kappa _\mathrm{c}\gamma _\mathrm{c}}\right)^{N_\mathrm{c}},$$
(87)
where we have introduced $`N_\mathrm{X}1/(\gamma _\mathrm{X}1)`$. We see that in chemical equilibrium, i.e. $`\stackrel{~}{\mu }^\mathrm{n}=\stackrel{~}{\mu }^\mathrm{c}`$, the two fluids share a common outer surface. The equilibrium proton fraction $`x_\mathrm{c}`$ can be expressed as

$$x_\mathrm{c}=\left(1+\frac{(\gamma _\mathrm{c}\kappa _\mathrm{c})^{N_\mathrm{c}}}{(\gamma _\mathrm{n}\kappa _\mathrm{n})^{N_\mathrm{n}}}\stackrel{~}{\mu }^{N_\mathrm{n}N_\mathrm{c}}\right)^1,$$
(88)
which shows that the proton fraction is constant whenever $`\gamma _\mathrm{n}=\gamma _\mathrm{c}`$, while the behaviour in the case of different indices falls into the two categories:

$`0\stackrel{\stackrel{~}{\mu }\mathrm{}}{}x_\mathrm{c}\stackrel{\stackrel{~}{\mu }0}{}1,\text{for}\gamma _\mathrm{n}<\gamma _\mathrm{c},`$ (89)
$`1\stackrel{\stackrel{~}{\mu }\mathrm{}}{}x_\mathrm{c}\stackrel{\stackrel{~}{\mu }0}{}0,\text{for}\gamma _\mathrm{n}>\gamma _\mathrm{c}.`$ (90)
### 5.2 Calculating the background models

The equilibrium background solution is determined by equation (20), together with the regularity requirement $`\stackrel{~}{\mu }^{}(0)=0`$, and the boundary condition of vanishing pressure at the surface, i.e.666Strictly speaking $`\mu `$ is determined only up to a constant, which is usually fixed such that $`\stackrel{~}{\mu }`$ vanishes together with the pressure. $`\stackrel{~}{\mu }(1)=0`$, where the surface of the static background star is situated at $`R=1`$ in the natural units defined in Sect. 2.2. Using the equation of state we can express $`\rho =\rho (\stackrel{~}{\mu })`$ in chemical equilibrium, and therefore equation (20) can be written as the following nonlinear eigenvalue problem,

$$\stackrel{~}{\mu }^{\prime \prime }+\frac{2}{r}\stackrel{~}{\mu }=\lambda \rho (\stackrel{~}{\mu }),$$
(91)
where the (dimensionless) eigenvalue $`\lambda `$ is given by

$$\lambda =4\pi G\frac{\rho _0}{\stackrel{~}{\mu }(\rho _0)}R^2.$$
(92)
The eigenvalue $`\lambda `$ represents the a–priori unknown radius of the star, and thereby the actual values of our natural units for a given central density $`\rho _0`$. The method used here to solve this equation is to iterate a linear eigenvalue problem that converges to the solution of (91). This can be done by solving in step $`k+1`$ the following linear eigenvalue problem,

$$\stackrel{~}{\mu }_{k+1}^{\prime \prime }+\frac{2}{r}\stackrel{~}{\mu }_{k+1}^{}=\lambda _{k+1}\left(\frac{\rho (\stackrel{~}{\mu }_k)}{\stackrel{~}{\mu }_k}\right)\stackrel{~}{\mu }_{k+1},$$
(93)
where $`\stackrel{~}{\mu }_k`$ is the solution of the of the previous step $`k`$. The equation in each step is solved using the spectral linear eigenvalue solver package LSB developed by L. Valdettaro and M. Rieutord. With a resolution of 40 Chebychev polynomials and a Gauss–Lobatto collocation method, this iteration converges to about machine precision (i.e. $`10^{15}`$ relative difference between successive steps) in about $`20`$ steps. Another practical advantage of this method is that we calculate the background quantities on the same Gauss–Lobatto grid we use for the numerical code for the eigenmodes. Even when using a different resolution for the eigenmode–calculation, the Chebychev expansion provides a canonical $`𝒞^{\mathrm{}}`$ interpolation which allows us to easily “re-grid” the background solution.

In principle we can solve the background for any given equation of state, but we restrict our attention in this work to the class of two–constituent polytropes (86).The reason for this choice is their analytic simplicity and because our main emphasis is to clarify the qualitative properties of superfluid neutron stars rather than to construct a physically very precise model. This would in any case be quite impossible in a Newtonian description because of the neglect of relativistic effects, in addition to the important uncertainties in our current knowledge of the equation of state of neutron stars.

In the numerical analysis we consider the three different background models defined in table 1, and which are represented in Fig. 1. These three models correspond to three different types of behaviour at the surface. In the case of Model I one can easily find the background solution analytically, as shown in Prix et al. (2002), which allows us to check the numerical method of calculating the background, and we find a maximal relative error of $`10^{14}`$ between the numerical and the analytic solution for model I. Model II represents a generic “stiff” model similar to those used in Comer et al. (1999) and Andersson et al. (2002), which has infinite density gradients at the surface. Model III is of a “soft” type with vanishing density gradients at the surface. These different types of behaviour at the surface are quite analogous to the case of the usual one–constituent polytropes for different polytropic indices.

We note that model I is the only non–stratified model, because it has a constant proton fraction $`x_\mathrm{c}`$ (and finite density gradients at the surface), while models II and III have a non–zero composition gradient, i.e. $`x_\mathrm{c}^{}0`$, as expected from (88).

For easier comparison of the frequencies given in the next section in units of $`\sqrt{4\pi G\rho _0}`$, we provide in table 2 the conversion factors into three important systems of units, namely the SI unit Hz, the “Cox” units $`\sqrt{GM/R^3}`$ (variants thereof, like those used by Lindblom & Mendell (1994), only differ by a constant factor) and the “geometric” units $`c^3/GM`$ typically used in general relativity (Comer et al. 1999; Andersson et al. 2002).

### 5.3 The two–fluid oscillation modes

The eigenmode equations (57)-(61) together with the boundary conditions of Sect. 4.2 form a linear eigenvalue system which we solve numerically using the spectral solver of the LSB–package. The convergence of the results was determined by increasing the resolution starting from $`40`$ Chebychev polynomials up to $`80`$, and we found the changes in frequency decrease very quickly to about $`10^9`$ (or better), which is why we give the frequencies with nine decimals in tables 3 to 5, corresponding roughly to the convergence achieved by the numerical method. This is for future reference and comparison, not because these frequencies represent a physically measurable prediction in any sense.

#### 5.3.1 Eigenmodes of “locally uncoupled” fluids

In this section we consider the case of zero entrainment, i.e. $`\alpha =0`$ and $`\epsilon _\mathrm{n}=\epsilon _\mathrm{c}=0`$. We refer to this situation as “locally uncoupled” fluids, as it is important to note that the two fluids are nevertheless coupled “globally” through the perturbation of the gravitational potential $`\delta \mathrm{\Phi }`$ and (61). We consider the cases of radial ($`l=0`$), dipolar ($`l=1`$) and quadrupolar ($`l=2`$) oscillations, which differ qualitatively in some properties and boundary conditions (see Sect. 4.2), while all higher $`l`$ cases are qualitatively very similar to $`l=2`$. The lowest eigen-frequencies for these three values of $`l`$ are shown in tables 3, 4 and 5 respectively. We label these modes in analogy to the one–fluid case as f- and p- modes, and group them in pairs where the lower frequency mode is labelled as “$`\mathrm{o}`$” and the higher frequency one as “$`\mathrm{s}`$”. The pairs of p–modes are indexed in the order of increasing frequency. We emphasise that this labelling is a pure convention, as one can generally *not* say that these modes would be either co– or counter-moving, or that the subscript would exactly reflect the number of radial nodes.

Let us first consider the special separable case of the non–stratified model I. The first three pairs of eigenfunctions are presented in Fig. 2, and we see that in the “$`\mathrm{o}`$” modes the two fluids are comoving, resulting in a non–zero $`\delta \mathrm{\Phi }`$, and they also remain in strict chemical equilibrium, i.e. $`\delta \beta =0`$. These “ordinary” modes are actually identical to the normal–fluid modes of the same background (see Sect. 5.4). In the case of the $`\mathrm{s}`$–type modes the two fluids are counter–moving in exactly such a way that the total density remains constant, i.e. $`\delta \rho =0`$ and therefore $`\delta \mathrm{\Phi }=0`$, while the fluids are driven out of chemical equilibrium, i.e. $`\delta \beta 0`$. The number of radial nodes in $`\delta v^r`$ is the same for the $`\mathrm{o}`$ and $`\mathrm{s}`$ modes, and corresponds exactly to their index. All these results confirm the analytic predictions for non–stratified models in Sect. 4.3.

However, it would be wrong to assume that these properties are generally true for superfluid oscillations. Stratification makes this picture more complex, even in the case of locally uncoupled fluids considered in this section. If we look at the first three pairs of eigenfunctions for the model II in Fig. 3, we see that the “$`\mathrm{o}`$” modes are not comoving at all (only the fo is nearly comoving), and they have non–zero $`\delta \beta `$, while the “$`\mathrm{s}`$” modes have non–zero $`\delta \rho `$ and $`\delta \mathrm{\Phi }`$. One can not say either that the relative amplitude of $`\delta \beta `$ would be different between the $`\mathrm{o}`$– and $`\mathrm{s}`$– cases, as Lindblom & Mendell (1994) wrongly induced from the properties of the fo–mode, the only eigenmode they presented.

In the case of the low–order modes presented in Fig. 2, the number of radial nodes still seems to correspond to the index, and the fluids are roughly in opposite phase in the $`\mathrm{s}`$– modes, while they are approximately in phase for the $`\mathrm{o}`$–modes. Even this, however, is not true in general, as can be seen in Fig. 4 which shows the higher order p6–modes. In this case neither $`\mathrm{o}`$– nor $`\mathrm{s}`$– are dominantly in phase or in opposite phase, and the two radial velocities have different numbers of radial nodes. Therefore neither the index nor the $`\mathrm{o}`$/$`\mathrm{s}`$ label bear any reliable information about the properties of the modes. This behaviour is possible because the that this eigenvalue problem is not of Sturm–Liouville type except in the non–stratified case.

An interesting fact to notice in tables 3 to 5 is that the fundamental modes (fo and fs) are the lowest frequency modes in the spectrum, in other words there are no g–modes present (which usually lie far below the f–mode) in these superfluid models. This confirms the numerical findings by Lee (1995) and the local analysis of Andersson & Comer (2001a).

The absence of g-modes can be made clearer when acoustic modes and surface gravity modes are filtered out. The latter modes are easily removed by suppressing surface motions and imposing therefore $`W_X=0`$ at the star surface. Acoustic modes, on the other hand, are filtered out by using the so-called anelastic approximation which makes an expansion in powers of the Brunt-Väisälä frequency (see Dintrans & Rieutord 2001; Rieutord & Dintrans 2002). Using this approximation mass conservation now reads

$$(\rho _n\delta v_\mathrm{n})=0\mathrm{and}(\rho _c\delta v_\mathrm{c})=0$$
(94)
instead of (45) and (46). In the case we have considered, i.e. that of no entrainment ($`\alpha =0`$), the equations of motions read:

$`(\delta \stackrel{~}{\mu }^\mathrm{n}+\delta \mathrm{\Phi })`$ $`=`$ $`i\omega \delta v_\mathrm{n},`$ (95)
$`(\delta \stackrel{~}{\mu }^\mathrm{c}+\delta \mathrm{\Phi })`$ $`=`$ $`i\omega \delta v_\mathrm{c},`$ (96)
$`^2\delta \mathrm{\Phi }`$ $`=`$ $`k_\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}+k_\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c},`$ (97)
Using (94), we eliminate the velocities and are left with the system:

$`\left(\rho _n(\delta \stackrel{~}{\mu }^\mathrm{n}+\delta \mathrm{\Phi })\right)`$ $`=`$ $`0,`$ (98)
$`\left(\rho _c(\delta \stackrel{~}{\mu }^\mathrm{c}+\delta \mathrm{\Phi })\right)`$ $`=`$ $`0,`$ (99)
$`^2\delta \mathrm{\Phi }`$ $`=`$ $`k_\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}+k_\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c},`$ (100)
where we see that the mode frequency has disappeared. Using the boundary conditions, it turns out that the only solution is $`\delta \stackrel{~}{\mu }^\mathrm{c}=\delta \stackrel{~}{\mu }^\mathrm{n}=\delta \mathrm{\Phi }=0`$, thus showing that no eigenmode exists when acoustic and surface gravity modes are filtered out.

However, the presence of g–modes due to chemical composition gradients in normal–fluid neutron star models has been pointed out by Reisenegger & Goldreich (1992), and their possibly observable excitation in a coalescing binary neutron star has been discussed by Reisenegger & Goldreich (1994) and Lai (1994). We will see in Sect. 5.4 that these predicted composition g–modes do indeed appear in the normal–fluid case. In principle the presence or absence of these modes could therefore be used as a possibly observable indicator for superfluidity in neutron stars.

#### 5.3.2 The effect of coupling by entrainment

In this section we study the dependence of the mode-frequencies and properties on the coupling by entrainment. Obviously, we only need to specify *one* entrainment function, $`\epsilon _\mathrm{c}`$ say, as $`\epsilon _\mathrm{n}`$ is then determined by (6). Because the uncertainties and differences of the “realistic” models for $`\epsilon _\mathrm{c}`$ provided by nuclear physics calculations so far are still considerable, we chose the simplest entrainment model, namely $`\epsilon _\mathrm{c}=\epsilon `$ being a constant. The value of this constant $`\epsilon `$ can be related to the proton effective mass $`m_\mathrm{p}^{}`$ (22), and is roughly constrained from the nuclear physics calculations (Chao et al. 1972; Sjöberg 1976; Baldo et al. 1992; Borumand et al. 1996) to lie in the range $`0.3\epsilon 0.7`$. We nevertheless consider the broader range between $`0.8\epsilon 0.8`$ to demonstrate the qualitative behaviour more clearly. This will also show that the “locally uncoupled” case $`\epsilon =0`$ (considered in the previous section) is not special, contrary to what one might have expected.

The results for the mode-frequencies as functions of $`\epsilon `$ for the three background models are represented in Fig. 5. In the case of the non–stratified model I, we observe the predicted (sect. 4.3) decoupling, and in particular the independence of the “ordinary”–type modes of entrainment. Because of this decoupling the respective frequencies of the two mode families can simply cross each other when $`\epsilon `$ is varied. In the generic stratified models (model II and III), the modes of the doublets are coupled and avoided crossings result when mode–frequencies come to close to each other, as also found recently by Andersson et al. (2002). In this process of avoided crossing the two modes seem to exchange some of their respective properties of being dominantly “co-” or “counter-moving”, as can be seen in Fig. 6), and they also can exchange their number of radial nodes, as we see in the avoided crossing of the p$`{}_{}{}^{\mathrm{s}}{}_{1}{}^{}`$ and p$`{}_{}{}^{\mathrm{o}}{}_{2}{}^{}`$ in Fig. 7.

Another important conclusion can be drawn from Fig. 5, namely that the “locally uncoupled” case $`\epsilon =0`$ discussed in the previous section does not represent a special case in any respect, because the two fluids are *always* coupled through $`\delta \mathrm{\Phi }`$. The effect of $`\epsilon `$ is simply to change the coupling, but no configuration is completely uncoupled. On can see in Fig. 5 that several avoided crossings happen practically at $`\epsilon =0`$, which is the case in particular for the p6-modes of model II presented in Fig. 4.

### 5.4 The one–fluid case: recovering the g–modes

Following the discussion in sect. 3.2, the one–fluid case is defined by $`\delta v_\mathrm{n}=\delta v_\mathrm{c}\delta v`$. We only have one Euler equation in this case, which in the harmonic decomposition (54) has the two components

$`x_\mathrm{n}\delta \stackrel{~}{\mu }_{}^{\mathrm{n}}{}_{}{}^{}+x_\mathrm{c}\delta \stackrel{~}{\mu }_{}^{\mathrm{c}}{}_{}{}^{}+\delta \mathrm{\Phi }^{}`$ $`=`$ $`{\displaystyle \frac{i\omega }{r}}W,`$ (101)
$`x_\mathrm{n}\delta \stackrel{~}{\mu }^\mathrm{n}+x_\mathrm{c}\delta \stackrel{~}{\mu }^\mathrm{c}+\delta \mathrm{\Phi }`$ $`=`$ $`i\omega V.`$ (102)
These two equations replace (57)–(60), while the remaining equations of this system are unchanged (subject to the substitutions $`W_\mathrm{n}=W_\mathrm{c}=W`$ and $`V_\mathrm{n}=V_\mathrm{c}=V`$). The eigen-frequencies of this system are shown in table 6, where we see the presence of composition g–modes as expected for all models with stratification. This is consistent with the prediction by Reisenegger & Goldreich (1992) and the numerical findings of Lee (1995).

We also see that in the non–stratified model I, the one–fluid frequencies and modes correspond exactly to the corresponding “ordinary”-type solutions of the two–fluid case (see table 5 and Fig. 2), as would be expected from the separability of the system as discussed in Sect. 4.3.

We note that the perfect fluid modes of stratified models generally have $`\delta \beta 0`$, because adiabatic oscillations generally drive fluid elements out of equilibrium, only in the non–stratified case (model I) is $`\delta \beta =0`$ strictly satisfied.

The absence/presence of $`g`$–modes in superfluid/normal fluid models might seem somewhat surprising but one can get a better intuitive understanding by considering the physical origin of these $`g`$–modes: a radially displaced fluid element will remain close to mechanical (pressure) equilibrium with its surroundings, but its respective values of $`\rho _\mathrm{n}`$ and $`\rho _\mathrm{c}`$ will generally differ from the surroundings (i.e. when $`x_\mathrm{c}^{}0`$) and therefore (via the equation of state) its total density will also be different, resulting in a buoyant restoring force and a corresponding oscillation mode (in unstable models this restoring force will actually drive the fluid element still further away from its initial position, leading to convection). In the simple (cold) superfluid models considered here, each fluid element of either fluid ($`\mathrm{n}`$ or $`\mathrm{c}`$) is only characterised by a *single* quantity, namely $`\rho _\mathrm{n}`$ or $`\rho _\mathrm{c}`$. Displacing an element of fluid $`\mathrm{n}`$, say, will therefore result not only in mechanical equilibrium ($`\stackrel{~}{\mu }^\mathrm{n}`$), but also in buoyant equilibrium. This can by seen by expressing its density at the new position as $`\rho _\mathrm{n}=\rho _\mathrm{n}(\stackrel{~}{\mu }^\mathrm{n},\rho _\mathrm{c})`$. The fluid $`\mathrm{c}`$ was not displaced, therefore not only $`\stackrel{~}{\mu }^\mathrm{n}`$ but also $`\rho _\mathrm{c}`$ of the fluid element are identical to the background values, and so is $`\rho _\mathrm{n}`$. If we had allowed for an additional comoving quantity like entropy $`s`$, we would expect to find g–modes driven by a stratification in $`s/\rho _\mathrm{n}`$.

It is intriguing to see that the absence of the $`g`$–modes in superfluid models is accompanied by an apparent doubling of acoustic modes, but it is not obvious to establish a link between these different classes of modes as we are currently not aware of a continuous transition from a two–fluid to a one–fluid model (either the two fluids are locked together or they are not).

## 6 Conclusions

In this paper we have tried to clarify the qualitative properties of the eigenmode spectrum of superfluid neutron stars, using a simple two–fluid model. We have shown the important — previously somewhat overlooked — role of stratification for these modes. The picture has been found to be more complex than previous studies have suggested, and some of the earlier conclusions have been shown to apply only for non–stratified models. In particular, one can not generally talk about two distinct families of “superfluid” and “ordinary” modes. The system of equations describing two–fluid modes can not be separated in the case of stratified stars, and its solutions have no direct correspondence to the eigenmodes of the one–fluid system. The two–fluid modes are generally neither co– nor counter-moving, rather *all* of them are characterised by non–zero amplitudes of relative velocity $`\delta \mathrm{\Delta }`$, deviation of chemical equilibrium $`\delta \beta `$ and total density perturbation $`\delta \rho `$. Also the order of the mode does not necessarily correspond to the number of radial nodes (as seen in Fig. 4), which is possible because the system is not of the Sturm–Liouville type. We have further confirmed earlier findings about the absence of g–modes in these superfluid models (Lee 1995; Andersson & Comer 2001a), as well as the appearance of avoided crossings between mode frequencies when changing the entrainment parameter (Andersson et al. 2002).

Given the radical difference and richer structure of the oscillations of superfluid neutron star models as compared to the simple perfect fluid models, we think that much future effort is needed to further clarify these properties and evaluate possibly observable consequences. The respective absence and presence of g–modes in these two different models is a striking example of such a potentially observable indicator of superfluidity in neutron stars. However, many more physical effects have to be taken into account in order to achieve a more realistic description of superfluid neutron stars, namely the inclusion of vortex–forces and beta reactions, both of which will lead to dissipation. Furthermore, an “envelope” or an elastic crust should be included, and maybe most importantly, the effects of rotation and magnetic field, which add new restoring forces and result in a much richer spectrum of modes. Eventually, for a realistic study of oscillations of superfluid neutrons stars, one needs to work in a generally relativistic framework, as pioneered by Comer et al. (1999) and Andersson et al. (2002). This step is crucially important also for the assessment of the gravitational radiation emitted by these modes, and their stability/instability via the CFS mechanism (Friedman & Schutz 1978).

###### Acknowledgements.

We thank N. Andersson, D. Langlois, D. Gondek, J.L. Zdunik, H. Beyer, G.L. Comer and B. Dintrans for very valuable discussions in the early stages of this work. We also thank D.I. Jones and N. Andersson for a careful reading of the manuscript. RP acknowledges support from the EU Programme ’Improving the Human Research Potential and the Socio-Economic Knowledge Base’ (Research Training Network Contract HPRN-CT-2000-00137).