A&A 393, 949-963 (2002)
DOI: 10.1051/0004-6361:20021049

Adiabatic oscillations of non-rotating superfluid neutron stars

R. Prix¹ - M. Rieutord^2,3

1 - Department of Mathematics, University of Southampton, Southampton SO17 1BJ, UK
2 - Laboratoire d'Astrophysique de Toulouse, Obs. Midi-Pyrénées, 14 avenue E. Belin, 31400 Toulouse, France
3 - Institut Universitaire de France

Received 2 May 2002/ Accepted 12 July 2002

Abstract
We present results concerning the linear (radial and non-radial) oscillations of non-rotating superfluid neutron stars in Newtonian physics. We use a simple two-fluid model to describe the superfluid neutron star, where one fluid consists of the superfluid neutrons, while the second fluid contains all the remaining constituents (protons, electrons). The two fluids are assumed to be "free'' in the sense of absence of vortex-mediated forces like mutual friction or pinning, but they can be coupled by the equation of state, in particular by entrainment. We calculate numerically the eigen-frequencies and -modes of adiabatic oscillations, neglecting beta-reactions that would lead to dissipation. We find a doubling of all acoustic-type modes (f-modes, p-modes), and confirm the absence of g-modes in these superfluid models. We show analytically and numerically that only in the case of non-stratified background models (i.e. with no composition gradient) can these acoustic modes be separated into two distinct families, which are characterized by either co- or counter-moving fluids respectively, and which are sometimes referred to as "ordinary'' and "superfluid'' modes. In the general, stratified case, however, this separation is not possible, and these acoustic modes can not be classified as being either purely "ordinary'' or "superfluid''. We show how the properties of the two-fluid modes change as functions of the coupling by entrainment. We find avoided mode-crossings for the stratified models, while the crossings are not avoided in the non-stratified, separable case. The oscillations of normal-fluid neutron stars are recovered as a special case simply by locking the two fluids together. In this effective one-fluid case we find the usual singlet f- and p-modes, and we also find the expected g-modes of stratified neutron star models.

Key words: stars: neutron - stars: oscillations

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.

Most studies of neutron star oscillations are still 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; 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 emphasizing 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 characterized 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 used 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 specialize 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. Section 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 summarize 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 protons, electrons and entropy (and generally further particles like muons etc.). We denote the particle number densities for neutrons, protons and electrons as $n_{\rm n}$, $n_{\rm p}$ and $n_{\rm e}$ respectively, and we use sfor the entropy density. The velocities of the two fluids are $\vec{v}_{\rm n}$ for the neutron fluid, and $\vec{v}_{\rm c}$ for the fluid of comoving constituents, the relative velocity $\vec{\Delta}$ between the two fluids is therefore

$$\vec{\Delta} \equiv \vec{v}_{\rm c} - \vec{v}_{\rm 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_{\rm e}= n_{\rm 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 ( ${\rm n}\rightleftharpoons {\rm p}+{\rm e}+ \bar{\nu}$) between the two fluids, the total mass in the reaction has to be conserved in a consistent Newtonian description. Therefore we set

$$m_{\rm b}\equiv m_{\rm n} \equiv m_{\rm p} + m_{\rm e}.$$ (3)

The respective mass densities of the two fluids can now be written as

$$\rho_{\rm n} = m_{\rm b}n_{\rm n},{\quad{\rm and}\quad} \rho_{\rm c} = m_{\rm b}n_{\rm p} ,$$ (4)

and the total mass density $\rho$ is simply $\rho=\rho_{\rm n} + \rho_{\rm c}$.

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

$${\rm d}{\cal E}= T ~{\rm d} s + \tilde{\mu}^{\rm c} ~{\rm d... ...{\rm n} ~{\rm d} \rho_{\rm n} +\alpha ~{\rm d}\vec{\Delta}^2 .$$ (5)

This differential defines the dynamic quantities, namely the temperature T, the specific chemical potentials $\tilde{\mu}^{\rm c}$, $\tilde{\mu}^{\rm n}$, and the entrainment $\alpha$ as the thermodynamic conjugates of the kinematic quantities s, $\rho_{\rm c}$, $\rho_{\rm n}$ and $\vec{\Delta}$. The specific chemical potentials $\tilde{\mu}$ are related to the more common definition of the chemical potentials $\mu$ via $\tilde{\mu}\equiv \mu/m_{\rm b}$. We note that $\mu^{\rm 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^{\rm c}=\mu^{\rm p} + \mu^{\rm e}$, where $\mu^{\rm p}$ and $\mu^{\rm 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 $\vec{\Delta}$between the two fluids, which characterizes 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 $\varepsilon_{\rm n}$ and $\varepsilon_{\rm c}$ as

$$\varepsilon_{\rm c} \equiv {2\alpha\over \rho_{\rm c}}, {\q... ... n}} = {\rho_{\rm c} \over \rho_{\rm n}} \varepsilon_{\rm c} .$$ (6)

The "generalized'' pressure P is introduced as the Legendre-conjugate of the energy density, namely by the usual thermodynamic relation

$${\cal E}+ P = \rho_{\rm n} \tilde{\mu}^{\rm n} + \rho_{\rm c} \tilde{\mu}^{\rm c} + s T,$$ (7)

which results in the total differential of the pressure function $P(\tilde{\mu}^{\rm c},\tilde{\mu}^{\rm n},T,\vec{\Delta}^2)$:

$${\rm d} P = \rho_{\rm n} ~ {\rm d}\tilde{\mu}^{\rm n} + \rh... ...{\mu}^{\rm c} + s~ {\rm d} T - \alpha~ {\rm d}\vec{\Delta}^2 .$$ (8)

This generalized pressure P can be seen to reduce to the usual definition of the pressure of a perfect fluid in the case of $\vec{\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{\mathit{\Gamma}}_{\rm s}= {\vec{f}}\cdot\vec{\Delta} + \beta {\mathit{\Gamma}}_{\rm n} ,$$ (9)

where $\vec{f}$ is the mutual force between the two fluids, and ${\mathit{\Gamma}}_{\rm s}$ and ${\mathit{\Gamma}}_{\rm n}$ are the creation rates of entropy and neutrons respectively, i.e.

$${\mathit{\Gamma}}_{\rm s}\equiv \partial_{\rm t} s + {\nabla}... ...} \rho_{\rm n} + {\nabla}\cdot(\rho_{\rm n} \vec{v}_{\rm n}) ,$$ (10)

while the proton creation rate ${\mathit{\Gamma}}_{\rm c}$ has to satisfy ${\mathit{\Gamma}}_{\rm c}=-{\mathit{\Gamma}}_{\rm n}$ for baryon conservation. The quantity $\beta$ in (9) characterizes the deviation from chemical equilibrium and its explicit expression is found as

$$\beta \equiv \tilde{\mu}^{\rm c} - \tilde{\mu}^{\rm n} - {1\over2}\left(1 - 2\varepsilon_{\rm n}\right)~{\vec{\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 $\tilde{\mu}^{\rm c} - \tilde{\mu}^{\rm n}$ in the case of relative motion. The mutual force density $\vec{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. ${\mathit{\Gamma}}_{\rm s}\ge0$, is satisfied by (9), we can write the neutron creation rate ${\mathit{\Gamma}}_{\rm n}$ and the mutual force $\vec{f}$ in the form

$${\mathit{\Gamma}}_{\rm n} = \Xi \beta, \hspace{1.7cm}\quad{\rm with} \quad \Xi\ge0 ,$$ (12)

$${\vec{f}} = \eta \vec{\Delta} + \vec{\kappa}\times \vec{\Delta} , \quad{\rm with} \quad \eta \ge 0 ,$$ (13)

where the non-negative functions $\Xi$ and $\eta$ govern the beta-reaction rate and the friction force, while the vector $\vec{\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 ${\rm n}\rightleftharpoons {\rm c}$, has $\Xi=0$, and free vortices correspond to $\eta=0$ and $\vec{\kappa}=0$.

The momentum equation for the superfluid neutrons is given by

$$(\partial_{\rm t} + \vec{v}_{\rm n}\cdot{\nabla})(\vec{v}_{\r... ... n} {\nabla}v_{\rm n}^i = {1\over\rho_{\rm n}}{\vec{f}},\quad$$ (14)

while the equation for the normal fluid reads as

 

$$(\partial_{\rm t} + \vec{v}_{\rm c}\cdot{\nabla})( \vec{v}_{\rm c... ...de{\mu}^{\rm c} +\Phi) - \Delta_i \varepsilon_{\rm c} {\nabla}v_{\rm c}^i \quad \quad$$

$$+ {s\over\rho_{\rm c}}\nabla T = -{1\over \rho_{\rm c}} {\vec{f}}... ...on_{\rm n} \right) {{\mathit{\Gamma}}_{\rm n} \over \rho_{\rm c}} \vec{\Delta}.$$ (15)

The gravitational potential $\Phi$ is related to the mass densities $\rho_{\rm n}$ and $\rho_{\rm c}$ via the Poisson equation

$$\nabla^2 \Phi = 4\pi G (\rho_{\rm n} + \rho_{\rm c}).$$ (16)

   
2.2 The static equilibrium background

We consider a static background star, so we set

$$\vec{v}_{\rm n}=\vec{v}_{\rm p}=\vec{\Delta}=0, {\quad{\rm an... ...{\rm n}={\mathit{\Gamma}}_{\rm c}={\mathit{\Gamma}}_{\rm s}=0,$$ (17)

which by (13) implies the vanishing of the mutual force, i.e. ${\vec{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 = \tilde{\mu}^{\rm c} - \tilde{\mu}^{\rm n} = 0.$$ (18)

These equilibrium conditions reduce the equation of motion (14) to

$${\nabla}\tilde{\mu}^{\rm n} = {\nabla}\tilde{\mu}^{\rm c} = - {\nabla}\Phi ,$$ (19)

and with (15) this also implies that the background star is in thermal equilibrium, i.e. ${\nabla}T=0$.

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

$$\tilde{\mu}''(r) + {2\over r} \tilde{\mu}'(r) = - 4\pi G \rho(r),$$ (20)

where we have introduced the equilibrium chemical potential $\tilde{\mu}\equiv \tilde{\mu}^{\rm n} = \tilde{\mu}^{\rm c}$, and the prime (') denotes the radial derivative ${\rm d}/{\rm d}r$.

The equation of state allows one to relate the equilibrium chemical potential $\tilde{\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 $\vec{\Delta}$, we can separate the "bulk'' equation of state from the entrainment by expanding ${\cal E}(s, \rho_{\rm n}, \rho_{\rm c}; \vec{\Delta}^2)$ in terms of $\vec{\Delta}$, i.e. by writing

$${\cal E}= {\cal E}^{(0)}(s, \rho_{\rm n}, \rho_{\rm c}) + ... ...\rm n},\rho_{\rm c})\vec{\Delta}^2 + {\cal O}(\vec{\Delta}^4),$$ (21)

where the quantities ${\cal E}^{(0)}$ and $\alpha^{(0)}$ are evaluated at zero relative velocity $\vec{\Delta}$. The background equation of state ${\cal E}^{(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_{\rm p}^*$ (which is generally a function of the densities), in the form

$$2\alpha= \rho_{\rm c} \left( 1 - {m_{\rm p}^*\over m_{\rm b}}\right)\cdot$$ (22)

The dimensionless entrainment functions $\varepsilon_{\rm n}$ and $\varepsilon_{\rm c}$ can then be expressed according to (6) as

$$\varepsilon_{\rm c} = 1 - {m_{\rm p}^* \over m_{\rm b}} ,{\... ... n} = {\rho_{\rm c}\over \rho_{\rm n}} ~ \varepsilon_{\rm c} .$$ (23)

   
3 Linearized 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. Linearizing the equations of motion (14) and (15) yields

$$\rho_{\rm n} \partial_{\rm t} ( \delta \vec{v}_{\rm n} + \va... ...bla}(\delta\tilde{\mu}^{\rm n} + \delta\Phi) =\delta{\vec{f}},$$ (24)

$$\rho_{\rm c} \partial_{\rm t}( \delta \vec{v}_{\rm c} - \vare... ...}^{\rm c} + \delta\Phi) =-\delta{\vec{f}} -s {\nabla}\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{\vec{f}} = \eta \delta\vec{\Delta} + \vec{\kappa}\times \delta\vec{\Delta},$$ (26)

and the linearized energy conservation (9) and (12) together with the condition of baryon conservation lead to

$$\partial_{\rm t}~ \delta s + {\nabla}\cdot(s ~ \delta\vec{v}_{\rm c}) = 0,$$ (27)

$$\partial_{\rm t}~ \delta \rho_{\rm n} + {\nabla}\cdot(\rho_{\... ...Xi ~(\delta\tilde{\mu}^{\rm c} - \delta\tilde{\mu}^{\rm n}),$$ (28)

$$\partial_{\rm t}~ \delta \rho_{\rm c} + {\nabla}\cdot(\rho_{\... ...\Xi ~(\delta\tilde{\mu}^{\rm c} - \delta\tilde{\mu}^{\rm n}).$$ (29)

The perturbed Poisson Eq. (16) reads (in natural units) as

$$\nabla^2 \delta\Phi = \delta \rho_{\rm n} + \delta \rho_{\rm c} .$$ (30)

The system is closed by specification of the mutual force functions $\eta$ and $\vec{\kappa}$, the transfusion function $\Xi$, and an equation of state which allows to express the dynamical quantities $\delta T$, $\delta\tilde{\mu}^{\rm n}$ and $\delta\tilde{\mu}^{\rm c}$ in terms of the kinematic variables $\delta s$, $\delta\rho_{\rm n}$ and $\delta\rho_{\rm 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 ${\rm 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 $\vec{v}_{\rm n} = \vec{v}_{\rm c}$, and therefore

$$\delta\vec{v}_{\rm n} = \delta\vec{v}_{\rm c} = \delta\vec{v}, \quad{\quad\Longrightarrow\quad} \delta\vec{\Delta} = 0.$$ (31)

By linking the two constituents together, the degrees of freedom have been reduced by three, and instead of the individual momentum Eqs. (24) and (25), now only the sum of momenta can be required to be conserved, i.e.

$$\rho ~\left( \partial_{\rm t} \delta\vec{v} + {\nabla}\delt... ...c} {\nabla}\delta\tilde{\mu}^{\rm c} + s{\nabla}\delta T = 0 .$$ (32)

We introduce the notation

$$\rho_{\rm c} = x_{\rm c}\rho,\quad \rho_{\rm n} = x_{\rm n}\rho,{\quad{\rm and}\quad} s = \tilde{s}\rho,$$ (33)

for the proton and neutron fractions $x_{\rm c}$ and $x_{\rm n}$, and the specific entropy $\tilde{s}$, which allows us to rewrite the one-fluid equation of motion (32) in the slightly more familiar form

$$\partial_{\rm t} \delta\vec{v} + {\nabla}\left( {\delta P \... ...a~ {\nabla}x_{\rm c}+ \delta T ~ {\nabla}\tilde{s}\right] = 0,$$ (34)

where we used the fact that the total pressure differential ${\rm d} P$ of (8) in this perfect fluid case reduces to

$${\rm d} P = \rho_{\rm n} ~ {\rm d}\tilde{\mu}^{\rm n} + \rho_{\rm n} ~ {\rm d}\tilde{\mu}^{\rm c} + s~ {\rm d} T.$$ (35)

We can now compare (34) 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

$$\partial_{\rm t} \delta\vec{v} + {\nabla}\left( {\delta P \... ...ht) + {\cal C}(r)~ ({\nabla}\cdot\delta\vec{v})~\vec{A} = 0,$$ (36)

where ${\cal C}(r)$ is a function of the background, and $\vec{A}$ is the so-called Schwarzschild discriminant that is responsible for the presence of g-modes. By comparing Eqs. (34) and (36), we see that $\vec{A}$ will be non-zero (indicating the presence of g-modes) whenever ${\nabla}\tilde{s}\not=0$ or ${\nabla}x_{\rm c}\not=0$. This reflects the well-known fact that any type of stratification, either in specific entropy $\tilde{s}$ or in the chemical composition $x_{\rm 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 Eqs. (24)-(30), we need a specific model for the mutual force $\vec{f}$ and the transfusion $\Xi$, in addition to an equation of state of the form ${\cal E}(\rho_{\rm n},\rho_{\rm c}, 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. $\Xi=0$) of free fluids (meaning ${\vec{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=0and T=0. The resulting simplified system of equations is

$$% \partial_{\rm t} \delta \rho_{\rm n} + {\nabla}\cdot\left(\rho_{\rm n} ~\delta\vec{v}_{\rm n}\right) = 0,$$ (37)

$$% \partial_{\rm t} \delta \rho_{\rm c} + {\nabla}\cdot\left(\rho_{\rm c} ~\delta\vec{v}_{\rm c}\right) = 0,$$ (38)

$$% \partial_{\rm t} \left[ (1-\varepsilon_{\rm n}) ~\delta\vec... ... {\nabla}\left(\delta\tilde{\mu}^{\rm n} + \delta\Phi\right),$$ (39)

$$\partial_{\rm t} \left[ (1-\varepsilon_{\rm c})~ \delta\vec{v... ... {\nabla}\left(\delta\tilde{\mu}^{\rm c} + \delta\Phi\right),$$ (40)

$$% \nabla^2 \delta\Phi = \delta \rho_{\rm n} + \delta \rho_{\rm c}.$$ (41)

We point out that this system of equations is identical to the one used in Andersson & Comer (2001a), which was obtained as the 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_{\rm X}$ to $\delta\tilde{\mu}^{\rm X}$ (with the constituent index notation ${\rm X}={\rm n},~{\rm c}$) to linear order, namely

$$\begin{array}{r l} \delta\rho_{\rm n} =& {\cal S}_{{\rm n}... ...S}_{{\rm c}{\rm c}}~ \delta\tilde{\mu}^{\rm c},\\ \end{array}$$ (42)

where the symmetric "structure matrix'' ${\cal S}_{{\rm X}{\rm Y}}$ is defined as

$$\left({\cal S}^{-1}\right)_{{\rm X}{\rm Y}} \equiv {\partia... ...ho_{\rm Y}},\qquad {\rm X},{\rm Y}\in \{{\rm n},{\rm c}\}\cdot$$ (43)

Due to the "dual'' role of the pressure P with respect to the energy density ${\cal E}$, as seen in (7) and (8), we can equivalently express ${\cal S}_{{\rm X}{\rm Y}}$ as

$${\cal S}_{{\rm X}{\rm Y}} = {\partial \rho_{\rm X}\over \pa... ...partial \tilde{\mu}^{\rm X}~\partial \tilde{\mu}^{\rm Y}}\cdot$$ (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 ${\cal S}_{{\rm n}{\rm c}}$ in (42), while the second is due to the entrainment terms  $\varepsilon_{\rm X}$.

The background quantities $\tilde{\mu}^{\rm X}$ can be seen in (19) to behave like the gravitational potential $\Phi$; this means in particular that their gradient is always finite, even at the surface. Therefore $\delta\tilde{\mu}^{\rm X}$ is finite everywhere, contrary to $\delta\rho_{\rm X}$ which can diverge at the surface when $\rho'\rightarrow-\infty$. This is seen from the relation $\Delta\rho = \delta\rho + \xi^r~\rho'(r)$ between the Lagrangian perturbation $\Delta\rho$ and the Eulerian $\delta\rho$, for a radial displacement $\xi^r$. On physical grounds $\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'\rightarrow-\infty$ and $\xi^r\not=0$ 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\rightarrow\infty$ 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 $\Delta\rho$ is still perfectly regular. If one wanted to impose that $\delta\rho$ should be bounded everywhere (as done by Lindblom & Mendell 1994), then this situation would be inverted and the Lagrangian quantity $\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\tilde{\mu}^{\rm X}$ instead of the potentially diverging $\delta\rho_{\rm X}$, by using (42) to substitute for $\delta\rho_{\rm X}$. We note that the coefficients ${\cal S}_{{\rm X}{\rm Y}}$ in this expression will generally diverge (or vanish) at the surface, depending on the equation of state, and which reflects the behaviour of $\delta\rho_{\rm X}$. The system of Eqs. (37)-(41) for eigenmode solutions of the form $\delta Q(\vec{x},t) = \delta Q(\vec{x})~{\rm e}^{i{\omega} t}$ now yields

$${\nabla}\cdot( \rho_{\rm n} \delta\vec{v}_{\rm n}) = - i{\ome... ... {\cal S}_{{\rm n}{\rm c}} \delta\tilde{\mu}^{\rm c} \right],$$ (45)

$${\nabla}\cdot( \rho_{\rm c} \delta\vec{v}_{\rm c}) = - i{\ome... ... {\cal S}_{{\rm c}{\rm c}} \delta\tilde{\mu}^{\rm c} \right],$$ (46)

$${\nabla}( \delta\tilde{\mu}^{\rm n} + \delta\Phi ) = - i{\ome... ...rm n} + \varepsilon_{\rm c} ~ \delta\vec{v}_{\rm c} \right],$$ (47)

$${\nabla}( \delta\tilde{\mu}^{\rm c} + \delta\Phi ) = - i{\ome... ...rm c} + \varepsilon_{\rm n} ~ \delta\vec{v}_{\rm n} \right],$$ (48)

$$\nabla^2 \delta\Phi = k_{\rm n} ~ \delta\tilde{\mu}^{\rm n} + k_{\rm c}~ \delta\tilde{\mu}^{\rm c},$$ (49)

where we have introduced the convenient "structure vector'' $k_{\rm X}$, which is defined as

$$% k_{\rm X}\equiv \sum_{{\rm Y}={\rm n},{\rm c}} S_{{\rm X}{\rm Y}},\qquad {\rm X}\in\{{\rm n},~{\rm c}\}\cdot$$ (50)

Using the definition (43) of the matrix ${\cal S}_{{\rm X}{\rm Y}}$together with the equilibrium condition ${\tilde{\mu}^{\rm n}}{}' = {\tilde{\mu}^{\rm c}}{}' = \tilde{\mu}'$we can now write

$$\rho_{\rm X}' = k_{\rm X}\tilde{\mu}'.$$ (51)

For a spherically symmetric background, we can separate the radial and angular dependence and obtain solutions with definite quantum numbers l and $\vert m\vert\le l$ using the ansatz

$$\delta\Phi(r,\theta,\varphi) = \delta\Phi(r) ~Y_l^m(\theta,\varphi),$$ (52)

$$\delta\tilde{\mu}^{\rm X}(r,\theta,\varphi) = \delta\tilde{\mu}^{\rm X}(r) ~Y_l^m(\theta,\varphi),$$ (53)

$$\delta\vec{v}_{\rm X}(r,\theta,\varphi) = {W_{\rm X}(r)\over ... ...{R}_l^m + V_{\rm X}(r)\vec{S}_l^m - i U_{\rm X}(r)\vec{T}_l^m,$$ (54)

where the spherical harmonics $Y_l^m(\theta,\varphi)$ are the eigenfunctions of $r^2 \nabla^2 Y_l^m = - l(l+1)~ Y_l^m$, and $\vec{R}$, $\vec{S}$ and $\vec{T}$ form the orthogonal harmonic basis, defined as

$$\vec{R}_l^m \equiv Y_l^m {\nabla}r,\quad \vec{S}_l^m \equ... ...abla}Y_l^m ,\quad \vec{T}_l^m \equiv {\nabla}\times \vec{R},$$ (55)

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_{\rm n} W_{\rm n}\right)' - {l(l+1)} \rho_{\rm n... ... + {\cal S}_{{\rm n}{\rm c}} \delta\tilde{\mu}^{\rm c}\right],$$ (56)

$$\left(r \rho_{\rm c} W_{\rm c}\right)' - {l(l+1)}\rho_{\rm c}... ...+ {\cal S}_{{\rm c}{\rm c}} \delta\tilde{\mu}^{\rm c} \right],$$ (57)

$$-r\left({\delta\tilde{\mu}^{\rm n}}' + \delta\Phi'\right) = i... ...{\rm n} + \varepsilon_{\rm n} W_{\rm c}) \right],\qquad\mbox{}$$ (58)

$$-r\left({\delta\tilde{\mu}^{\rm c}}' + \delta\Phi'\right) = i... ...n_{\rm c}) W_{\rm c} + \varepsilon_{\rm c} W_{\rm n}) \right],$$ (59)

$$-\delta\tilde{\mu}^{\rm n} - \delta\Phi = i{\omega} \left[ (1... ...on_{\rm n}) V_{\rm n} + \varepsilon_{\rm n} V_{\rm c})\right],$$ (60)

$$-\delta\tilde{\mu}^{\rm c} - \delta\Phi = i{\omega} \left[ (1... ...on_{\rm c}) V_{\rm c} + \varepsilon_{\rm c} V_{\rm n})\right],$$ (61)

$$\left(r^2\delta\Phi'\right)' - {l(l+1)}~\delta\Phi = r^2~(k_... ...ta\tilde{\mu}^{\rm n} + k_{\rm c} ~\delta\tilde{\mu}^{\rm c}).$$ (62)

The axial velocity component $U_{\rm 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 Eqs. (60) and (61) only hold for l>0, because $\vec{S}_0^0 = 0$.

   
4.2 Boundary and regularity conditions

4.2.1 At the center

It can be shown that the representation (52)-(54) of a regular physical quantity requires the following asymptotic behaviour of the radial functions as $r \rightarrow 0$,

$$\delta\tilde{\mu}^{\rm X}\sim \delta\Phi \sim W_{\rm X}\sim... ...m X}\sim {\cal O}(r^l),\qquad {\rm X}\in \{{\rm n},~{\rm c}\},$$ (63)

where ${\cal O}(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. (58) and (59) we see that $W_{\rm X}$ must vanish as ${\cal O}(r)$ at the origin, while from Eqs. (56) and (57) we can derive the further regularity requirement

$$\left(1+ r{\rho_{\rm X}'\over\rho_{\rm X}}\right) W_{\rm X}+ r W_{\rm X}' - l(l+1) V_{\rm X}\sim {\cal O}(r^2).$$ (64)

Regularity of solutions of Poisson's Eq. (62) requires

$$\delta\Phi \sim {\cal O}(r^2)~\quad{\rm for}\quad l> 0,$$ (65)

while for radial oscillations (l=0) only $\delta\Phi' \sim {\cal O}(r)$ is required. These constraints are automatically satisfied by (63) for $l\ge2$, but they are stronger requirements than (63) 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 $\Phi$, which results (e.g. see Ledoux & Walraven 1958) in the boundary condition

$$\delta\Phi'(R) + {l+1\over R}\delta\Phi(R) = -\xi^r~ \rho(R-),$$ (66)

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_{r\rightarrow R-}\rho(r)$. In the present work we will only consider stars with vanishing density at the surface i.e. $\rho(r)\rightarrow0$. The conservation Eqs. (56) and (57) contain a (regular) singularity at the surface r=R, because $\rho'/\rho$ diverges when $\rho\rightarrow0$. We first rewrite these equations as

$$(r W_{\rm X})' - l(l+1) V_{\rm X}= -{r\over \rho_{\rm X}}\l... ...\cal S}_{{\rm X}{\rm Y}} \delta\tilde{\mu}^{\rm Y}\right].\;\;$$ (67)

Regularity of the solution therefore requires the following asymptotic behaviour at the surface:

$$\rho_{\rm X}' W_{\rm X}+ i{\omega} r \sum {\cal S}_{{\rm X}{\rm Y}} \delta\tilde{\mu}^{\rm Y}\sim {\cal O}(\rho_{\rm X}).$$ (68)

   
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}{r l r l} \delta\vec{\Delta} \equiv & \delta\v... ...^{\rm c} + x_{\rm n}~\delta\tilde{\mu}^{\rm n} .\\ \end{array}$$ (69)

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, characterized 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 (69), 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 (50) of $k_{\rm n}$ and $k_{\rm c}$, and defining $k\equiv k_{\rm n} + k_{\rm c}$, we can express the density perturbations $\delta\rho_{\rm X}$ as

$$\begin{array}{r l} \delta\rho_{\rm n} = & k_{\rm n} \delta... ...\rm c}} + x_{\rm n}k_{\rm c}\right)\delta\beta,\\ \end{array}$$ (70)

and therefore

$$\delta\rho = k \delta\tilde{\mu}- (x_{\rm c}k_{\rm n} - x_{\rm n}k_{\rm c})~ \delta\beta.$$ (71)

In the case of a spherically symmetric background considered here, we can use (51) to express

$$x_{\rm c}k_{\rm n} - x_{\rm n}k_{\rm c} = - k{\rho\over\rho'} ~x_{\rm c}',$$ (72)

and using these variables and relations, we now rewrite the eigenmode Eqs. (37)-(41) in the form

$$k ~\partial_{\rm t} \delta\tilde{\mu}+ {\nabla}\cdot(\rho~\de... ... ~\left[k{\rho\over\rho'}~\partial_{\rm t} \delta\beta\right],$$ (73)

$$\partial_{\rm t} \delta\vec{v} + {\nabla}(\delta\tilde{\mu}+ ... ...a\Phi) = \;\;~x_{\rm c}' ~\left[ \delta\beta~\vec{e}^r\right],$$ (74)

$$\nabla^2 \delta\Phi - k ~ \delta\tilde{\mu}= -x_{\rm c}' ~\left[k{\rho\over\rho'}~\delta\beta\right],$$ (75)

$${\cal K}~ \partial_{\rm t} \delta\beta+ {\nabla}\cdot\left( \... ...}~\partial_{\rm t} \delta\tilde{\mu}+ \rho~\delta v^r \right],$$ (76)

$$\left(1-{2\alpha\rho\over\rho_{\rm n} \rho_{\rm c}}\right) \p... ...al_{\rm t}~ {\delta\vec{\Delta}} + {\nabla}\delta\beta= \;\;0,$$ (77)

where $\vec{e}^r$ is the radial basis vector, and we introduced the abbreviation ${\cal K}\equiv x_{\rm c}^2 k_{\rm n} + x_{\rm n}^2 k_{\rm c} - {\cal S}_{{\rm n}{\rm c}}$.

We see that the "ordinary''-type of motion ( $\delta\tilde{\mu}$, $\delta\vec{v}$, $\delta\Phi$) does not decouple from the "superfluid''-type variables ( $\delta\vec{\Delta}$, $\delta \beta$) whenever there is stratification, i.e. when $x_{\rm c}'\not=0$! This can be understood as follows: while a non-zero relative velocity $\delta\vec{\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 (69) 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_{\rm c}$ to a constant, we can separate the equations into two decoupled sets. One system describes "ordinary'' modes, namely

$$k \partial_{\rm t} \delta\tilde{\mu}+ {\nabla}\cdot( \rho\delta\vec{v} ) = 0,$$ (78)

$$\partial_{\rm t} \delta\vec{v} + {\nabla}(\delta\tilde{\mu}+ \delta\Phi) = 0,$$ (79)

$$\nabla^2 \delta\Phi - k ~ \delta\tilde{\mu}= 0,$$ (80)

which are seen to be independent of the entrainment $\varepsilon$ as well as of the coupling through the "bulk'' equation of state, i.e. ${\cal S}_{{\rm n}{\rm c}}$. The second system of equations governs the "superfluid'' modes and reads as

$$\left(k - {{\cal S}_{{\rm n}{\rm c}}\over x_{\rm n}x_{\rm c}}... ...m t} \delta\beta+ {\nabla}\cdot(\rho~\delta\vec{\Delta}) = 0,$$ (81)

$$\left(1-{2\alpha\rho\over\rho_{\rm n}\rho_{\rm c}}\right)~ \partial_{\rm t}~ \delta\vec{\Delta} + {\nabla}\delta\beta= 0.$$ (82)

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. ${\cal S}_{{\rm n}{\rm c}}$, but they are completely decoupled from the gravitational perturbation $\delta\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\vec{v},\delta\tilde{\mu}, \delta\Phi, 0, 0)$ and "superfluid'' modes $(0,0,0,\delta\vec{\Delta},\delta\beta)$. One of our numerical models (see next section) has a constant proton fraction $x_{\rm 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 "generalized polytrope'', defined as

$${\cal E}^{(0)}(\rho_{\rm n}, \rho_{\rm c}) = {\k_{\rm n}}~\... ...^{\gamma_{\rm n}}+ {\k_{\rm c}}~ \rho_{\rm c}^{\gamma_{\rm c}}$$ (83)

which simply consists of the sum of two ordinary polytropes. For regularity of the chemical potentials (5) in the limit $\rho_{\rm n}\rightarrow0$ and $\rho_{\rm c}\rightarrow0$, the polytropic indices must satisfy ${\gamma_{\rm n}}\ge1$ and ${\gamma_{\rm c}}\ge1$. This equation of state allows the explicit inversion

$$\rho_{\rm n}(\tilde{\mu}^{\rm n}) = \left( {\tilde{\mu}^{\rm ... ...rm c} \over {\k_{\rm c}}{\gamma_{\rm c}}} \right)^{N_{\rm c}},$$ (84)

where we have introduced $N_{\rm X}\equiv 1/({\gamma_{\rm X}}-1)$. We see that in chemical equilibrium, i.e. $\tilde{\mu}^{\rm n}=\tilde{\mu}^{\rm c}$, the two fluids share a common outer surface. The equilibrium proton fraction $x_{\rm c}$ can be expressed as

$$x_{\rm c}= \left( 1 + { ({\gamma_{\rm c}}{\k_{\rm c}})^{N_{... ...rm n}} } ~ \tilde{\mu}^{N_{\rm n} - N_{\rm c}} \right)^{-1} ,$$ (85)

which shows that the proton fraction is constant whenever ${\gamma_{\rm n}}={\gamma_{\rm c}}$, while the behaviour in the case of different indices falls into the two categories:

$$0 \stackrel{\tilde{\mu}\rightarrow\infty}{\longleftarrow} x_{\rm ... ...w 0}{\longrightarrow} 1,\quad{\rm for}\quad {\gamma_{\rm n}}< {\gamma_{\rm c}},$$ (86)

$$1 \stackrel{\tilde{\mu}\rightarrow\infty}{\longleftarrow} x_{\rm ... ...w 0}{\longrightarrow} 0,\quad{\rm for}\quad {\gamma_{\rm n}}> {\gamma_{\rm c}}.$$ (87)

   
5.2 Calculating the background models

The equilibrium background solution is determined by Eq. (20), together with the regularity requirement $\tilde{\mu}'(0)=0$, and the boundary condition of vanishing pressure at the surface, i.e. $\tilde{\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(\tilde{\mu})$ in chemical equilibrium, and therefore Eq. (20) can be written as the following nonlinear eigenvalue problem,

$$\tilde{\mu}'' + {2\over r} \tilde{\mu}' = \lambda ~\rho(\tilde{\mu}),$$ (88)

where the (dimensionless) eigenvalue $\lambda$ is given by

$$\lambda = 4\pi G {\rho_0 \over \tilde{\mu}(\rho_0)} R^2 .$$ (89)

The eigenvalue $\lambda$ represents the a-priori unknown radius of the star, and thereby the actual values of the 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 (88). This can be done by solving in step k+1 the following linear eigenvalue problem,

$$\tilde{\mu}''_{k+1} + {2\over r}\tilde{\mu}'_{k+1} = -\lamb... ...(\tilde{\mu}_k)\over \tilde{\mu}_k}\right)~ \tilde{\mu}_{k+1},$$ (90)

where $\tilde{\mu}_k$ is the solution 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. $\sim$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 which we use in the numerical code for the eigenmodes. Even when using a different resolution for the eigenmode-calculation, the Chebychev expansion provides a canonical ${\cal C}^\infty$ interpolation which allows us to easily "re-grid'' the background solution.

	Figure 1: Density profiles and proton fraction $x_{\rm c}$ of the background models I, II and III as defined in Table 1. $\rho _0$ is the central density. In this figure we have divided the neutron density by 5 in order to obtain similar magnitudes for the different curves.
Open with DEXTER

	Figure 2: The first three quadrupolar (l=2) eigenmode doublets (${\rm o}$)/(${\rm s}$) for the background model I. In the case of the (${\rm o}$)-modes, the two fluids are comoving to numerical precision, therefore the $\delta v_{\rm n}$ and $\delta v_{\rm c}$ curves can not be seen separately.
Open with DEXTER

	Figure 3: The first three quadrupolar (l=2) eigenmode doublets (${\rm o}$)/(${\rm s}$) for the background model II. The value of $\delta \beta$ in the lower left figure has been multiplied by 100, in order to make it visible.
Open with DEXTER

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 (83).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.

 

 

Table 1: Parameters describing the background models I, II and III, based on the two-constituent polytrope (83). Further shown are the central mass density $\rho _0$, the central proton fraction  $x_{\rm c}(0)$, the total mass M and the radius R. ^a The units of the coefficients  $k_{\rm X}$ are $c^2{\rho _{\rm nuc}}^{1-{\gamma _{\rm X}}}$, where $\rho _{\rm nuc}=1.66\times 10^{14}$ g cm^-3.

	${\gamma_{\rm n}}$	${\gamma_{\rm c}}$	${\k_{\rm n}}^a$	${\k_{\rm c}}^a$	${\rho_0\over\rho_{\rm nuc}}$	$x_{\rm c}(0)$	${M\over M_{\odot}}$	${R\over{\rm km}}$
I	2.0	2.0	0.01	0.09	10.83	0.10	1.414	10.7
II	2.5	2.1	0.01	0.20	3.20	0.10	1.440	14.3
III	1.9	1.7	0.01	0.09	18.65	0.10	1.412	9.4

We note that model I is the only non-stratified model, i.e. it has a constant proton fraction $x_{\rm c}$ (and finite density gradients at the surface), while models II and III have a non-zero composition gradient, i.e. $x_{\rm c}'\not=0$, as expected from (85).

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 (1976) units $\sqrt{G M/R^3}$ (variants thereof, like those used by Lindblom & Mendell 1994 only differ by a constant factor) and the "geometric'' units c³/G M typically used in general relativity (Comer et al. 1999; Andersson et al. 2002).

 

 

Table 2: Conversion factors for our frequency units ( $\sqrt {4\pi G\rho _0}$) into different systems of units (SI, "Cox'', and "geometrized'' units) for models I, II and III.

	I	II	III
kHz	38.8318047	21.1116405	50.9444446
$\sqrt{G M/R^3}$	3.14159265	2.61829313	3.40801371
c3/GM	0.270521502	0.149723677	0.354013053

   
5.3 The two-fluid oscillation modes

The eigenmode Eqs. (56)-(62) 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-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 $\varepsilon_{\rm n}=\varepsilon_{\rm 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\Phi$ and (62). 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 "o'' and the higher frequency one as "s''. The pairs of p-modes are indexed in the order of increasing frequency. We emphasize 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 "${\rm o}$'' modes the two fluids are comoving, resulting in a non-zero $\delta\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 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\Phi=0$, while the fluids are driven out of chemical equilibrium, i.e. $\delta\beta\not=0$. The number of radial nodes in $\delta\vec{v}^r$ is the same for the o and s modes, and corresponds exactly to their index. All these results confirm the analytic predictions for non-stratified models in Sect. 4.3.

 

 

Table 3: Frequency spectrum of radial eigenmodes (l=0) for models I, II and III in natural units ${\omega }_0=\sqrt {4\pi G \rho _0}$.

l=0	${\omega}_{\rm I}$	${\omega}_{\rm II}$	${\omega}_{\rm III}$
$f^{\rm o}$	0.616 801 012	0.860 501 159	0.539 820 916
$f^{\rm s}$	0.825 395 141	1.004 218 360	0.713 202 951
$p_1^{\rm o}$	1.272 153 763	1.650 440 676	1.114 032 342
$p_1^{\rm s}$	1.398 557 067	1.780 808 966	1.186 350 906
$p_2^{\rm o}$	1.855 852 617	2.326 264 710	1.582 750 279
$p_2^{\rm s}$	1.949 822 942	2.573 698 536	1.675 075 045
$p_3^{\rm o}$	2.418 457 671	2.985 429 110	2.018 290 454
$p_3^{\rm s}$	2.493 326 179	3.352 297 316	2.169 556 932
$p_4^{\rm o}$	2.970 977 248	3.638 960 946	2.445 611 455
$p_4^{\rm s}$	3.033 169 591	4.122 179 788	2.658 934 809

	Figure 4: The p6 eigenmode doublet (o)/(s) (for l=2) for the background model II. The radial velocity $\delta v^r_{\rm n}$ has been multiplied by 10 for better visibility. Neither mode is co- or counter-moving, and the number of radial nodes differs between the two fluids.
Open with DEXTER

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 "o'' modes are not comoving at all (only the $f^{\rm o}$ is nearly comoving), and they have non-zero $\delta \beta$, while the "s'' modes have non-zero $\delta\rho$ and $\delta\Phi$. One can not say either that the relative amplitude of $\delta \beta$ would be different between the o- and s-cases, as Lindblom & Mendell (1994) wrongly induced from the properties of the $f^{\rm o}$-mode, the only eigenmode they presented.

In the case of the low-order modes presented in Fig. 3, the number of radial nodes still seems to correspond to the index, and the fluids are roughly in opposite phase in the s-modes, while they are approximately in phase for the o-modes. Even this, however, is not true in general, as can be seen in Fig. 4 which shows the higher order p₆-modes. In this case neither o- nor 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 o/s label bear any reliable information about the properties of the modes. This behaviour is possible because this eigenvalue problem is not of Sturm-Liouville type except in the non-stratified case.

 

 

Table 4: Frequency spectrum of dipolar eigenmodes (l=1) for models I, II and III in natural units ${\omega }_0=\sqrt {4\pi G \rho _0}$. ^a The zero frequency modes correspond to the analytic result of a constant displacement field $\vec{\xi}$. This solution can not be produced by our code because the equations have been expressed in terms of velocities instead of displacements.

l=1	${\omega}_{\rm I}$	${\omega}_{\rm II}$	${\omega}_{\rm III}$
$f^{\rm o}$	0a	0a	0a
$f^{\rm s}$	0.389 835 134	0.440 176 989	0.387 577 390
$p_1^{\rm o}$	0.898 011 966	1.191 707 859	0.798 353 611
$p_1^{\rm s}$	1.040 570 747	1.280 907 920	0.895 985 737
$p_2^{\rm o}$	1.518 621 841	1.930 140 941	1.317 638 486
$p_2^{\rm s}$	1.621 023 028	2.105 369 668	1.384 086 979
$p_3^{\rm o}$	2.099 218 641	2.609 579 518	1.770 672 472
$p_3^{\rm s}$	2.179 811 073	2.908 118 483	1.887 282 252
$p_4^{\rm o}$	2.662 623 840	3.274 392 031	2.206 080 421
$p_4^{\rm s}$	2.729 000 908	3.692 629 614	2.385 023 708

 

 

Table 5: Frequency spectrum of quadrupolar eigenmodes (l=2) for models I, II and III in natural units ${\omega }_0=\sqrt {4\pi G \rho _0}$.

l=2	${\omega}_{\rm I}$	${\omega}_{\rm II}$	${\omega}_{\rm III}$
$f^{\rm o}$	0.390 550 961	0.424 294 338	0.376 662 787
$f^{\rm s}$	0.526 990 499	0.604 904 572	0.516 636 947
$p_1^{\rm o}$	1.101 827 434	1.423 800 575	0.980 798 266
$p_1^{\rm s}$	1.206 881 695	1.514 568 856	1.039 415 772
$p_2^{\rm o}$	1.723 444 868	2.164 156 378	1.478 705 720
$p_2^{\rm s}$	1.806 873 473	2.380 237 199	1.557 015 636
$p_3^{\rm o}$	2.310 315 782	2.858 007 106	1.932 213 886
$p_3^{\rm s}$	2.379 342 400	3.200 741 221	2.072 313 682
$p_4^{\rm o}$	2.879 785 468	3.533 905 837	2.372 551 143
$p_4^{\rm s}$	2.938 448 588	3.997 594 591	2.576 350 147

Another interesting fact to notice in Tables 3-5 is that the fundamental modes ($f^{\rm o}$ and $f^{\rm s}$) 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

$${\nabla}\cdot(\rho_{\rm n}\delta\vec{v}_{\rm n}) = 0\quad {\rm and}\quad {\nabla}\cdot(\rho_{\rm c}\delta\vec{v}_{\rm c}) = 0$$ (91)

instead of (45) and (46). In the case we have considered, i.e. that of no entrainment ($\alpha=0$), the equations of motion read:

$${\nabla}( \delta\tilde{\mu}^{\rm n} + \delta\Phi ) = - i{\omega} \delta\vec{v}_{\rm n} ,$$ (92)

$${\nabla}( \delta\tilde{\mu}^{\rm c} + \delta\Phi ) = - i{\omega} \delta\vec{v}_{\rm c} ,$$ (93)

$$\nabla^2 \delta\Phi = k_{\rm n} ~ \delta\tilde{\mu}^{\rm n} + k_{\rm c}~ \delta\tilde{\mu}^{\rm c}.$$ (94)

Using (91) we eliminate the velocities and are left with the system:

$${\nabla}\cdot \left(\rho_{\rm n}{\nabla}(\delta\tilde{\mu}^{\rm n} + \delta\Phi ) \right)= 0,$$ (95)

$${\nabla}\cdot \left(\rho_{\rm c}{\nabla}(\delta\tilde{\mu}^{\rm c} + \delta\Phi ) \right)= 0,$$ (96)

$$\nabla^2 \delta\Phi = k_{\rm n} ~ \delta\tilde{\mu}^{\rm n} + k_{\rm c}~ \delta\tilde{\mu}^{\rm c},$$ (97)

where we see that the mode frequency has disappeared. Using the boundary conditions, it turns out that the only solution is $\delta\tilde{\mu}^{\rm c}=\delta\tilde{\mu}^{\rm n}=\delta\Phi=0$, thus showing that no eigenmode exists when acoustic and surface gravity modes are filtered out.

Figure 5: Eigenmode frequencies ${\omega }$ (for l=2) as functions of entrainment $\varepsilon$ for models I, II and III. The full and dashed lines represent the frequencies of the modes denoted respectively by o and s in Table 5. The avoided crossings marked by the full and the dashed box are represented in terms of the eigenfunctions in Figs. 6 and 7. We note that for models II and III all crossings shown are avoided.

Open with DEXTER

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.

Figure 6: Avoided crossing between the $p_1^{\rm o}$ and the $p_1^{\rm s}$ modes of Model II (for l=2), indicated by the full box in Fig. 5. The three columns show the corresponding eigenfunctions for three different values of entrainment. The small arrows indicate the "evolution'' of the eigenmode when increasing  $\varepsilon$. The $\delta v_{\rm n}^r$ curves in the two lower left figures have been magnified by a factor of 10.

Open with DEXTER

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, $\varepsilon_{\rm c}$ say, as $\varepsilon_{\rm n}$ is then determined by (6). Because the uncertainties and differences of the "realistic'' models for $\varepsilon_{\rm c}$ provided by nuclear physics calculations so far are still considerable, we chose the simplest entrainment model, namely $\varepsilon_{\rm c}=\varepsilon$ being a constant. The value of this constant $\varepsilon$ can be related to the proton effective mass $m_{\rm p}^*$ by (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\leq\varepsilon\leq 0.7$. We nevertheless consider the broader range between $-0.8\leq\varepsilon\leq 0.8$ to demonstrate the qualitative behaviour more clearly. This will also show that the "locally uncoupled'' case $\varepsilon=0$ (considered in the previous section) is not special, contrary to what one might have expected.

	Figure 7: Avoided crossing between the $p_1^{\rm s}$ and the $p_2^{\rm o}$ modes of Model II (l=2), indicated by the hatched box in Fig. 5. The three columns show the corresponding eigenfunctions for three different values of entrainment. The small arrows indicate the "evolution'' of the eigenmode when increasing  $\varepsilon$.
Open with DEXTER

The results for the mode-frequencies as functions of $\varepsilon$ 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 $\varepsilon$ 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 too 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_1^{\rm s}$ and $p_2^{\rm o}$ in Fig. 7.

Another important conclusion can be drawn from Fig. 5, namely that the "locally uncoupled'' case $\varepsilon=0$ discussed in the previous section does not represent a special case in any respect, because the two fluids are always coupled through $\delta\Phi$. The effect of $\varepsilon$ 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 $\varepsilon=0$, which is the case in particular for the p₆-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\vec{v}_{\rm n} = \delta\vec{v}_{\rm c} \equiv \delta\vec{v}$. We only have one Euler equation in this case, which in the harmonic decomposition (55) has the two components

$$x_{\rm n}~ {\delta\tilde{\mu}^{\rm n}}' + x_{\rm c}~{\delta\tilde{\mu}^{\rm c}}' + \delta\Phi' -{i{\omega}\over r}~ W,$$ = (98)

$$x_{\rm n}~ \delta\tilde{\mu}^{\rm n}~ + x_{\rm c}~\delta\tilde{\mu}^{\rm c}~ + \delta\Phi~ -i{\omega} ~V.$$ = (99)

These two equations replace (58)-(61), while the remaining equations of this system are unchanged (subject to the substitutions $W_{\rm n}=W_{\rm c}=W$ and $V_{\rm n}=V_{\rm 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).

 

 

Table 6: The eigenvalue spectrum (for l=2) for the models I, II, and III. The non-stratified model I has no g-modes.

l=2	${\omega}_{\rm I}$	${\omega}_{\rm II}$	${\omega}_{\rm III}$
...	-	...	...
g4	-	0.012 105 268	0.011 489 709
g3	-	0.015 003 335	0.014 157 505
g2	-	0.019 814 997	0.018 492 575
g1	-	0.029 631 110	0.026 880 058
f	0.390 550 961	0.424 310 492	0.376 717 911
p1	1.101 827 434	1.477 988 230	0.988 324 062
p2	1.723 444 868	2.348 478 094	1.533 337 250
p3	2.310 315 782	3.163 853 143	2.050 348 013
p4	2.879 785 468	3.954 289 860	2.552 745 872
...	...	...	...

	Figure 8: The first few (l=2) g-, f- and p-modes for model II. We note that the g-modes are characterized by a very small radial velocity $\delta v^r$ at the surface, while it is very large in the case of p-modes.
Open with DEXTER

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\not=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_{\rm n}$ and $\rho_{\rm c}$ will generally differ from the surroundings (when $x_{\rm c}'\not=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 (${\rm n}$ or ${\rm c}$) is only characterized by a single quantity, namely $\rho_{\rm n}$ or $\rho_{\rm c}$. Displacing an element of fluid ${\rm n}$, say, will therefore result not only in mechanical equilibrium ( $\tilde{\mu}^{\rm n}$), but also in buoyant equilibrium. This can by seen by expressing its density at the new position as $\rho_{\rm n}=\rho_{\rm n}(\tilde{\mu}^{\rm n}, \rho_{\rm c})$. The fluid ${\rm c}$ was not displaced, therefore not only $\tilde{\mu}^{\rm n}$ but also $\rho_{\rm c}$ of the fluid element are identical to the background values, and so is $\rho_{\rm 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_{\rm 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 - and 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 characterized by non-zero amplitudes of relative velocity $\delta\vec{\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 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 absence of g-modes in the superfluid models is in strong contrast with the normal fluid models and 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 (Chandrasekhar 1970; 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).

References

 

• Andersson, N., & Comer, G. 2001a, MNRAS, 328, 1129 In the text NASA ADS  
• Andersson, N., & Comer, G. 2001b, Phys. Rev. Lett., 87, 241101 In the text  
• Andersson, N., Comer, G., & Langlois, D. 2002, unpublished, preprint [gr-qc/0203039] In the text  
• Andreev, A., & Bashkin, E. 1975, JETP, 42, 164 In the text  
• Baldo, M., Cugnon, J., Lejeune, A., & Lombardo, U. 1992, Nucl. Phys. A, 536, 349 In the text  
• Borumand, M., Joynt, R., & Kluzniak, W. 1996, Phys. Rev. C, 54, 2745 In the text NASA ADS  
• Carter, B. 1989, in Lecture Notes in Mathematics, vol. 1385, Relativistic Fluid Dynamics (Noto, 1987), ed. A. Anile, & M. Choquet-Bruhat (Springer-Verlag, Heidelberg), 1 In the text  
• Carter, B., Langlois, D., & Sedrakian, D. M. 2000, A&A, 361, 795 In the text NASA ADS  
• Chandrasekhar, S. 1970, Phys. Rev. Lett., 24, 611 In the text NASA ADS  
• Chao, N.-C., Clark, J., & Yang, C.-H. 1972, Nucl. Phys. A, 179, 320 In the text  
• Comer, G., Langlois, D., & Lin, L. M. 1999, Phys. Rev. D, 60, 104025 In the text  
• Cox, J. P. 1976, ARA&A, 14, 247 In the text NASA ADS  
• Cox, J. P. 1980, Theory of Stellar Pulsation (Princeton University Press) In the text  
• Dintrans, B., & Rieutord, M. 2001, MNRAS, 324, 635 In the text NASA ADS  
• Epstein, R. I. 1988, ApJ, 333, 880 In the text NASA ADS  
• Friedman, J. L., & Schutz, B. F. 1978, ApJ, 222, 281 In the text NASA ADS  
• Lai, D. 1994, MNRAS, 270, 611 In the text NASA ADS  
• Landau, L., & Lifshitz, E. 1959, Course of Theoretical Physics, vol. 6, Fluid Mechanics (Pergamon, Oxford) In the text  
• Langlois, D., Sedrakian, D. M., & Carter, B. 1998, MNRAS, 297, 1189 In the text NASA ADS  
• Ledoux, P., & Walraven, T. 1958, in Handbuch der Physik, ed. S. Flügge, vol. 51 (Berlin: Springer-Verlag), 353 In the text  
• Lee, U. 1995, A&A, 303, 515 In the text NASA ADS  
• Lindblom, L., & Mendell, G. 1994, ApJ, 421, 689 In the text NASA ADS  
• Lindblom, L., & Mendell, G. 1995, ApJ, 444, 804 In the text NASA ADS  
• Lindblom, L., & Mendell, G. 2000, Phys. Rev. D, 61, 104003 In the text  
• Link, B., Epstein, R. I., & Lattimer, J. 2000, in Stellar Astrophysics, Proceedings of the Pacific Rim Conference, ed. L. Cheng, H. Chau, K. Chan, & K. Leung (Kluwer Academic Publishers, Netherlands), 117 In the text  
• Prix, R. 2002, unpublished, preprint [physics/0209024] In the text  
• Prix, R., Comer, G., & Andersson, N. 2002, A&A, 381, 178 In the text NASA ADS  
• Reisenegger, A., & Goldreich, P. 1992, ApJ, 395, 240 In the text NASA ADS  
• Reisenegger, A., & Goldreich, P. 1994, ApJ, 426, 688 In the text NASA ADS  
• Rieutord, M. 1987, Geophys. Astrophys. Fluid Dyn., 39, 163 In the text  
• Rieutord, M., & Dintrans, B. 2002, unpublished, preprint [astro-ph/0206357] In the text  
• Sedrakian, A., & Wasserman, I. 2000, Phys. Rev. D, 63, 024016 In the text  
• Sjöberg, O. 1976, Nucl. Phys. A, 265, 511 In the text  
• Unno, W., Osaki, Y., Ando, H., Saio, H., & Shibahashi, H. 1989, Nonradial Oscillations of Stars, 2nd ed. (University of Tokyo Press) In the text