A&A 381, 178-196 (2002)
DOI: 10.1051/0004-6361:20011499

Slowly rotating superfluid Newtonian neutron star model with entrainment

R. Prix1 - G. L. Comer2 - N. Andersson1

1 - Department of Mathematics, University of Southampton, Southampton SO17 1BJ, UK
2 - Department of Physics, Saint Louis University, PO Box 56907, St. Louis, MO 63156-0907, USA

Received 11 July 2001 / Accepted 25 October 2001

We develop a formalism that can be used to model slowly rotating superfluid Newtonian neutron stars. A simple two-fluid model is used to describe the matter, where one fluid consists of the superfluid neutrons that are believed to exist in the inner crust and core of mature neutron stars, while the other fluid is a charge neutral conglomerate of the remaining constituents (crust nuclei, core superconducting protons, electrons, etc.). We include the entrainment effect, which is a non-dissipative interaction between the two fluids whereby a momentum induced in one of the fluids will cause part of the mass of the other fluid to be carried along. The equations that describe rotational equilibria (i.e. axisymmetric and stationary configurations) are approximated using the slow-rotation approximation; an expansion in terms of the rotation rates of the two fluids where only terms up to second-order are kept. Our formalism allows the neutrons to rotate at a rate different from that of the charged constituents. For a particular equation of state that is quadratic in the two mass-densities and relative velocities of the fluids, we find an analytic solution to the slow-rotation equations. This result provides an elegant generalisation to the two-fluid problem of the standard expressions for the one-fluid polytrope ${\cal E} \propto \rho^2$. The model equation of state includes entrainment and is general enough to allow for realistic values for, say, mass and radius of the star. It also includes a mixed term in the mass densities that can be related to "symmetry energy'' terms that appear in more realistic equations of state. We use the analytic solution to explore how relative rotation between the two fluids, the "symmetry energy'' term, and entrainment affect the neutron star's local distribution of particles, as well as global quantities as the Kepler limit, ellipticity, and moments of inertia.

Key words: stars: neutron - stars: rotation - hydrodynamics

1 Introduction

Pulsars were first discovered over 34 years ago (Hewish et al. 1968). Since then, nearly 1300 such rapidly rotating, highly magnetised neutron stars have been found, and as pointed out by Lorimer in his recent Living Reviews in Relativity (Lorimer 2001), 700 of these were found in the last 4 years alone. Without doubt, the overwhelming majority of these pulsars are old and cold (with core temperatures below 109 K). According to equation of state calculations, their interiors will contain superfluid neutrons, superconducting protons, a plasma of highly degenerate and ultra-relativistic electrons, and perhaps other more exotic particles (pions, hyperons, etc.) deep in their cores. In recent years, there has been a continued effort to model superfluid dynamics in neutron stars in both the Newtonian (Epstein 1988; Mendell & Lindblom 1991; Mendell 1991a, 1991b; Lindblom & Mendell 1994, 1995, 2000; Lee 1995; Prix 1999; Sedrakian & Wasserman 2000; Andersson & Comer 2001b) and general relativistic regimes (Carter 1989; Carter & Langlois 1995a, 1995b, 1998; Langlois et al. 1998; Comer et al. 1999; Andersson & Comer 2001a). This paper is aimed at improving our understanding of local and global properties of rotating Newtonian superfluid neutron stars.

The strongest evidence for superfluidity in the inner crust and core of a neutron star (see, for example, Sauls 1989 and references therein) is provided by the well-known glitch phenomenon (the occasional sudden spin-up of relatively young pulsars). Our confidence in an explanation based on the transfer of angular momentum between two loosely coupled components is bolstered by the fact that the neutrons and protons are described using the same many-body theory of Fermi liquids and BCS mechanism that has been so successful at describing superconductors (Pines & Nozières 1966). This being the case, one would expect superfluidity in neutron stars to share many of the well-established properties of laboratory superfluids and superconductors. For instance, it is known that in a mixture of two interpenetrating fluids there is a coupling that arises whereby the momentum of one of the liquids is not simply proportional to that liquid's velocity, rather, it is a linear combination of the velocities of both fluids. This is the so-called entrainment effect and it implies that when one liquid starts to flow then it will necessarily induce a momentum in the other constituent. Entrainment between protons and neutrons is a key component of models of neutron star superfluidity, and it is something that we will focus on in this study.

We develop a formalism for describing slowly rotating Newtonian superfluid neutron stars which allows the neutrons and protons to rotate at different rates. Our analysis is based on a Newtonian model that is the non-relativistic limit of a comprehensive model developed by Carter, Langlois and their collaborators for the general relativistic regime (Carter 1989; Carter & Langlois 1995a, 1995b, 1998; Langlois et al. 1998). Our model is simplified in the sense that it describes a superfluid neutron star in terms of only two fluids (we refer the reader to, for example, Comer et al. 1999 for justification). One fluid is composed of the superfluid neutrons, existing in the inner crust and core, and the other fluid is a charge neutral conglomerate of the remaining constituents (i.e. crust nuclei, core protons, and the crust and core electrons) that we will loosely refer to as the "protons.''

A further simplification of our model regards the vortices of the superfluid. Because the density of neutron vortices is expected to be about $10^2{-}10^5~ {\rm cm}^{-2}$ for typical pulsar rotation rates, it is reasonable to adopt a smooth averaged description of vorticity on a macroscopic scale. The resultant model then resembles a perfect fluid. There is, however, one important difference that arises because of the possible interaction between the vortex lattice and the second fluid. This is a dissipative effect usually referred to as "mutual friction''. In general, this mechanism tends to drive the two fluids towards co-rotation, but there are two extreme cases where a stationary description of a two-fluid star with different rotation rates would still be possible. The first case is (obviously) when the interaction is very small (the "free'' vortex limit). Somewhat surprisingly, similar results apply in the case when the interaction is very strong (resembling the "pinned'' vortex case in a solid crust, cf. Langlois et al. 1998). Although we could, in principle, allow for this second case in the present context (as has been done in an earlier study by Prix 1999), we will refrain from doing so and assume that vortex friction is negligible on the time-scales we are considering. This is done because the emphasis of the present work is on the detailed role of the entrainment effect. Including vortex friction would unnecessarily complicate the discussion and distract the attention away from the new piece of physics we have incorporated.

We will generalize the Chandrasekhar-Milne slow-rotation method (Chandrasekhar 1933; Milne 1923) to our two-fluid model. That is, we make an expansion in terms of the fluid rotation rates keeping only terms up to second order. We include entrainment, and consider the most general case where the neutrons do not corotate with the protons. This relative rotation of neutrons to protons is also present in the general relativistic slow rotation scheme of Andersson & Comer (2001a) (who provide an extended discussion as to why such a generalization is necessary), and the Newtonian slow rotation scheme of Prix (1999). However, although the formalism developed by Andersson & Comer does allow for entrainment, their numerical study did not include it, and Prix excluded entrainment altogether. As we wish to explore the effects of entrainment on rotational equilibria, and the current best model for entrainment is for the Newtonian regime (Borumand et al. 1996), it is natural for us to take up where Prix (1999) left off, and consider Newtonian models.

As an application of the slow-rotation formalism, we will consider a particular form of the equation of state (EOS) that is the most general quadratic form in the mass-densities and in the relative velocity. This EOS contains a mixed term in the neutron and proton mass densities that can be related to so-called "symmetry energy'' terms in realistic equations of state (i.e. terms that vanish when there are equal numbers of protons and neutrons). Remarkably, we find an analytic solution to the slow-rotation equations for this EOS. This solution has enough free parameters that realistic values for the neutron star mass, say, can be obtained. Using this solution we will show how relative rotation, entrainment, and the "symmetry energy'' term affect the neutron star's local distribution of particles, as well as global quantities like the ellipticity, the Kepler limit, and the moments of inertia.

Before concluding this introduction we should note one of the main motivations behind the present work. Once we have developed a framework for obtaining stationary superfluid models, we want to study their dynamics as manifested by the various modes of pulsation. This issue is of particular interest as it is known that several different modes of a rotating neutron star are generically unstable due to the emission of gravitational radiation (Friedman & Schutz 1978; Andersson & Kokkotas 2001). One can even speculate that pulsation modes unique to superfluids will be excited during a pulsar glitch and that the resultant gravitational radiation may be detected by a future generation of advanced detectors (Andersson & Comer 2001c). Our future aim is to investigate these various possibilities by using the models developed here as background for a linear perturbation calculation of the relevant pulsation modes of a rotating superfluid neutron star.

2 Canonical two-fluid hydrodynamics

In this section we present the canonical Newtonian description of a non-dissipative interacting two-fluid system, which can be derived either as the non-relativistic limit (Andersson & Comer 2001b) of the canonical covariant description developed by Carter & Langlois (1998), or directly from an analogous Newtonian variational principle (Prix 2000, 2001). Since we are mainly interested in describing a superfluid neutron star core we label the two fluids by indices "${\rm n}$'' and "${\rm p}$'', representing the neutrons and protons, respectively. The fundamental variables of our description are the respective particle number densities $n_{\rm n}$ and $n_{\rm p}$ and the corresponding particle currents $\vec{n}_{\rm n}$ and $\vec{n}_{\rm p}$. We will only consider situations where the particle numbers are conserved individually, so that we have

\begin{displaymath}\partial_t n_{\rm n} + \nabla_{i} n_{\rm n}^{i} = 0 , \quad {...
...} \quad
\partial_t n_{\rm p} + \nabla_{i} n_{\rm p}^{i} = 0 ,
\end{displaymath} (1)

where we sum over repeated spatial indices i,j = 1,2,3. The respective velocities $\vec{v}_{\rm n}$, $\vec{v}_{\rm p}$ and the mass densities $\rho_{\rm n}$ and $\rho_{\rm p}$ follow from

 \begin{displaymath}\vec{n}_{\rm n} = n_{\rm n} \vec{v}_{\rm n} , \quad
\vec{n}_{\rm p} = n_{\rm p} \vec{v}_{\rm p} ,
\end{displaymath} (2)


 \begin{displaymath}\rho_{\rm n} = m^{\rm n} n_{\rm n} , \quad
\rho_{\rm p} = m^{\rm p} n_{\rm p} , \quad
\rho = \rho_{\rm n} + \rho_{\rm p} ,
\end{displaymath} (3)

where $m^{\rm n}$ and $m^{\rm p}$ are the respective (fixed) masses per particle of the two fluids, and $\rho$ is the total mass density.

The internal energy density or "equation of state'' ${\cal E}(n_{\rm n}, n_{\rm p}, \vec{n}_{\rm n}, \vec{n}_{\rm p})$ of the two-fluid system must satisfy Galilean invariance (and isotropy, as we want to describe isotropic fluids), and therefore has to be of the form

 \begin{displaymath}{\cal E}= {\cal E}(n_{\rm n}, n_{\rm p}, \Delta^2) , \quad{\r...
\vec{\Delta} \equiv \vec{v}_{{\rm n}} - \vec{v}_{{\rm p}} .
\end{displaymath} (4)

This energy function defines the respective chemical potentials $\mu^{\rm n}$, $\mu^{\rm p}$ and the "entrainment function'' $\alpha$ as the conjugate variables to $n_{\rm n}$, $n_{\rm p}$ and $\Delta^2$, namely

 \begin{displaymath}{\rm d}{\cal E}= \mu^{\rm n}\,{\rm d}n_{\rm n} + \mu^{\rm p}\,{\rm d}n_{\rm p} + \alpha\, {\rm d}\Delta^2 ,
\end{displaymath} (5)

which can be regarded as the "first law of thermodynamics'' for this system. The canonical description of the two-fluid system is based on a convective variational principle (Prix 2000, 2001) analogous to the method used by Carter (1985) in the general relativistic context. The dynamics of the system is therefore described by a Lagrangian density ${\cal L}(n_{\rm n}, n_{\rm p}, \vec{n}_{\rm n}, \vec{n}_{\rm p})$, which has the usual form of "${\cal L}=$ kinetic energy - potential energy''. In the present context this means that

 \begin{displaymath}{\cal L}= {1 \over 2} \rho_{\rm n} \vec{v}_{\rm n}^{\,2} + {1...
\vec{v}_{\rm p}^{\,2} - \left({\cal E}+ \rho \Phi\right) ,
\end{displaymath} (6)

where $\Phi$ is the gravitational potential, which is related to the total mass density $\rho$ by the Poisson equation

 \begin{displaymath}\nabla^2 \Phi = 4\pi G \rho .
\end{displaymath} (7)

Variation of the Lagrangian density ${\cal L}$ defines the momenta per (fluid) particle, $\vec{p}^{\,{\rm n}}$ and $\vec{p}^{\,{\rm p}}$, and the "rest frame chemical potentials,'' $-p^{\rm n}_0$ and $-p^{\rm p}_0$, as the conjugate variables to the currents $\vec{n}_{\rm n}$, $\vec{n}_{\rm p}$ and particle densities $n_{\rm n}$, $n_{\rm p}$, i.e.

\begin{displaymath}{\rm d}{\cal L}= \vec{p}^{\,{\rm n}} \cdot {\rm d}\vec{n}_{\r...
p^{\rm n}_0 {\rm d}n_{\rm n} + p^{\rm p}_0 {\rm d}n_{\rm p}.
\end{displaymath} (8)

Using the explicit expression (6) for the Lagrangian density together with the "first law'' (5), we obtain the following expressions for $p_0^{\rm n}$ and $p_0^{\rm p}$,
$\displaystyle -p_0^{\rm n}$=$\displaystyle \mu^{\rm n} + \vec{v}_{\rm n} \cdot\vec{p}^{\,{\rm n}} - {1\over2} m^{\rm n}
\vec{v}^{\,2}_{\rm n} + m^{\rm n} \Phi ,$  
$\displaystyle -p_0^{\rm p}$=$\displaystyle \mu^{\rm p} + \vec{v}_{\rm p} \cdot\vec{p}^{\,{\rm p}} - {1\over2} m^{\rm p}
\vec{v}^{\,2}_{\rm p} + m^{\rm p} \Phi ,$ (9)

and the conjugate momenta $\vec{p}^{\,{\rm n}}$ and $\vec{p}^{\,{\rm p}}$ are related to the particle currents as
$\displaystyle \vec{p}^{\,{\rm n}}$ = $\displaystyle {\cal K}^{{\rm n}{\rm n}}\,\vec{n}_{\rm n} +
{\cal K}^{{\rm n}{\rm p}}\,\vec{n}_{\rm p} \,$  
$\displaystyle \vec{p}^{\,{\rm p}}$ = $\displaystyle {\cal K}^{{\rm n}{\rm p}}\,\vec{n}_{\rm n} +
{\cal K}^{{\rm p}{\rm p}}\,\vec{n}_{\rm p} ,$ (10)

where the symmetric "entrainment matrix'' ${\cal K}^{XY}$ has the following components:
    $\displaystyle {\cal K}^{{\rm n}{\rm n}} = {1\over n_{\rm n}^2}\left(m^{\rm n} n...
...m p}{\rm p}} = {1\over n_{\rm p}^2}\left(m^{\rm p} n_{\rm p} - 2\alpha\right) ,$  
$\displaystyle \quad$   $\displaystyle {\cal K}^{{\rm n}{\rm p}} = {2\alpha\over n_{\rm n} n_{\rm p}} \cdot$ (11)

The general form of the entrainment relation (10), which states that the momenta are linear combinations of the currents, is a consequence of the Galilean invariance of the internal energy ${\cal E}$. The most important difference from the case of a single fluid is that the particle momentum of each of the two fluids is in general not aligned with its respective current. From (10) we see that this deviation is caused by the off-diagonal matrix element ${\cal K}^{{\rm n}{\rm p}}$ of (11), which is proportional to the entrainment function $\alpha$ and therefore expresses the dependence of the internal energy density ${\cal E}$ on the relative velocity $\Delta$, as seen in (5). Intuitively this means that a particle current $\vec{n}_{\rm n}$ of the neutrons impinges some momentum on the protons and vice versa, due to the interaction between the two fluids. In the absence of such an interaction, i.e. for $\alpha=0$, the general entrainment relation (10) simply reduces to the usual single-fluid form, and we have $\vec{p}^{\,{\rm n}} = m^{\rm n} \vec{v}_{\rm n}$ and $\vec{p}^{\,{\rm p}} = m^{\rm p} \vec{v}_{\rm p}$.

The equations of motion for the two fluids, obtained from a Newtonian convective variational principle (see Prix 2000, 2001; Andersson & Comer 2001b), thus have the following form:

$\displaystyle n_{\rm n} \left(\partial_t \,\vec{p}^{\,{\rm n}} - \vec{\nabla} p...
...j} \left( \nabla_{j}\, \vec{p}^{\,{\rm n}} - \vec{\nabla} p^{\rm n}_{j}
\right)$ $\textstyle \,=\,$ 0 ,  
$\displaystyle n_{\rm p} \left(\partial_t \, \vec{p}^{\,{\rm p}} - \vec{\nabla} ...
...^{j} \left(\nabla_{j}\, \vec{p}^{\,{\rm p}} - \vec{\nabla} p^{\rm p}_{j}\right)$ $\textstyle \,=\,$ 0 . (12)

At this point, the two equations are "formally'' uncoupled, but the interaction between the two fluids enters through relations (9) and (10), which relate the kinematical quantities $n_{\rm n}$, $n_{\rm p}$, $\vec{n}_{\rm n}$ and $\vec{n}_{\rm p}$ to their respective conjugate dynamical quantities $p_0^{\rm n}$, $p_0^{\rm p}$, $\vec{p}^{\,{\rm n}}$ and $\vec{p}^{\,{\rm p}}$. In fact, the fluids will be coupled even in the case of two "non-interacting'' fluids, i.e. for an equation of state of the form ${\cal E}= {\cal E}_{\rm n}(n_{\rm n}) + {\cal E}_{\rm p}(n_{\rm p})$. This is simply because the two fluids live in the same gravitational potential, cf. Prix (1999).

3 Entrainment and effective masses

In the framework of condensed matter physics the description of interacting independent constituents, for example electrons moving through a "background'' of ions, is usually formulated in terms of an "effective'' or "dynamical'' mass (Ziman 1965), rather than in the entrainment formalism presented in the previous section. One reason for this might be that in the usual contexts one of the two constituents (the "background'') can usually be assumed to be at rest with respect to the observer. It is then convenient to formulate the dynamics of the second constituent in a form that resembles the vacuum expressions, with the free particle mass being replaced by an "effective mass'' that accounts for the interaction with the background.

The general definition of the effective mass is based on the particle energy spectrum, $E(\vec{p})$, say, of particles with momentum $\vec{p}$ moving relative to the non-vacuum background. The velocity $\vec{v}$ of the particles as a function of their momentum is then given by $\vec{v} = {\partial E
/ \partial \vec{p}}$, and the acceleration $\dot{\vec{v}}$ (where the dot represents the total time derivative) is linked to the force $\dot{\vec{p}}$ via a relation of the form

 \begin{displaymath}\dot{v}^i = m^{-1}_{ij} \dot{p}_j , \quad{\rm where} \quad
...}_{ij} \equiv {\partial^2 E \over \partial p_i \partial p_j} ,
\end{displaymath} (13)

which defines the effective mass tensor mij. In the two-fluid problem we are interested in, the background can usually be assumed to be isotropic, and therefore the particle energy spectrum will be of the form E = E(p2), which implies that

 \begin{displaymath}\vec{v} = \left(2{\partial E \over \partial p^2}\right) \, \vec{p} .
\end{displaymath} (14)

In this case the effective mass tensor of (13) can be expressed as

\begin{displaymath}m^{-1}_{ij} = \left(2 {\partial E \over \partial p^2}\right) ...
...artial^2 E \over \partial p^2 \partial p^2}\right) p_i p_j
\end{displaymath} (15)

In the limit of small momentum, where the energy will be mainly quadratic in p, we can neglect the second term and define an effective mass scalar m* in the usual way. This leads to

\begin{displaymath}{1 \over m^*} \equiv 2 {\partial E \over \partial p^2} \, \qu...
...ghtarrow \quad E(p) \approx {p^2 \over 2 m^*} +
{\rm const}.
\end{displaymath} (16)

We see from (14) that an equivalent form of this definition is

 \begin{displaymath}\vec{p} = m^* \, \vec{v} ,
\end{displaymath} (17)

which clarifies the link between the effective mass and the more general formalism of entrainment used in the present work. The key observation to make is that the effective mass description is formulated in the "rest-frame'' of the background. If we choose the background to be the "neutrons,'' we could define the neutron rest-frame by $\vec{n}_{\rm n} = 0$, and therefore the entrainment relation (10) in this frame reads

 \begin{displaymath}\vec{p}^{\,{\rm p}} = \left({\cal K}^{{\rm p}{\rm p}}n_{\rm p}\right) \, \vec{v}_{\rm p} .
\end{displaymath} (18)

By comparison with (17) we can relate the proton effective mass $m^{\rm p*}$ to the entrainment matrix by $m^{{\rm p}*} =
{\cal K}^{{\rm p}{\rm p}}n_{\rm p}$. Using the explicit expression (11) for ${\cal K}^{{\rm p}{\rm p}}$ in terms of the entrainment function $\alpha$, we find

 \begin{displaymath}2\alpha= n_{\rm p} ( m^{\rm p} - m^{{\rm p}*}) = \varepsilon \rho_{\rm p} ,
\end{displaymath} (19)

where we introduced the dimensionless "entrainment coefficient'' $\varepsilon$ as[*]

 \begin{displaymath}\varepsilon \equiv {m^{\rm p} - m^{{\rm p}*}\over m^{\rm p}}\cdot
\end{displaymath} (20)

Therefore the entrainment matrix ${\cal K}^{XY}$ is determined completely in terms of the effective mass of a single particle species, for example the protons in the present formulation. Of course, we also see (from (19) by symmetry) that the effective masses of the two species are not independent, because they obviously have to satisfy the consistency relation

 \begin{displaymath}m^{\rm n} - m^{{\rm n}*} = {n_{\rm p} \over n_{\rm n}} ( m^{\rm p} - m^{{\rm p}*} ).
\end{displaymath} (21)

However, as pointed out by Carter (2001, private communication), the above definition of the effective mass is not quite unique because we have to specify what we mean by the "rest-frame'' of the neutrons. Indeed, we have defined it by setting $\vec{n}_{\rm n} = 0$. However, we see from (10) that in this frame we generally have $\vec{p}^{\,{\rm n}}\not=0$. Therefore another equally viable choice of the "neutron rest-frame'' would be given by setting $\vec{p}^{\,{\rm n}}=0$. This would then lead to $\vec{n}_{\rm n}\not=0$. Based on this second choice, the entrainment relation (10) becomes

\begin{displaymath}\vec{p}^{\,{\rm p}} = \left( {\cal K}^{{\rm p}{\rm p}} - {({\...
... {\cal K}^{{\rm n}{\rm n}}}\right) n_{\rm p} \vec{v}_{\rm p} ,
\end{displaymath} (22)

leading to an alternative definition of the effective proton mass $m^{{\rm p}\char93 }$, which is not equivalent to that of Eq. (18). The corresponding relation between $\alpha$ and $m^{{\rm p}\char93 }$ is

 \begin{displaymath}2 \alpha= n_{\rm p} (m^{\rm p} - m^{{\rm p}\char93 }) \left\{...
...rm p} -
m^{{\rm p}\char93 }) \over m^{\rm n}} \right\}^{-1} ,
\end{displaymath} (23)

while the consistency relation instead of (21) now reads

 \begin{displaymath}m^{\rm n} - m^{{\rm n}\char93 } = {n_{\rm p} \over n_{\rm n}}...
...93 }/m^{\rm n} \over m^{{\rm p}\char93 }/m^{\rm p}} \right\} ,
\end{displaymath} (24)

corresponding to the choice of effective mass definition taken in the original approach by Andreev & Bashkin (1975). How does this ambiguity in the definition of the effective mass affect our attempt to model entrainment in superfluid neutron stars? Fortunately, it turns out that if we determine the entrainment $\alpha$either via (19) or (23), then the difference between the two (equally valid) definitions is numerically small in two interesting special cases:
neutron star matter: $m^{\rm p} - m^{{\rm p}\char93 }\sim m^{\rm n}$ and $n_{\rm p} \ll n_{\rm n}$,
electrons in a metal: n- = n+ and $m_- - m_-^\char93 \ll m_+$,
where in the second case n- and n+ represent the number densities of negative and positive charges respectively (i.e. electrons and protons), and m- and m+ the corresponding mass per particle. In these two cases Eq. (23) becomes approximately equal to the previous expression (19), and the two possible definitions lead to very similar results. The consistency relations (21) and (24), however, differ significantly in the case of neutron star matter (up to about a factor of 2), but unfortunately none of the current results in the literature on effective masses in neutron star matter seem to satisfy either of the two relations, making a well-defined identification even more troublesome.

There exist in the nuclear physics literature a few calculations of the proton effective mass at neutron star densities. The results are very much equation of state dependent, but one can extract some useful constraints on the range of values. Chao et al. (1972) find, for the range $\rho_{\rm nuc} \leq \rho \leq 2.5 \rho_{\rm nuc}$ (where $\rho_{\rm nuc}$ represents nuclear density, i.e. $\rho_{\rm nuc} = 2.7 \times 10^{14}$ g/cm3), that $m^*_{\rm p}/m_{\rm p} \approx 0.5{-}0.6$; Sjöberg (1976) finds, for $\rho_{\rm nuc} \leq \rho$, that $m^*_{\rm p}/m_{\rm p} \approx 0.3{-}0.6$; and more recently, Baldo et al. (1992) find, for the same range as Chao et al., that $m^*_{\rm p}/m_{\rm p} \approx 0.7$. Thus, we see that the proton effective mass apparently can range over the values $0.3 \leq m^*_{\rm p}/m_{\rm p} \leq 0.7$. The calculations also show that the effective proton mass varies slowly with the density. This means that one can, as a rough but reasonable approximation, assume that the effective proton mass remains constant throughout the core of a neutron star. This is, in fact, the case for the analytic solution that will be discussed later.

4 Stationary two-fluid neutron star models

Stationary two-fluid configurations are characterised by $\partial_t \vec{p}^{\,{\rm n}}=0$ and $\partial_t \vec{p}^{\,{\rm p}}=0$, and therefore the equations of motion (12) take the form

$\displaystyle \nabla_{i} \left[\mu^{\rm n} + m^{\rm n} \Phi -{1\over2} m^{\rm n...
...n}_{j} \nabla_{i} v_{\rm n}^{j} + v_{\rm n}^{j} \nabla_{j} p^{\rm n}_{i}\right)$ $\textstyle \,=\,$ 0 ,  
$\displaystyle \nabla_{i} \left[\mu^{\rm p} + m^{\rm p} \Phi -{1\over2} m^{\rm p...
...p}_{j} \nabla_{i} v_{\rm p}^{j} + v_{\rm p}^{j} \nabla_{j} p^{\rm p}_{i}\right)$ $\textstyle \,=\,$ 0 , (25)

where we have inserted the explicit expressions (9) for $p^{\rm n}_0$ and $p^{\rm p}_0$. We consider stationary configurations where both fluids are rotating uniformly with rotation rates $\Omega_{\rm n}$ and  $\Omega_{\rm p}$ about a common axis (otherwise the system would not be stationary),

 \begin{displaymath}v_{\rm n}^{i} = \Omega_{\rm n} \varphi^{i} , {\quad{\rm and}\quad}
v_{\rm p}^{i} = \Omega_{\rm p} \varphi^{i} ,
\end{displaymath} (26)

where $\varphi^{i}$ is the assumed axial symmetry vector with norm $(\varphi^{i} \varphi_{i})^{1/2} = \varpi$, and $\varpi$ is the distance from the rotation axis. In spherical coordinates $(r,\,\theta,\,\varphi)$ we, of course, have $\varpi = r\sin\theta$. As both fluids are flowing in the same direction $\varphi^{i}$, we can write

\begin{displaymath}v_{\rm p}^{i} = {\Omega_{\rm p} \over \Omega_{\rm n}} \,v_{\rm n}^{i} ,
\end{displaymath} (27)

which allows us to express the momenta $\vec{p}^{\,{\rm n}}$ and $\vec{p}^{\,{\rm p}}$ of (10) simply as
$\displaystyle \vec{p}^{\,{\rm n}}$ $\textstyle \,=\,$ $\displaystyle M^{\rm n} \vec{v}_{\rm n} , \quad {\rm with} \quad
M^{\rm n} \equ...
...} + {\cal K}^{{\rm n}{\rm p}} n_{\rm p} {\Omega_{\rm p} \over
\Omega_{\rm n}} ,$  
$\displaystyle \vec{p}^{\,{\rm p}}$ $\textstyle \,=\,$ $\displaystyle M^{\rm p} \vec{v}_{\rm n} , \quad {\rm with} \quad
M^{\rm p} \equ...
...{\cal K}^{{\rm p}{\rm p}} n_{\rm p} {\Omega_{\rm p} \over
\Omega_{\rm n}} \cdot$ (28)

Because of the axial symmetry, scalars like $M^{\rm n}$ and $M^{\rm p}$ cannot depend on the angle $\varphi$, and therefore we have

 \begin{displaymath}v_{\rm n}^{j} \nabla_{j} M^{\rm n} = 0 , {\quad{\rm and}\quad}
v_{\rm n}^{j} \nabla_{j} M^{\rm p} = 0 .
\end{displaymath} (29)

Furthermore, since we are assuming uniform rotations, cf. (26), we can deduce the identity

 \begin{displaymath}v_{\rm n}^{j} \nabla_{j} v_{{\rm n}{i}} = - {1\over 2} \nabla...
- \nabla_{i}\left({\varpi^2\over2} \Omega_{\rm n}^2\right) .
\end{displaymath} (30)

This means that the second parenthesis in the equations of motion (25) vanishes. Using (28), (29) and (30), we find
$\displaystyle p^{\rm n}_{j} \nabla_{i} v_{\rm n}^{j} + v_{\rm n}^{j} \nabla_{j} p^{\rm n}_{i}$ $\textstyle \,=\,$ $\displaystyle \hspace{1.3em} M^{\rm n} \left({1\over2}\nabla_{i} v_{\rm n}^2 + v_{\rm n}^{j} \nabla_{j}
v_{{\rm n}{i}}\right) = 0 ,$  
$\displaystyle p^{\rm p}_{j} \nabla_{i} v_{\rm p}^{j} + v_{\rm p}^{j} \nabla_{j} p^{\rm p}_{i}$ $\textstyle \,=\,$ $\displaystyle {\Omega_{\rm p}\over\Omega_{\rm n}}M^{\rm p} \left({1 \over 2}\nabla_{i} v_{\rm n}^2 + v_{\rm n}^{j}
\nabla_{j} v_{{\rm n}{i}}\right) = 0 .$  

This reduces Eqs. (25) to two first integrals of motion[*], namely
$\displaystyle \mu^{\rm n} + m^{\rm n} \Phi - {\varpi^2\over2} m^{\rm n} \Omega_{\rm n}^2$ $\textstyle \,=\,$ $\displaystyle {\cal C}^{\rm n} ,$  
$\displaystyle \mu^{\rm p} + m^{\rm p} \Phi - {\varpi^2\over2} m^{\rm p} \Omega_{\rm p}^2$ $\textstyle \,=\,$ $\displaystyle {\cal C}^{\rm p} ,$ (31)

where the constants ${\cal C}^{\rm n}$ and ${\cal C}^{\rm p}$ are in general determined by the rotation rates $\Omega_{\rm n}$ and $\Omega_{\rm p}$. These first integrals, together with Poisson's Eq. (7), an equation of state ${\cal E}(n_{\rm n},\,n_{\rm p},\,\Delta^2)$ and appropriate boundary conditions completely specify the solution. Note that these first integrals are also the Newtonian limits of the general relativistic results obtained by Andersson & Comer (2001a).

It is interesting to note that the stationary two-fluid model with different rotation rates $\Omega_{\rm n}\not=\Omega_{\rm p}$ does not allow for chemical equilibrium between the two fluids (as has already been shown in the relativistic framework by Andersson & Comer 2001a). If we assume the neutrons to be superfluid, and the entropy being carried by the "proton'' fluid, then the condition of chemical equilibrium can be found (Langlois et al. 1998; Carter et al. 2000; Prix 2001) to be

$\displaystyle \left[\mu^{\rm p} + m^{\rm p} \Phi - {1\over2}m^{\rm p} v_{\rm p}...
... - \left[ \mu^{\rm n}
+ m^{\rm n} \Phi - {1\over2} m^{\rm n} v_{\rm n}^2\right]$      
$\displaystyle = \left( m^{\rm n} v_{\rm n}^j - {2\alpha\over n_{\rm n}} \Delta^j \right) \Delta_j,$     (32)

where $\Delta_j$ is the relative velocity defined in (4). The conditions of stationarity and uniform rotations $\Omega_{\rm n}$, $\Omega_{\rm p}$ have been seen to result in the two first integrals (31), so that the requirement of chemical equilibrium (32) now reads

 \begin{displaymath}{\cal C}^{\rm p} - {\cal C}^{\rm n} = \varpi^2 (\Omega_{\rm n...
...lpha\over n_{\rm n}}(\Omega_{\rm n} - \Omega_{\rm p}) \right),
\end{displaymath} (33)

which makes it obvious that the only stationary configurations compatible with chemical equilibirum are co-rotating, i.e. $\Omega_{\rm n}=\Omega_{\rm p}$. If we allowed for "chemical ($\beta$) reactions'' ${\rm n}\rightleftharpoons{\rm p}$ (characterising the "transfusive'' type of models considered in Langlois et al. 1998), then these reactions would inevitably dissipate away the energy associated with the relative velocity $\Delta^i$ until the star reaches a chemical equilibrium state where both fluids are co-rotating. In order to have different rotation rates in a stationary state, one therefore has to exclude $\beta$-reactions altogether, or at least assume that their characteristic timescale is much longer than the timescale at which we want to consider the configuration to be approximately stationary, for example the dynamical timescale (of order ms). In practice this assumption seems physically well justified by the slowness of these reactions (Haensel 1992), which operate on timescales of months to years for mature neutron stars.

It will be convenient in the following to write our equations in dimensionless form by expressing all quantities in their "natural units.'' We choose these to be the radius R of the (non-rotating) star for lengthscales, its central density $\rho_{0}$ for densities, and $1/\sqrt{4\pi G \rho_0}$ for timescales. For simplicity of notation we will continue to use the same symbols for the dimensionless quantities, so Poisson's Eq. (7) now reads

 \begin{displaymath}\nabla^2 \Phi(r,\theta) = \rho(r,\,\theta) ,
\end{displaymath} (34)

whereas the form of the first integrals (31) is unchanged. We note that the "natural rotation rate'' $\sqrt{4\pi G \rho_0}$ is slightly greater than the maximal rotation rate, the Kepler limit, $\Omega_{\rm K}$ at which mass shedding will occur.

For realistic equations of state the mass shedding limit is well approximated by the empirical expression $\Omega_{\rm K} \approx
\sqrt{4\pi G \bar{\rho}}/3$, where $\bar{\rho}$ is the mean density. Therefore the dimensionless rotation rates $\Omega_{\rm n}$ and $\Omega_{\rm p}$ are generally smaller than 1, and for typical rotation rates $\Omega$ of most observed pulsars we will have

 \begin{displaymath}\Omega\ll 1.
\end{displaymath} (35)

Note that the difference between the "natural rotation rate'' and the mass shedding limit is further amplified by the fact that $\rho_0 \sim
5{-}10 ~\bar{\rho}$ for a typical neutron star equation of state. That being said, we must still proceed somewhat cautiously near the Kepler limit. A comparison between slow-rotation results and numerical calculations for rapidly rotating stars show that, while the slow-rotation approximation can accurately describe the fastest observed pulsars it deteriorates significantly near the Kepler limit. This is illustrated in Fig. 1, where we compare the slow-rotation results for the equatorial and polar radii of ordinary, one-fluid stars (governed by an equation of state of the form ${\cal E} \propto \rho^2$) to accurate numerical results. The latter were obtained using the LORENE code, based on pseudo-spectral techniques, which was developed by the Meudon Numerical Relativity group (Bonazzola et al. 1993; Gourgoulhon et al. 1999) and which can be used to build rapidly rotating, Newtonian and general relativistic neutron stars. We see that the slow-rotation approximation works quite well up to and including the fastest known pulsar, but starts to fail (by some $15\%$ to $20\%$) as the Kepler limit is approached. These arguments motivate the range of applicability of the "slow-rotation expansion'', which we will use in the following sections.
\end{figure} Figure 1: Comparison of one-fluid slow-rotation configurations with numerical results obtained using the LORENE code. $R_{\rm equ}$ and $R_{\rm pole}$ are the star's equatorial and polar radius repectively. The stellar model is a N=1 polytrope with mass $M=1.4~M_{\odot}$ and (static) radius R=10 km. The slow-rotation approximation leads to a significant underestimate of the rotational flattening near the Kepler limit.
Open with DEXTER

5 The slow-rotation approximation

5.1 Slow-rotation expansion of the two-fluid model

In order to approximate the solution to (31) and (34), we apply a method initially due to Chandrasekhar (1933) and Milne (1923). This method is based on the assumption that the rotation is slow enough that the configuration is only slightly oblate, in such a way that one can express it in terms of a Taylor expansion around the non-rotating spherical solution. In the following we consider such a "slowly-rotating'' two-fluid star. In the previous section we have seen that this corresponds to $\Omega_{\rm n}$ and $\Omega_{\rm p}$ being small compared to the natural rotation rate $\sqrt{4\pi G \rho_0}$. In the dimensionless form of the equations, the small parameters of the expansion are simply the rotation rates, as $\Omega_{\rm n}\ll1$ and $\Omega_{\rm p}\ll1$. We can therefore write the solution as a Taylor expansion in $\Omega_{\rm n}$ and $\Omega_{\rm p}$. In doing this we will neglect all terms beyond second order.

Under these assumptions any scalar physical quantity Q of the rotating star can be written

$\displaystyle Q(r,\theta; \Omega_{\rm n}, \Omega_{\rm p})$ $\textstyle \,=\,$ $\displaystyle \left.Q\right\vert _{\Omega_{\rm n}=\Omega_{\rm p}=0}$  
    $\displaystyle + \left.{\partial Q \over \partial \Omega_{\rm n}}\right\vert _0 ...
\left.{\partial Q \over \partial \Omega_{\rm p}}\right\vert _0 \Omega_{\rm p}$  
    $\displaystyle + {1\over2}\left.{\partial^2 Q \over \partial \Omega_{\rm n}^2}\r...
...ega_{\rm n} \partial
\Omega_{\rm p}}\right\vert _0 \Omega_{\rm n}\Omega_{\rm p}$  
    $\displaystyle + {1 \over 2} \left.{\partial^2 Q \over \partial \Omega_{\rm p}^2}\right\vert _0
\Omega_{\rm p}^2 + {\cal O}(\Omega^3) .$ (36)

Furthermore, because we are actually expanding in the "vector'' $(\Omega_{\rm n}, \Omega_{\rm p})$, it will be convenient to introduce the notation ${\underline{\Omega}}$ for the "constituent vector'' with components $\Omega_X$, where $X = {\rm n},{\rm p}$ is the constituent index. Thus we have

\begin{displaymath}{\underline{\Omega}}= \left(\Omega_{\rm n}, \, \Omega_{\rm p}...
... i.e.} \quad
\left({\underline{\Omega}}\right)_X = \Omega_X .
\end{displaymath} (37)

Terms of odd-order in ${\underline{\Omega}}$ in the expansion (36) will vanish identically because the configuration has to be invariant under a simultaneous inversion of both rotation rates, i.e. for $\Omega_{\rm n} \rightarrow -\Omega_{\rm n}$ and $\Omega_{\rm p} \rightarrow -\Omega_{\rm p}$.

We know that the static solution $Q\vert _{{\underline{\Omega}}=0}$ is purely spherical and denote it Q(0)(r), i.e.

\begin{displaymath}Q(r,\theta; {\underline{\Omega}}=0) = Q^{(0)}(r).
\end{displaymath} (38)

With these definitions, the second-order expansion (36) of Q can be written in a more compact form as

 \begin{displaymath}Q(r,\theta; {\underline{\Omega}}) = Q^{(0)}(r) + \delta Q(r,\theta) + {\cal O}({\underline{\Omega}}^4) \,
\end{displaymath} (39)


\begin{displaymath}\delta Q = \Omega_X Q^{XY} \Omega_Y,
\end{displaymath} (40)

where we automatically sum over repeated constituent indices (this convention is assumed from now on), and $Q^{XY}(r,\theta)$ is the symmetric expansion coefficient matrix

\begin{displaymath}Q^{XY}(r,\theta) \equiv \left.{\partial^2 Q \over \partial \O...
\partial \Omega_Y} \right\vert _{{\underline{\Omega}}=0} .
\end{displaymath} (41)

It has been stated several times in the literature (Monaghan & Roxburgh 1965; Hartle 1967; Smith 1975, 1976; Tassoul 1978) that the slow-rotation expansion of the density $\rho(r,\theta)$ breaks down when approaching the surface of the non-rotating star, i.e. for $r\rightarrow 1$, because there $\rho^{(0)}(1)=0$. However, it is evident from (39) that the validity of this expansion is completely independent of the values of Q(0)(r). The only condition for its validity is that the neglected terms ${\cal O}({\underline{\Omega}}^4)$ are effectively small compared to the ${\cal O}({\underline{\Omega}}^2)$ terms included in the analysis.

We now expand all physical quantities of (31) and (34) up to second-order in ${\underline{\Omega}}$:

    $\displaystyle \Phi(r,\theta) = \Phi^{(0)}+ \delta\Phi , \qquad
\delta\Phi = \Omega_X \,\Phi^{XY}\, \Omega_Y,$ (42)
    $\displaystyle \rho(r,\,\theta) = \rho^{(0)}+ \;\delta \rho , \qquad
\;\delta\rho = \Omega_X \,\rho^{XY}\, \Omega_Y,$ (43)
    $\displaystyle n_A(r,\theta) = n_A^{(0)}\;+ \delta n_A, \qquad
\delta n_A= \Omega_X \,n_A^{XY}\,\Omega_Y,$ (44)
    $\displaystyle \mu^A(r,\theta) = \mu^{A{(0)}} + \delta\mu^A, \qquad
\delta\mu^A= \Omega_X \,\mu^{A,\,XY}\,\Omega_Y,$ (45)
    $\displaystyle {\cal C}^A= {\cal C}^{A{(0)}} + \;\delta{\cal C}^A,\qquad
\delta{\cal C}^A= \Omega_X \,{\cal C}^{A,\,XY} \,\Omega_Y.$ (46)

The total mass density $\rho$ is, of course, linked to the individual particle densities by (3), i.e. $\rho = m^A\,n_A$. Therefore we have the relations

\begin{displaymath}\rho^{(0)}= m^A\, n_A^{(0)},{\quad{\rm and}\quad}
\rho^{XY} = m^A\, n_A^{XY}.
\end{displaymath} (47)

Inserting the various expansions into (31) and (34), the zeroth-order equations, which determine the non-rotating configuration, are found to be
$\displaystyle \mu^{A{(0)}} + m^A\Phi^{(0)}$ $\textstyle \,=\,$ $\displaystyle {\cal C}^{A{(0)}} ,$ (48)
$\displaystyle \nabla^2 \Phi^{(0)}(r)$ $\textstyle \,=\,$ mAnA(0)(r) , (49)

while the second-order coefficients of (42)-(46) satisfy
$\displaystyle \mu^{A,\,XY} + m^A\Phi^{XY} - {\varpi^2\over2}
\delta^{A,\,XY}$ $\textstyle \,=\,$ $\displaystyle {\cal C}^{A,\,XY} ,$ (50)
$\displaystyle \nabla^2 \Phi^{XY}$ $\textstyle \,=\,$ mAnAXY . (51)

In Eq. (50) we have introduced the constant matrices $\delta^{A,\,XY}$, which are defined as

 \begin{displaymath}\delta^{{\rm n},\,XY} \equiv \left(\begin{array}{rr} m^{\rm n...
0 & ~~0 \\
0 & ~~m^{\rm p} \end{array}\right) ,
\end{displaymath} (52)

These allow us to write the second-order terms $m^{\rm n}\Omega_{\rm n}^2$ and $m^{\rm p}\Omega_{\rm p}^2$ in (31) as

\begin{displaymath}m^{\rm n}\Omega_{\rm n}^2 = \Omega_X \,\delta^{{\rm n},\,XY} \, \Omega_Y\,
\end{displaymath} (53)


\begin{displaymath}m^{\rm p}\Omega_{\rm p}^2 = \Omega_X \,\delta^{{\rm p},\,XY} \, \Omega_Y.
\end{displaymath} (54)

5.2 Reduction to a single equation

We can find an algebraic relation between the chemical potential coefficients $\mu^{A,\,XY}$ and the density coefficients nAXY, which will allow us to considerably simplify the problem. This relation is found by expanding the chemical potential  $\mu^A(n_{\rm n}, n_{\rm p},
\Delta^2)$ defined in (5) in its arguments up to second-order in ${\underline{\Omega}}$. Using (44) we find

$\displaystyle \mu^A(n_{\rm n},\, n_{\rm p},\, \Delta^2)$ $\textstyle \,=\,$ $\displaystyle \mu^A(n_{\rm n}^{(0)}, n_{\rm p}^{(0)}, 0)$  
    $\displaystyle + \Omega_X \, \left(\left.{\partial \mu^A\over \partial
n_Z}\right\vert _0 n_Z^{XY} \right.$  
    $\displaystyle + \left.\left.{\partial \mu^A\over \partial
\Delta^2}\right\vert _0 \varpi^2 \Delta^{XY}\right) \Omega_Y+
{\cal O}({\underline{\Omega}}^4) ,$ (55)

where we have defined the constant matrix $\Delta^{XY}$ as

 \begin{displaymath}\Delta^{XY} \equiv \left(\begin{array}{rr} 1 & ~-1 \\
- 1 & ~1 \end{array} \right).
\end{displaymath} (56)

This definition allows us to write the relative velocity squared, namely $\Delta^2=\varpi^2 \left(\Omega_{\rm n} - \Omega_{\rm p}\right)^2$, in the form

\begin{displaymath}\Delta^2 = \varpi^2 \,\Omega_X\,\Delta^{XY} \, \Omega_Y.
\end{displaymath} (57)

The partial derivatives in (55) are evaluated for the non-rotating configuration, and thus they are assumed to be known functions, depending only on the static solution for a given equation of state ${\cal E}(n_{\rm n},\,n_{\rm p},\,\Delta^2)$. Thus we define the "density structure function'' $\S_{XY}(r)$ as

 \begin{displaymath}\S_{XY}(r)\equiv \left(\left.{\partial \mu^X \over \partial
\partial n_X \partial n_Y}\right\vert _0 \right)^{-1} ,
\end{displaymath} (58)

which is symmetric (and we assume invertible with inverse $\left(\S^{-1}\right)^{XY}$), and the "entrainment structure function'' $\beta^X(r)$ as

 \begin{displaymath}\beta^X(r) \equiv \left.{\partial \mu^X \over \partial
...= \left.{\partial \alpha\over \partial
n_X } \right\vert _0,
\end{displaymath} (59)

where we used the definition (5) of the chemical potentials $\mu^X$ and the entrainment function $\alpha$. Comparing (55) to the expansion (45) we can identify
    $\displaystyle \mu^{A{(0)}} = \mu^A(n_{\rm n}^{(0)}, n_{\rm p}^{(0)}, 0),$ (60)
    $\displaystyle \mu^{A,\,XY} = \left(\S^{-1}\right)^{AB} \,n_B^{XY} +
\varpi^2 \beta^A\Delta^{XY},$ (61)

and we have thus arrived at an algebraic relation between the $\mu^{A,\,XY}$ and the nBXY.

Inserting this relation into the first integral (50), we obtain an explicit expression for the density coefficients nAXY in terms of the gravitational potential coefficients $\Phi^{XY}$. Specifically, we have

nAXY $\textstyle \,=\,$ $\displaystyle \S_{AB} \Big[ C^{B,\,XY} + {\varpi^2 \over 2}
\left(\delta^{B,\,XY} - 2 \beta^B\Delta^{XY}\right)$  
    $\displaystyle -m^B\Phi^{XY} \Big].$ (62)

Introducing the "derived'' background functions
    $\displaystyle E_A^{XY}(r) \equiv {1 \over 3} \S_{AB}(r) \left(
\delta^{B,\,XY} - 2 \beta^B(r) \Delta^{XY}\right) ,$  
    $\displaystyle k_A(r) \equiv \S_{AB}(r) m^B,$ (63)

Eq. (62) can be written

 \begin{displaymath}n_A^{XY} = \S_{AB}(r) {\cal C}^{B,\,XY} + {3 \varpi^2 \over 2}
E_A^{XY}(r) - k_A(r) \,\Phi^{XY} ,
\end{displaymath} (64)

which is a reformulation of the first integrals of motion (50) in terms of the second-order coefficients.

The total mass density coefficient $\rho^{XY}$ can now be written

$\displaystyle \rho^{XY}$ $\textstyle \,=\,$ mAnAXY  
  $\textstyle \,=\,$ $\displaystyle m^A\S_{AB}(r) {\cal C}^{B,\,XY} + {3\varpi^2 \over 2}
E^{XY}(r) - k(r) \Phi^{XY} ,$ (65)

where we have introduced the further abbreviations
$\displaystyle E^{XY}(r)\,$ $\textstyle \equiv$ $\displaystyle \, m^AE_A^{XY}(r)$  
$\displaystyle \,$ = $\displaystyle \, {1\over3} m^A\S_{AB}(r) \left(\delta^{B,\,XY}
- 2\beta^B(r)\,\Delta^{XY} \right) \,$ (66)

\begin{displaymath}k(r) \equiv m^Ak_A(r) = m^A\S_{AB}(r) \,m^B, \end{displaymath}

which are "known'' functions of the non-rotating star's configuration, determined in terms of the "structure functions'' $\S_{XY}(r)$ and $\beta^X(r)$, together with the constants  $\delta^{A,\,XY}$ and $\Delta^{XY}$ and the constants of integration ${\cal C}^{A,\,XY}$. These constants are determined by the boundary conditions. Inserting (65) into Poisson's Eq. (50) reduces the problem to a single partial differential equation for the coefficients $\Phi^{XY}(r,\theta)$, namely

 \begin{displaymath}\nabla^2 \Phi^{XY} + k(r) \Phi^{XY}={3\varpi^2\over2}
E^{XY}(r) + m^A\S_{AB}(r) {\cal C}^{B,\,XY}.
\end{displaymath} (67)

Given the non-rotating background, a solution to this equation fully specifies a slowly-rotating, two-fluid configuration.

We note that a necessary condition for this method to work, i.e. for (67) to be well-defined, is that the structure function $\S_{XY}(r)$ is regular everywhere inside the star. However, it is well-known that in the case of a single fluid the standard polytropic equations of state, i.e.  ${\cal E}=\kappa n^{1+1/N}$, leads either to a vanishing (for N>1) or an infinite density gradient (for N<1) at the surface. This behaviour could make (58) singular at the surface. In order to clarify this point, we can express the definition (58) equivalently as

 \begin{displaymath}{n_X^{{(0)}}}'(r) = \S_{XY} \, {\mu^{Y{(0)}}}'(r),
\end{displaymath} (68)

where a prime denotes the radial derivative ${\rm d}/{\rm d}r$. It is seen from (48) that the chemical potentials $\mu^{X{(0)}}(r)$ and their derivatives are regular everywhere, even at the surface. Therefore we can conclude that the structure function $\S_{XY}$ behaves like the density gradients nX(0)' at the surface. The slow rotation expansion used in the present work is therefore well-defined for stars with finite or zero density gradients at the surface, but is not applicable to stars with an infinite density gradient at the surface, similar to the N<1 case for the polytropes, as the structure function SXY then diverges at the surface.

Inserting (48) into (68), we obtain expressions for the derived structure functions kA(r)(63) and k(r), which indicate their physical meaning. Using (63) we have

 \begin{displaymath}{n_A^{(0)}}' = - k_A\, {\Phi^{(0)}}' \quad
\Longrightarrow\quad k_A= - {{\rm d}n_A^{(0)}\over {\rm d}\Phi^{(0)}},
\end{displaymath} (69)

and from the definition (66) it follows that

 \begin{displaymath}{\rho^{(0)}}' = - k \,{\Phi^{(0)}}' \quad
\Longrightarrow\quad k = - {{\rm d}\rho^{(0)}\over {\rm d}\Phi^{(0)}} \cdot
\end{displaymath} (70)

From these relations we see that the "structure'' functions kA(r)reflect the change of the individual number densities with respect to the gravitational potential, while the function k(r) reflects the corresponding change in the total mass density. These relations proved very useful in deriving the results discussed in Sect. 7.

5.3 Separation of variables

So far, the slow-rotation approximation has allowed us to reduce the initial problem (31) and (34) to a single Eq. (67) for a single unknown function $\Phi^{XY}(r,\theta)$. We now separate the variables r and $\theta$ by expanding $\Phi^{XY}$ in the orthogonal basis of Legendre polynomials $P_l(\theta)$. These are the eigenfunctions of the angular part of the Laplace operator, $\nabla^2_\theta$, satisfying

\begin{displaymath}\nabla^2 P_l(\theta) = -{l (l + 1) \over r^2} P_l(\theta) .
\end{displaymath} (71)

Therefore, writing

 \begin{displaymath}\Phi^{XY}(r,\,\theta) = \sum_{l=0}^\infty \Phi^{XY}_l(r)
\end{displaymath} (72)

leads to

 \begin{displaymath}\nabla^2 \Phi^{XY} = \sum_{l=0}^\infty P_l(\theta) \, {\cal D}_l
\Phi^{XY}_l(r) ,
\end{displaymath} (73)

where the differential operator ${\cal D}_l$ is defined as

\begin{displaymath}{\cal D}_l \equiv {{\rm d}^2 \over {\rm d}r^2} + {2 \over r} {{\rm d} \over {\rm d}r} - {l (l + 1)
\over r^2},
\end{displaymath} (74)

which represents the radial part of the Laplacian with an additional "angular momentum'' part proportional to l(l + 1).

Using (73) together with the identity

\begin{displaymath}{3 \varpi^2 \over 2} = r^2 \left[1 - P_2(\theta)\right] ,
\end{displaymath} (75)

we arrive at a series of ordinary differential equations for the coefficients $\Phi^{XY}_l(r)$ in (72). We have
    $\displaystyle {\cal D}_0 \Phi^{XY}_0 + k\,\Phi^{XY}_0 = r^2\,E^{XY}
+ m^A\S_{AB}{\cal C}^{B,\,XY},$ (76)
    $\displaystyle {\cal D}_2 \Phi^{XY}_2 + k\,\Phi^{XY}_2 = - r^2
\, E^{XY} ,$ (77)
    $\displaystyle {\cal D}_l \Phi^{XY}_l + k \, \Phi^{XY}_l = 0 ,
\quad {\rm for} \quad l \ge 4$ (78)

where only terms of even l will be non-vanishing because of the expected equatorial symmetry of the solution. Furthermore, we will see in the following section that the boundary conditions are homogeneous, and therefore the homogeneous Eqs. (78) have only trivial solutions, i.e.

\begin{displaymath}l \ge 4: \quad \Phi^{XY}_l(r) = 0 .
\end{displaymath} (79)

The only non-vanishing contributions to $\Phi^{XY}(r,\,\theta)$ therefore consist of the two functions $\Phi^{XY}_0(r)$ and $\Phi^{XY}_2(r)$, and the problem has been reduced to that of solving the two differential matrix Eqs. (76) and (77). This is, of course, in complete analogy with the standard result of the single fluid slow-rotation approximation.

Once the solutions $\Phi^{XY}_l$ to the differential Eqs. (76) and (77) have been obtained for a given equation of state (subject to the boundary conditions to be discussed in the next section), the individual density coefficients nAXY are determined from (64).

6 Boundary conditions and integration constants

6.1 Boundary conditions

As Eqs. (76)-(78) are linear, second-order differential equations, two boundary conditions are necessary for each equation. In the case of the l = 0 Eq. (76) we need two further conditions in order to fix the two constants of integration ${\cal C}^{{\rm n},\,XY}$ and ${\cal C}^{{\rm p},\,XY}$. First, we require that the solution must be regular at the centre of the star. At r=0 the differential operator ${\cal D}_l$ is singular, which leads to the first boundary conditions

l = 0 : $\displaystyle \quad {\Phi^{XY}_0}'(0) = 0 ,$  
$\displaystyle l \ge 2$ : $\displaystyle \quad \Phi^{XY}_l(0) = 0, \, {\quad{\rm and}\quad}{\Phi^{XY}_l}'(0) =
0,$ (80)

where the primes represent the radial derivative ${\rm d}/{\rm d}r$. Secondly, we require continuity of the gravitational potential and its derivative across the star's surface. That is, we have to match the interior solution $\Phi(r,\,\theta)$ to the exterior gravitational potential $\psi(r,\,\theta)$ at the surface of the rotating star $R(\theta)$:

 \begin{displaymath}\left.\Phi\right\vert _{R(\theta)} = \left.\psi\right\vert _{...
...ght\vert _{R(\theta)} = \left.\psi'\right\vert _{R(\theta)} .
\end{displaymath} (81)

Here, the radial derivative deviates from the normal derivative with respect to the actual surface $R(\theta)$ only by terms ${\cal O}({\underline{\Omega}}^4)$, which have been neglected.

We write the slow-rotation expansion of the exterior potential $\psi$ in the usual way as

\begin{displaymath}\psi(r,\,\theta) = \psi^{(0)}(r) + \delta \psi(r,\,\theta)+
{\cal O}({\underline{\Omega}}^4) ,
\end{displaymath} (82)


\begin{displaymath}\delta\psi = \Omega_X \, \psi^{XY}\Omega_Y,
\end{displaymath} (83)

and expand the second-order coefficient in Legendre polynomials:

 \begin{displaymath}\psi^{XY}(r,\,\theta) = \sum_{l=0}^\infty \psi^{XY}_l(r)
P_l(\theta) .
\end{displaymath} (84)

Because $\psi$ is the solution to the Laplace equation $\nabla^2 \psi
= 0$ (normalised by $\lim_{r\rightarrow\infty}\psi=0$), we know the radial eigenfunctions for all l, namely

 \begin{displaymath}\psi^{XY}_l(r) = {\kappa^{XY}_l \over r^{l+1}} ,
\end{displaymath} (85)

where the $\kappa^{XY}_l$ are constant coefficients. The continuity condition (81) for the static solution implies

 \begin{displaymath}\Phi^{(0)}(1) = \psi^{(0)}(1) \, {\quad{\rm and}\quad}
{\Phi^{(0)}}'(1) = {\psi^{(0)}}'(1) ,
\end{displaymath} (86)

and because the static potential $\Phi^{(0)}$ at the surface also satisfies the Laplace equation, i.e. $\left.\nabla^2
\Phi^{(0)}\right\vert _{r=1}=0$, the second derivatives must also be identical:

 \begin{displaymath}{\Phi^{(0)}}''(1) = {\psi^{(0)}}''(1) .
\end{displaymath} (87)

Up to second-order in ${\underline{\Omega}}$ the surface of the rotating star is

 \begin{displaymath}R(\theta) = 1 + \delta R(\theta) + {\cal O}({\underline{\Omega}}^4) ,
\end{displaymath} (88)

where we have used R(0)=1. This allows us to write the second-order expansion of the internal and external potentials at the matching surface $R(\theta)$ as
$\displaystyle \left.\Phi\right\vert _{R(\theta)}$ $\textstyle \,=\,$ $\displaystyle \Phi^{(0)}(1) + {\Phi^{(0)}}'(1) \,
\delta R + \delta \Phi(1,\theta) ,$ (89)
$\displaystyle \left.\psi\right\vert _{R(\theta)}$ $\textstyle \,=\,$ $\displaystyle \psi^{(0)}(1) + {\psi^{(0)}}'(1) \,
\delta R + \delta \psi(1,\theta) ,$ (90)

and of the radial derivatives:
$\displaystyle \left.\Phi'\right\vert _{R(\theta)}$ $\textstyle \,=\,$ $\displaystyle {\Phi^{(0)}}'(1) +
{\Phi^{(0)}}''(1) \,\delta R + \delta \Phi'(1,\theta) ,$ (91)
$\displaystyle \left.\psi'\right\vert _{R(\theta)}$ $\textstyle \,=\,$ $\displaystyle {\psi^{(0)}}'(1) +
{\psi^{(0)}}''(1) \,\delta R + \delta \psi'(1,\theta) .$ (92)

Therefore the matching conditions (81) together with (86) and (87) simply reduce to

 \begin{displaymath}\delta \Phi(1,\,\theta) = \delta \Psi(1,\,\theta) , {\quad{\r...
\delta \Phi'(1,\,\theta) = \delta \Psi'(1,\,\theta) ,
\end{displaymath} (93)

which shows that the resulting boundary condition is - since there are no terms involving $\delta R$ - independent of the actual second-order surface of matching. It is thus not necessary to determine the perturbed surface of the star explicitly (by the condition $\left.\rho\right\vert _{R(\theta)}=0$). Using (84), (85) and (72) for the second-order terms $\delta\psi$ and $\delta\Phi$ in (93) directly results in the second boundary condition

 \begin{displaymath}{\Phi^{XY}_l}' (1) + (l+1) \,\Phi^{XY}_l(1) = 0 , \quad
{\rm for} \quad l \ge 0 .
\end{displaymath} (94)

6.2 Fixing the integration constants

As stated earlier, the l=0 equation requires two additional conditions in order to fix the constants of integration ${\cal C}^{{\rm n},\,XY}$ and ${\cal C}^{{\rm p},\,XY}$. These conditions are crucial as they determine which type of rotating star sequence (as a function of the rotation rates $\Omega_{\rm n}$ and $\Omega_{\rm p}$) the solution will describe. Among the many possible sequences, two seem particularly useful and will be described here. The first is the fixed central density sequence (FCD), which is probably the most straightforward and therefore most commonly considered choice.

The FCD sequence is characterized by the simple condition (for $A=
{\rm n},\,{\rm p}$)

\begin{displaymath}\left.n_A\right\vert _{r=0} = \left. n_A^{(0)}\right\vert _{r=0}
\quad \Longrightarrow \quad n^{XY}_{A,\,0}(0) = 0 .
\end{displaymath} (95)

This condition obviously implies that the respective total masses will change as functions of the rotation rates. This sequence therefore does not describe the same physical star at different rotation rates, but it has the advantage of leading to a very simple condition. From (62) we find that we should require

 \begin{displaymath}{\cal C}^{A,\,XY} = m^A\Phi^{XY}_0(0).
\end{displaymath} (96)

The FCD sequence is, however, not particularly relevant from a physical point of view. For example, if we consider an isolated neutron star that is spinning down due to a magnetic torque we would expect the central density to increase as the star becomes less oblate. The physically most interesting sequence of rotating stars thus corresponds to requiring fixed masses (FM). In other words, one has to impose the condition of constant respective total masses (or equivalently, total particle numbers) for each of the two fluids if the sequence is to describe the same physical star at different rotation rates[*]. This means that (for $A=
{\rm n},\,{\rm p}$)

\begin{displaymath}\int n_A(r,\,\theta) \, {\rm d}V = \int n_A^{(0)}(r) \, {\rm d}V .
\end{displaymath} (97)

To second-order in the rotation rates, this condition reduces to

 \begin{displaymath}\int_0^1 r^2 n_{A,\,0}^{XY}(r) \, {\rm d}r = 0 .
\end{displaymath} (98)

Inserting the explicit expression (64) for the $n_{A,\,0}^{XY}$ component yields the integral condition (for $A, B
= {\rm n},\,{\rm p}$)

\begin{displaymath}\int_0^1 r^2 \left[ \S_{AB} {\cal C}^{B,\,XY} + r^2 E_A^{XY} - k_A
\Phi_0^{XY}(r) \right] \, {\rm d}r = 0 ,
\end{displaymath} (99)

which allows us to determine the constants ${\cal C}^{A,\,XY}$ for the fixed mass sequence. However, in practice one only needs to evaluate this integral for one of the two fluids, because the second fixed mass condition can be replaced by the equivalent but simpler requirement of fixed total mass, i.e. the integral over $\rho = m^X\, n_X$, which analogously to (98) leads to

\begin{displaymath}\int_0^1 r^2 m^An_{A,\,0}^{XY} \, {\rm d}r= 0 .
\end{displaymath} (100)

The l=0 component of Poisson's Eq. (51), i.e.

\begin{displaymath}{1 \over r^2} \left(r^2 {\Phi_0^{XY}}'(r)\right)' = m^A
n_{A,\,0}^{XY} ,
\end{displaymath} (101)

allows us to reduce this condition to

\begin{displaymath}{\Phi_0^{XY}}'(1) = 0 , \quad \Longleftrightarrow \quad
\Phi_0^{XY}(1) = 0 ,
\end{displaymath} (102)

where the second (equivalent) expression has been obtained from the second boundary condition (94).

7 Ellipticities, moments of inertia, and Kepler rotation rate

Before proceeding to discuss results for a specific model equation of state it is useful to digress on what kind of information we want to extract from the calculation. There are obviously many alternative ways of describing rotating configurations. We will focus our attention on the ellipticities, the moments of inertia and the Kepler rotation rate. The first and second are interesting because they highlight how rotation affects the shape and mass-distribution of the star, while the last describes the limit of rotation at which mass-shedding at the equator sets in. We can explicitly express the respective fluid surfaces  $R_A(\theta)$ to second-order in the rotation rates. Starting from the obvious definition

\begin{displaymath}\left. n_A(r,\,\theta)\right\vert _{r=R_A(\theta)} = 0 ,
\end{displaymath} (103)

and writing the second-order surfaces as

\begin{displaymath}R_A(\theta) = R_A^{(0)}+ \delta R_A(\theta) + {\cal O}({\underline{\Omega}}^4) ,
\end{displaymath} (104)

we can express the second-order correction terms as

 \begin{displaymath}\delta R_A(\theta) = \left.- {\delta n_A(r, \theta) \over
{n_A^{(0)}}'(r)} \right\vert _{R_A^{(0)}} ,
\end{displaymath} (105)

where primes again denote radial derivatives. We define the ellipticity of a surface in terms of the radii of the equator $R_{\rm equ}$ and the pole $R_{\rm pole}$ as

\begin{displaymath}e \equiv {R_{\rm equ} - R_{\rm pole} \over R_{\rm equ}} \cdot
\end{displaymath} (106)

To second-order, this means that

 \begin{displaymath}e_A= {3 \over 2 R_A^{(0)}} \left\vert {\delta n_{A,2} \over
{n_A^{(0)}}'}\right\vert _{R_A^{(0)}}.
\end{displaymath} (107)

Here and in the following we will assume that the two surfaces coincide in the static case, i.e. $R_{\rm n}^{(0)}= R_{\rm p}^{(0)}= 1$, in order to simplify the discussion.

The respective moments of inertia IA to second-order can be written as

\begin{displaymath}I_A= I^{(0)}_A+ \delta I_A+ {\cal O}({\underline{\Omega}}^4) ,
\end{displaymath} (108)

where the second order correction term can be seen to be

\begin{displaymath}\delta I_A= m \int_{V^{(0)}_A} \varpi^2 \delta n_A(r,\theta) \,
{\rm d}V,
\end{displaymath} (109)

which reduces to the explicit integral

 \begin{displaymath}\delta I_A= {8 \pi m \over 3} \int_0^1 \left(\delta n_{A,0}(r) -
{1 \over 5} \delta n_{A,2}(r)\right) r^4 {\rm d}r.
\end{displaymath} (110)

The Kepler rotation rate $\Omega_{\rm K}$ is reached when one of the two fluids spins at the rate of a test particle in Keplerian orbit around the equator of the star, i.e. at $\theta = \pi /2$ and $r=R_A(\pi/2)$, where we assume that Arepresents the "outer'' fluid at the equator. This means that

\begin{displaymath}\left(\partial_r \Phi(r,\pi/2) - r \Omega_{\rm K}^2 \right)_{r=R_A(\pi/2)} = 0
\end{displaymath} (111)

Due to the rotation-induced change in the shape of the star, the Kepler rate also has a ${\cal O}({\underline{\Omega}}^2)$ correction, and so we write

 \begin{displaymath}\Omega_{{\rm K},{\it A}}^2 = \Omega_{(0)}^2 + \delta\Omega_{{\rm K},{\it A}}^2 + {\cal O}({\underline{\Omega}}^4) ,
\end{displaymath} (112)

where $\Omega_{(0)}$ is the Kepler rate around a spherical star, i.e.

\begin{displaymath}\Omega_{(0)}^2 \equiv {\Phi^{(0)}}'(1) = {4\over3}\pi G\bar{\rho} .
\end{displaymath} (113)

Using the second-order expansions for $R_A(\theta)$ and $\Phi$ this becomes

 \begin{displaymath}\delta\Omega_{{\rm K},{\it A}}^2 = - 3 {\Phi^{(0)}}'(1) \,\de...
...) +
\delta \Phi'(1,\pi/2) + {\cal O}({\underline{\Omega}}^4),
\end{displaymath} (114)

and further using (105) and (69), we can reformulate this as

 \begin{displaymath}\delta \Omega_{{\rm K},{\it A}}^2 = - {3 \over k_A(1)} \delta n_A(1,\pi/2) + \delta
\end{displaymath} (115)

The Kepler rotation rate $\Omega_{\rm K}$ determines the maximal rotation rate of the respective fluids before mass-shedding will occur at the equator, i.e. if A denotes the "outer'' fluid at the equator, then

\begin{displaymath}\Omega_A\le \Omega_{{\rm K},{\it A}}.
\end{displaymath} (116)

Thus (115) is a quadratic expression that determines the respective maximal rotation rates of the two fluids. Given the results shown in Fig. 1 we would expect the above equations to determine the Kepler limit to within 10-20%.

8 An analytic solution

We have now completed the description of our slow-rotation formalism for two-fluid systems. Given any suitable equation of state (which must provide all the relevant parameters, e.g. pertaining to entrainment) Eqs. (76) and (77) can be solved to produce a rotating configuration. For a typical equation of state (EOS) the calculation will obviously require numerics. However, although we have written a code that solves this problem given a general EOS, we will not discuss such numerical results here. There are two reasons why we do not consider a tabulated realistic EOS in this paper. Firstly, it is well known that the use of realistic EOS in the Newtonian context is somewhat dubious. The main reason is that the masses and radii of stellar models determined for a specified central density differ greatly in Newtonian theory and General Relativity. Secondly, we need the tabulated EOS data to include also the parameters needed for our entrainment model. Unfortunately, the required parameters, like the effective proton mass, are often not included among the published data, which complicates the use of most of the current EOS. Instead we focus our attention on a somewhat surprising fact: it is possible to find an analytic solution to our equations, including the general case where the two fluids are rotating at different rates, for a particular class of EOS. This model EOS corresponds to constant structure functions $\S_{XY}$ and $\beta^X$ (and therefore also constant coefficients k, EXY...). As we will argue below, the resultant model is reasonably realistic and we expect it to prove useful in future studies of the dynamics of superfluid neutron stars. Before deriving the analytic solution, we shall analyze this particular class of EOS and assess its physical relevance.

8.1 The "analytic'' equation of state

Within the present slow-rotation approximation any equation of state can (quite generally) be expressed as

$\displaystyle {\cal E}(n_{\rm n},\, n_{\rm p},\, \Delta^2)$ $\textstyle \,=\,$ $\displaystyle {\cal E}^{(0)}(n_{\rm n},\, n_{\rm p}) +
\alpha^{(0)}(n_{\rm n},\, n_{\rm p})\,\Delta^2$  
    $\displaystyle + {\cal O}({\underline{\Omega}}^4) ,$ (117)

since the relative velocity $\Delta$ is ${\cal O}({\underline{\Omega}})$. We recall that the derived "structure functions'' EXY, k etc. of (63) and (66) are determined in terms of the two basic structure functions, the "density structure'' $\S_{XY}(r)$, cf. (58), and the "entrainment structure'' $\beta^X(r)$, cf. (59). In the slow-rotation expansion (117), they become

 \begin{displaymath}\S_{XY} = \left({\partial^2 {\cal E}^{(0)}\over \partial n_X
... \left({\partial \alpha^{(0)}\over \partial n_X}\right)_0\cdot
\end{displaymath} (118)

The condition of constant structure functions determines the following class of equations of state[*]:

 \begin{displaymath}{\cal E}= {1\over2}n_X\,A^{XY} \, n_Y+ \Delta^2 B^X n_X ,
\end{displaymath} (119)

with constant coefficients AXY (symmetric) and BX. Therefore (118) leads to the identifications

\begin{displaymath}\S_{XY} = \left(A^{-1}\right)_{XY} , {\quad{\rm and}\quad}
\beta^X = B^X .
\end{displaymath} (120)

This class of EOS is the natural generalization to the two-fluid case of the "analytic polytrope'' for single fluids ( ${\cal E} \propto \rho^2$) for which the rotational corrections can be found analytically (Chandrasekhar & Lebovitz 1962). This model EOS also contains the simple case of the sum of two such polytropes without entrainment, i.e. ${\cal E}\sim A^{{\rm n}{\rm n}} n_{\rm n}^2 + A^{{\rm p}{\rm p}} n_{\rm p}^2$, for which the corresponding analytical slow-rotation solution has already been found by Prix (1999).

The major novelty of (119) concerns the inclusion of entrainment, and it is obviously important to evaluate how well this model approximates "realistic entrainment''. From the relation (19) between entrainment and effective masses we see that the equation of state (119) is simply characterized by a constant effective mass $m^{\rm p*}$. We noted earlier that nuclear physics calculations of the effective proton mass imply that this is not an unreasonable approximation, even in the extreme context of neutron stars. This suggests that, not only does (119) have the attractive feature of leading to an analytical solution for slowly rotating configurations, it also appears to provide a reasonable approximation which should correctly reflect the main qualitative features of neutron star matter (in particular the entrainment).

Another novel feature of (119) is the off-diagonal term  $A^{{\rm n}{\rm p}}$, which leads to the term $S_{{\rm n}{\rm p}}$. It will turn out that the combination

 \begin{displaymath}\sigma \equiv {\S_{{\rm n}{\rm p}} \over \S_{{\rm p}{\rm p}}}...
...ial n_{\rm p} \over \partial
\mu_{\rm n}/\partial n_{\rm n}}
\end{displaymath} (121)

plays an especially important role in determining the structure of our neutron star model. For a reasonable EOS, the matrix $\S_{XY}$ has to be positive definite, which implies the constraint

 \begin{displaymath}\S_{{\rm n}{\rm p}}^2 < \S_{{\rm n}{\rm n}} \S_{{\rm p}{\rm p}} .
\end{displaymath} (122)

This then determines the range of acceptable values for $\sigma$. For realistic situations the proton fraction will be $x_{\rm p} < 0.5$, in which case one finds $\S_{{\rm p}{\rm p}} < \S_{{\rm n}{\rm n}}$. Therefore we can be sure that in this case (122) is satisfied if $\sigma\in[-1,1]$. We have checked the reasonableness of this range by determining $\sigma$ for a more realistic EOS due to Prakash et al. (1988). For the range of densities and proton fractions to be considered here, we find that $\sigma$ typically decreases monotonically from values slightly larger than one at normal nuclear matter density $\rho_{\rm nuc}$, to values around - 0.4 at about $10
\rho_{\rm nuc}$, with the zero occurring at around $2 \rho_{\rm nuc}$. Thus, it seems that restricting a constant $\sigma$ to the range $\sigma\in[-1,1]$ is reasonable. Finally, we note that $\sigma$ is intimately related to the so-called "symmetry energy'' term in the Prakash et al. (1988) EOS, which is designed to vanish whenever there are equal numbers of neutrons and protons.

8.2 The static solution

We begin by discussing the static configuration, which is the solution to (48) and (49). The chemical potentials (5) for the "analytic'' equation of state (119) are given by

\begin{displaymath}\mu^{X{(0)}} = A^{XY} n_Y^{(0)},
\end{displaymath} (123)

which vanish at the surface of the star, and therefore the static constants of integration ${\cal C}^{X{(0)}}$ are determined from (48) at the surface r=1, i.e.

\begin{displaymath}{\cal C}^{X{(0)}} = m^X\, \Phi^{(0)}(1)\cdot
\end{displaymath} (124)

Equation (48) can now be expressed in the form

 \begin{displaymath}n_X^{(0)}(r) = k_X \left[\Phi^{(0)}(1) - \Phi^{(0)}(r)\right] ,
\end{displaymath} (125)

and therefore the ratio of the respective densities in the static configuration with this equation of state is seen to be a constant, i.e.

 \begin{displaymath}{n_{\rm p}^{(0)}(r)\over n_{\rm n}^{(0)}(r)} = {k_{\rm p} \over k_{\rm n}} \cdot
\end{displaymath} (126)

This implies that the two fluids share a common outer surface at r=1. Inserting the densities (125) into the second static Eq. (49) results in the Lane-Emden equation for the total mass density $\rho^{(0)}$:

\begin{displaymath}\nabla^2 \rho^{(0)}(r) + k\, \rho^{(0)}(r) = 0,
\end{displaymath} (127)

which has the well known solution

 \begin{displaymath}\rho^{(0)}(r) = \rho_0{\sin(r\sqrt{k})\over r\sqrt{k} } ,
\quad{\rm where}\quad \rho_0 = 1 .
\end{displaymath} (128)

The requirement that the density vanishes at the surface r=R (where here R=1) determines the constant k, which is related via (66) to the symmetric structure matrix $\S_{XY}$, namely

 \begin{displaymath}k = \S_{XY} m^X m^Y= \pi^2 .
\end{displaymath} (129)

Integrating the density (128), we obtain the following relation between the total mass M, radius R and the central density $\rho_{0}$ for this EOS,

\begin{displaymath}M = {4 \over \pi} \rho_0 R^3 .
\end{displaymath} (130)

8.3 Analytic solution for the slow-rotation coefficients

Given the "analytic equation of state'' (119), the structure functions are constant, and so are the coefficients in the differential Eqs. (76) and (77). Therefore one can write their general regular solutions as

$\displaystyle \Phi_0^{XY}(r)$ $\textstyle \,=\,$ $\displaystyle {\cal A}^{XY}_0\, {J_{1/2}(r\sqrt{k}) \over
\sqrt{r}} + {E^{XY}\over k}\left( r^2 - {6\over k}\right)$  
    $\displaystyle +{m^A\S_{AB}{\cal C}^{B,\,XY}\over k} ,$  
$\displaystyle \Phi_2^{XY}(r)$ $\textstyle \,=\,$ $\displaystyle {\cal A}^{XY}_2\,{J_{5/2}(r \sqrt{k}) \over
\sqrt{r}} - {E^{XY}\over k} r^2,$ (131)

where J1/2(x) and J5/2(x) are the standard Bessel functions. From now on we will set $k=\pi^2$ in accordance with the constraint (129). The constants ${\cal A}_0^{XY}$, ${\cal
A}_2^{XY}$, and the ${\cal C}^{A,\,XY}$ are determined by the boundary conditions discussed in Sect. 6.

The ${\mathsfsl l}$ = 2 solution

The solution for the l = 2 component is independent of the chosen stellar sequence and is determined by the boundary condition (94) alone. This yields

 \begin{displaymath}\Phi^{XY}_2(r) = - {E^{XY} \over \pi^2}\, \left( r^2 - {5 \over
\sqrt{2}} {J_{5/2}(r\pi)\over \sqrt{r}} \right)\cdot
\end{displaymath} (132)

Inserting this into the l=2 component of (65) gives the corresponding total mass density coefficient

\begin{displaymath}\rho^{XY}_2(r) = - E^{XY} {5 \over \sqrt{2}}\,{J_{5/2}(r\pi) \over
\sqrt{r}} ,
\end{displaymath} (133)

while the respective particle density coefficients are found from (64)

\begin{displaymath}n^{XY}_{A,\,2} = {k_AE^{XY} \over \pi^2} \left(r^2 - {5 \over...
...qrt{2}} {J_{5/2} (r\pi) \over \sqrt{r}}\right) - E^{XY}_Ar^2 .
\end{displaymath} (134)

${\mathsfsl l}$ = 0: FCD sequence

The Fixed Central Density solution for the l=0 component is determined by (94) together with (96). In this case one finds that

\begin{displaymath}\Phi^{XY}_0(r) = {E^{XY} \over \pi^2} \left({3 \sqrt{2} \over...
...\pi) \over \sqrt{r}} + r^2 + {6 \over \pi^2} -
\end{displaymath} (135)

From the l = 0 component of (65) we have

\begin{displaymath}\rho^{XY}_0(r) = - {E^{XY}\over\pi^2} \left( 3\sqrt{2}\,
{J_{1/2}(r\pi) \over \sqrt{r}} - 6 \right) ,
\end{displaymath} (136)

and from (64) it follows that
$\displaystyle n^{XY}_{A,\,0}$ $\textstyle \,=\,$ $\displaystyle - {k_AE^{XY} \over \pi^2} \left({3 \sqrt{2}
\over\pi^2} {J_{1/2}(r\pi) \over \sqrt{r}} + r^2 - {6 \over \pi^2}
    + EXYAr2 . (137)

${\mathsfsl l}$ = 0: FM sequence

The Fixed Mass solution for the l = 0 component is determined by (94) together with (98). One finds

\begin{displaymath}\Phi^{XY}_0(r) = {E^{XY} \over \pi^2} \left(\sqrt{2}\,
{J_{1/2}(r\pi) \over \sqrt{r}} + r^2 - 1\right) ,
\end{displaymath} (138)

and from the l = 0 component of (65) it follows that

\begin{displaymath}\rho^{XY}_0(r) = - {E^{XY} \over \pi^2} \left(\pi^2\sqrt{2}
{J_{1/2}(r\pi) \over \sqrt{r}} - 6\right) ,
\end{displaymath} (139)

while (64) leads to
$\displaystyle n^{XY}_{A,\,0}$ $\textstyle \,=\,$ $\displaystyle - {k_AE^{XY} \over \pi^2}\left(\sqrt{2}
{J_{1/2}(r\pi) \over \sqrt{r}} + r^2 - {6 \over \pi^2} - {3 \over
    $\displaystyle + E^{XY}_A\left( r^2 - {3 \over 5}\right)\cdot$ (140)

8.4 Physical parameters of the two-fluid star

The analytic solution (132)-(140) has been expressed entirely in terms of the matrices EAXY and EXY, and we will now discuss the explicit form of these coefficients in terms of the physical parameters describing the configuration of the two-fluid star.

The symmetric matrix $\S_{XY}$ has three degrees of freedom. One of these is already determined by the static radius constraint (129) and another is associated with the $\sigma$ parameter discussed earlier. A further constraint comes from the proton fraction $x_{\rm p}$, which we define as

\begin{displaymath}x_{\rm p} \equiv {n_{\rm p}^{(0)}\over n_{\rm n}^{(0)}+ n_{\rm p}^{(0)}} \in [0,1] ,
\end{displaymath} (141)

and which is constant throughout the static star due to (126). In order to simplify expressions, we take the proton and neutron masses to be approximately equal, i.e.

\begin{displaymath}m^{\rm n} \approx m^{\rm p} \approx m .
\end{displaymath} (142)

Using (129) and (126), the kX are expressible as

\begin{displaymath}k_{\rm n} = {\pi^2 \over m} (1 - x_{\rm p}) , {\quad{\rm and}\quad}
k_{\rm p} = {\pi^2 \over m} x_{\rm p} .
\end{displaymath} (143)

The "density structure'' SXY can be expressed explicitly in terms of these parameters as

\begin{displaymath}\S_{XY}={\pi^2 \over m^2 (1 + \sigma)}\left(
x_{\rm p} \sigma & ~~~x_{\rm p} \\
\end{array} \right).
\end{displaymath} (144)

From this expression we see that, even though $\sigma\rightarrow
-1$ seemed to be allowed from the "thermodynamical'' point of view, it leads to a singular behaviour in the slow-rotation expansion. Of course, given the fact that the Prakash et al. (1988) equation of state indicated that $\sigma>-0.4$ for all reasonable neutron star densities, cf. Sect. 8.1, we do not expect this singularity to have any physical significance.

As we have seen earlier, cf. (19), the entrainment $\alpha$ can be completely described by the effective proton mass  $m^{\rm p*}$, which for our EOS (119) is a constant, and so we can choose without loss of generality

\begin{displaymath}\beta^{\rm n} = 0 , {\quad{\rm and}\quad}2 \beta^{\rm p} = m \varepsilon ,
\end{displaymath} (145)

where we used dimensionless entrainment coefficient $\varepsilon$ defined earlier in (20).

With the "structure functions'' $\S^{XY}$ and $\beta^X$ expressed completely in terms of the proton fraction $x\rm _p$, the "symmetry energy'' dependence $\sigma$ and the entrainment coefficient  $\varepsilon$, the explicit expressions for EAXY are

$\displaystyle E_{\rm n}^{XY}$ $\textstyle \,=\,$ $\displaystyle {\pi^2 \over 3m (1 + \sigma)}$  
  $\textstyle \times$ $\displaystyle \left(\begin{array}{cc} \left\{1 - x_{\rm p} +
\sigma(1 - 2 x_{\r...
...\sigma \varepsilon & x_{\rm p} \sigma (1 -
\varepsilon) \\  \end{array}\right),$ (146)


 \begin{displaymath}E_{\rm p}^{XY} = {\pi^2 x_{\rm p} \over 3m (1 + \sigma)}
\varepsilon & (1 - \varepsilon) \\
\end{array}\right) ,
\end{displaymath} (147)

while the coefficient EXY of (66) is found to be

 \begin{displaymath}E^{XY} = {\pi^2 \over 3}
1 - x_{\r...
...on & x_{\rm p} (1 - \varepsilon) \\
\end{array}\right) \cdot
\end{displaymath} (148)

It is interesting to note that in the case of co-rotation the corresponding terms in the analytic solution become

\begin{displaymath}\Omega_X E_A^{XY} \Omega_Y= {k_A\over 3 } \Omega^2 ,
\end{displaymath} (149)


\begin{displaymath}\Omega_X E^{XY} \Omega_Y= {\pi^2 \over 3} \Omega^2 ,
\end{displaymath} (150)

and are therefore seen to be independent of not only the entrainment $\varepsilon$, but also of $\sigma$.

8.5 Explicit results for ellipticity, moments of inertia, and Kepler rotation

Using the analytic solution of Sect. 8.3 with the physical parameters of Sect. 8.4, we obtain the following explicit expression for the ellipticities (107)

\begin{displaymath}e_A= {3\over2}\left( 1 - {15 \over \pi^2}\right) \Omega_X E^{XY}
\Omega_Y- {3\pi^2 \over 2 k_A} \Omega_X E_A^{XY} \Omega_Y,
\end{displaymath} (151)

where we have set RA(0)= 1.

For the moments of inertia, integrating (110) for the analytic solution leads to the explicit result

\begin{displaymath}{\delta I_A\over I_A^{(0)}} = a \, \Omega_X E^{XY} \Omega_Y+
{3\,b \over k_A} \, \Omega_X E_A^{XY} \Omega_Y,
\end{displaymath} (152)

with coefficients

\begin{displaymath}a = {9\over \pi^2 - 6}
\left( 3 - {\pi^2 \over 5} - {\pi^4 \over 175}\right),
\end{displaymath} (153)

\begin{displaymath}b = {3 \pi^6 \over 175 (\pi^2 - 6) }\cdot
\end{displaymath} (154)

At the formal level, this result is equivalent to that found in earlier work without entrainment (Prix 1999). However, the two results differ at a deeper level since the actual matrix elements of EXY and EAXY contain information about entrainment and the "symmetry energy'' contribution $\sigma$.

For the analytical solution we can express the correction to the Kepler rotation rate (114) as

$\displaystyle %
\delta\Omega_{{\rm K},{\it A}}^2$$\textstyle \,=\,$ $\displaystyle {6 \over \pi^2}\left({1 \over 5} - {3 \over \pi^2}
\right) \Omega_X E^{XY} \Omega_Y- {27 \over 10
k_A} \Omega_X E_A^{XY} \Omega_Y,$  

where A denotes the "outer'' fluid at the equator.

If we take the simple case of the two fluids co-rotating at the maximal rotation rate, i.e. $\Omega_{\rm n} = \Omega_{\rm p} = \Omega_{\rm K}$, this simply leads to

\begin{displaymath}\Omega_{\rm K}= {\Omega_{(0)}\over \sqrt{{6 \over \pi^2} + {3...
...0.69 \, \Omega_{(0)}\approx 0.80\, \sqrt{\pi G
\bar{\rho}} ,
\end{displaymath} (156)

independent of the parameters of this EOS. This is about 20%-25% too high compared to the numerical result $\Omega_{\rm K}\approx0.64\sqrt{\pi G \bar{\rho}}$. This discrepancy originates in the underestimate of the equatorial radius when approaching $\Omega_{\rm K}$, which is illustrated in Fig. 1. The "failure'' of the slow-rotation expansion to accurately determine the Kepler rotation limit has already been pointed out by, for example, Salgado et al. (1994).

9 Exploring the roles of relative rotation, entrainment and "symmetry energy''

In this section we will use the analytic solution to the superfluid slow-rotation problem to explore the effects of relative rotation, entrainment, and "symmetry energy'' on the distribution of matter, the Kepler limit, ellipticity, and moments of inertia for a fixed mass star. The results we present are for a particular stellar model with mass $1.4~M_\odot$ and radius 10 km. The proton fraction is taken to be $x_{\rm p} = 0.1$ and we choose the neutron rotation rate to be (except in Figs. 4 and 7) $\Omega _{\rm n} = 0.1$. This value corresponds to (in our units) the rotation rate of the fastest known pulsar (1.6 ms). The three parameters that will be varied are $\sigma$, $\varepsilon$, and the relative rotation rate $\Omega _{\rm p}/\Omega _{\rm n}$. For the first two parameters we will consider only the values $\sigma = -0.5,0,0.5$ and $\varepsilon = 0,0.4,0.7$, while the relative rotation rate will be assumed to lie in the range $-2 \leq
\Omega_{\rm p}/\Omega_{\rm n} \leq 2$ (which allows the protons to be counter-rotating with respect to the neutrons). The values chosen for $\sigma$ are in accordance with the discussion in Sect. 8.1. We also recall that the best current estimates for entrainment imply a range of $0.4 \leq \varepsilon \leq 0.7$. This means that our chosen models correspond to the expected upper and lower limits, as well as the no-entrainment limit (which provides a useful reference).

\par\includegraphics[width=16cm,height=6cm,clip]{Prix1666f2.eps} \end{figure} Figure 2: Plots of the radial profiles of the proton density corrections, normalized by the static number density, in the equatorial plane ( $\theta = \pi /2$) and along the polar axis ( $\theta = 0$) for $\sigma = 0$ (left graph), $\sigma =0.5$ (right graph), $\varepsilon = 0,0.4,0.7$, and $\Omega _{\rm n} = 0.1$ and $\Omega_{\rm p} =
0.25~ \Omega_{\rm n}$.
Open with DEXTER

9.1 Effects on the "local'' structure

First we focus on the "local'' structure by examining how the particles get redistributed throughout the star because of rotation. Figure 2 shows the rotationally induced corrections to the proton particle number densities in the equatorial plane and along the polar axis, for $\sigma = 0$ and $\sigma =0.5$. In Fig. 3 we show the corresponding corrections to the neutron number densities. For the neutrons we show plots for $\sigma = 0$ only because other values lead to similar results. An inset is provided to show that the effect of entrainment is small but not completely negligible. Recall that we are considering the fixed mass solution, so the density at the centre of the star will decrease as the star is spun up. This also explains why the density correction is negative along the polar axis. Furthermore, the density corrections in the equatorial plane are positive near the stellar surface since rotation causes the matter to bulge.

\par\includegraphics[width=13.4cm,clip]{Prix1666f3.eps} \end{figure} Figure 3: Plots of the radial profiles of the neutron density corrections, normalized by the static number density, in the equatorial plane ( $\theta = \pi /2$) and along the polar axis ( $\theta = 0$) for $\sigma = 0$, $\varepsilon = 0,0.4,0.7$, and $\Omega _{\rm n} = 0.1$ and $\Omega_{\rm p} =
0.25~ \Omega_{\rm n}$. In this, and the following figures, the dotted lines correspond to $\varepsilon = 0.7$, dashed have $\varepsilon = 0.4$, whereas the solid lines are for $\varepsilon = 0$.
Open with DEXTER

From Fig. 2 we notice that the effect of changing $\sigma$ and $\varepsilon$ is more pronounced for the protons. That this should be the case is easy to understand since the proton fraction is small and only $10\%$ of the mass of the star is in the protons. The relative effect of a changing $\varepsilon$ is more apparent in the proton plots, then, because the absolute magnitude of the proton number density corrections are always $10\%$ of those of the neutrons. From (140), (146) and (147) we see that the main difference between the two cases is that while changes in $\varepsilon$ and $\sigma$ affect $E_{\rm p}^{XY}$ at leading order ( ${\sim} x_{\rm p}$), they only lead to higher order corrections to  $E_{\rm n}^{XY}$(since there are also terms of O(1) in (146)). We can also see from the right-hand panel of Fig. 2 that changing $\sigma$ from 0 to 0.5 results in substantial modifications. This shows clearly that the "symmetry energy'' can play an important role in determining the rotational configuration of a two-fluid star.

\par\includegraphics[width=16cm,clip]{Prix1666f4.eps} \end{figure} Figure 4: Neutron (thick line), proton (thin line) and static configuration (dashed line) isodensity curves in the meridional plane, for $\sigma = 0,- 0.5$, $\varepsilon = 0,0.7$, and $\Omega _{\rm n} = 0.15$ and $\Omega_{\rm p} =
0.25~ \Omega_{\rm n}$.
Open with DEXTER

\par\includegraphics[width=16.8cm,clip]{Prix1666f5.eps} \end{figure} Figure 5: The neutron and proton ellipticities as functions of the relative rotation $\Omega _{\rm p}/\Omega _{\rm n}$, for $\sigma = -0.5,0,0.5$, $\varepsilon = 0,0.4,0.7$, and $\Omega _{\rm n} = 0.1$.
Open with DEXTER

\par\includegraphics[width=16.9cm,clip]{Prix1666f6.eps} \end{figure} Figure 6: Neutron (${\rm n}$) and proton (${\rm p}$) moments of inertia for $\sigma = -0.5,0,0.5$, $\varepsilon = 0,0.4,0.7$, and $\Omega _{\rm n} = 0.1$.
Open with DEXTER

Figure 4 provides a different perspective on the rotationally-induced redistribution of the particles. Here we show isodensity surfaces in the $(r,\theta)$-plane, for $\varepsilon = 0$, $\sigma = 0$ (in the left panel) and $\varepsilon = 0.7$, $\sigma = - 0.5$ (in the right panel). Note that we have considered slightly larger values of the rotation rates (just below the Kepler limit) in order to exaggerate certain effects. For a given density, we compare the isodensity curves for a non-rotating star to the neutron and proton isodensity curves for the rotating model. In both panels we see that the neutron and proton curves actually intersect each other near the surface of the star. That is, along the equator the neutrons actually extend further than the protons, while the protons extend further along the rotation axis. This means that, in this case the rotational configuration of the protons is, in fact, prolate. We note that this effect has been exaggerated in the panels, because of the higher rotation rates, but it happens also for lower rotation rates (as used in the earlier figures). Near the center of each panel, we see that the neutron and proton surfaces no longer intersect. This is no doubt due to the fact that the centrifugal forces are smaller closer to the center of the star.

9.2 Effects on the "global'' structure

Now we focus our attention on the roles of $\sigma$, $\varepsilon$, and the relative rotation rate $\Omega _{\rm p}/\Omega _{\rm n}$ in determining macroscopic properties of the star; in particular, the proton and neutron ellipticities, their respective moments of inertia, and the Kepler limit. As before, we will keep the mass of the star fixed.

Figure 5 illustrates the neutron and proton ellipticities as functions of the relative rotation rate, and Fig. 6 gives their moments of inertia (normalized by their static values). There is an obvious quadratic behavior in each plot due simply to the slow-rotation expansion. As well, the intersection of all curves at $\Omega_{\rm p}/\Omega_{\rm n} = 1$ occurs because the protons then co-rotate with the neutrons and the system is behaving as a single fluid. Notice that entrainment has the largest influence when the neutrons and protons counter-rotate. This is easily understood as a consequence of the basic fact that the entrainment parameter represents the way that the equation of state (as represented by the energy functional) depends on $\vert\vec{v}_{\rm n} - \vec{v}_{\rm p}\vert^2$. The protons are in general much more affected by changes in the various parameters than the neutrons, again due to the fact that the neutrons carry the bulk of the mass of the star. Perhaps most interesting is the effect of both the entrainment parameter $\varepsilon$ and "symmetry energy'' parameter $\sigma$ in determining the minima of the curves. For both the protons and the neutrons, we see that an increase in $\varepsilon$ for a fixed value of $\sigma$ leads to a deeper value for the minimum. Decreasing the value of $\sigma$ causes the minima to become even deeper. In particular, we see from the left-most panel in Fig. 5 that the minima have become deep enough that the protons can be prolate (i.e. have negative ellipticity) even though they rotate in the same direction as the neutrons. Finally, we note that as $\sigma$ is decreased, the neutron ellipticities go from having minima in the right-most and center panels, to having maxima in the left-most panel. That is, if the absolute value of the relative rotation could be made large enough the neutron fluid could also become prolate.

\par\includegraphics[width=17.3cm,clip]{Prix1666f7.eps} \end{figure} Figure 7: Plots of the neutron (n) and proton (p) Kepler limits as functions of the relative rotation $\Omega _{\rm n}/\Omega _{\rm p}$, for $\sigma = -0.5,0,0.5$ and $\varepsilon = 0,0.4,0.7$.
Open with DEXTER

It is interesting to compare the curves for the ellipticities in Fig. 5 and the corresponding moments of inertia in Fig. 6, for example for $\sigma=0,\,\varepsilon=0.4$ (middle panel, dashed line) for the protons. For $\Omega_{\rm p}=0$ we observe that the proton fluid is nearly spherical ( $e_{\rm p}\approx0$), and still its moment of inertia is higher than of the static (spherical) configuration $\delta I_{\rm p}>0$, and this is simply because the mass distribution has been shifted further away from the rotation axis even though the surface itself is (nearly) unchanged in this case. This can also be clearly seen in the density correction in Fig. 2, for example for $\sigma=0,\,\varepsilon=0.7$ and $\Omega_{\rm p}/\Omega_{\rm n}=0.25$ (left panel, dotted line), which nearly vanishes at the surface r=1, and is to be compared to the corresponding global quantities in Figs. 5 and 6.

Finally, results for the Kepler, or mass-shedding, limit are shown in Fig. 7. To understand these results one must appreciate that there is a subtle difference from the single-fluid case: In our case the two fluids can rotate independently at different rates. Thus, one of the fluids typically extends beyond the other, in particular at the equator. Since the Kepler limit is defined by the outermost fluid at the equator, we can use Eq. (112) in the following way: when the neutrons are outermost, set $A = {\rm n}$ and $\Omega_{\rm n} = \Omega_{K,{\rm n}}$ and solve the resulting quadratic for $\Omega_{K,{\rm n}}$ as a function of the ratio $\Omega _{\rm n}/\Omega _{\rm p}$, and vice versa in the case when the protons extend beyond the neutrons. In Fig. 7 we show the resultant solutions over the entire range of the relative rotation rate. The Kepler rate is easy to determine, however, because it is given by the neutron curves when $\Omega_{\rm n}/\Omega_{\rm p} \ge 1$, and the proton curves when $\Omega_{\rm n}/\Omega_{\rm p} \le 1$. Of course, the various curves always intersect when $\Omega_{\rm n}/\Omega_{\rm p} = 1$, the case that corresponds to corotation of the two fluids. For the case of $\sigma = \epsilon = 0$, we find results in good qualitative agreement with the relativistic study of Andersson & Comer (2001a). In particular, we see that the Kepler limit changes little when $\Omega_{\rm n} \ge
\Omega_{\rm p}$. As explained by Andersson & Comer (2001a), this is due to the fact that the neutrons make up $90\%$ of the mass of the star, and the star is behaving like a single-fluid star with a small proton component. When $\Omega_{\rm p} \ge \Omega_{\rm n}$ the Kepler limit increases as $\Omega_{\rm n}$ is decreased. Again, this is natural because the neutrons still dominate the mass of the star, and the Kepler rate is simply approaching the non-rotating star limit.

10 Conclusions

We have developed a formalism for modeling slowly-rotating Newtonian superfluid neutron stars incorporating entrainment. We have used a two-fluid description, where one fluid is the superfluid neutrons and the other is a charge-neutral conglomerate of the remaining constituents. A detailed discussion of the relation between entrainment and nuclear physics calculations (i.e. equations of state) was given. Using an equation of state that is quadratic in both the mass-densities and relative velocities of the fluids, we found that an analytic solution to the slow-rotation equations could be obtained. This solution is the natural extension to the two-fluid case of the ${\cal E} \propto \rho^2$ polytrope in the single fluid case (which has proven to be very useful for understanding the properties of ordinary fluid neutron stars). We used the analytic solution to explore effects due to relative rotation, entrainment, and "symmetry energy'' on both the "local'' and "global'' properties of a fixed-mass star. An unexpected result was that the "symmetry energy'' parameter had as much impact on the rotational equilibria as the entrainment parameter.

Our ultimate goal is to study the modes of oscillation of both Newtonian and general relativistic slowly rotating superfluid neutron stars. We believe that the formalism and analytic solution presented here will be valuable in reaching this goal. In particular, the inclusion of entrainment is absolutely necessary in determining how the dominant dissipative mechanism (the so-called mutual friction) of rotating superfluid neutron stars affects the gravitational radiation emitted from unstable modes.

We would like to thank the members of the Observatory of Paris-Meudon Numerical Relativity group for providing us with their LORENE code, and J. Novak for fruitful conversations. NA is Leverhulme Prize fellow and acknowledges support from PPARC via grant number PPA/G/S/1998/00606.

RP and NA acknowledge support from the EU Programme "Improving the Human Research Potential and the Socio-Economic Knowledge Base'' (Research Training Network Contract HPRN-CT-2000-00137).

GC gratefully acknowledges partial support from a Saint Louis University SLU2000 Faculty Research Leave Award as well as EPSRC in the UK via grant number GR/R52169/01 (to NA), and the warm hospitality of the Center for Gravitation and Cosmology of the University of Wisconsin-Milwaukee and the University of Southampton where part of this research was carried out.

Appendix A: An alternative derivation of the first integrals of motion

A more elegant way of finding the first integrals of motions (31) consists of using the identity

\begin{displaymath}v^{j}\nabla_{j} \,p_{i} = {\cal L}_v p_{i} - p_{j} \nabla_{i} v^{j},
\end{displaymath} (A.1)

in terms of the Lie derivative ${\cal L}_v$, which allows one to rewrite Euler's Eqs. (12) in the form
$\displaystyle \partial_t p^{\rm n}_{i}$ + $\displaystyle {\cal L}_{v_{\rm n}} p^{\rm n}_{i} - \nabla_{i} \left(p^{\rm n}_0+
v_{\rm n}^{j} p^{\rm n}_{j} \right) = 0,$ (A.2)
$\displaystyle \partial_t p^{\rm p}_{i}$ + $\displaystyle {\cal L}_{v_{\rm p}} p^{\rm p}_{i} - \nabla_{i} \left(p^{\rm p}_0 +
v_{\rm p}^{j} p^{\rm p}_{j} \right) = 0.$ (A.3)

Stationarity ( $\partial_t p_{i}^X=0$), uniform rotation ( $v_X^{i} = \Omega_X\varphi^{i}$) and axial symmetry ( ${\cal L}_\varphi p_{i}^X=0$) then directly results in the integrals
$\displaystyle \nabla_{i} \left(p^{\rm n}_0 + v_{\rm n}^{j} p^{\rm n}_{j} \right)$ $\textstyle \,=\,$ 0, (A.4)
$\displaystyle \nabla_{i} \left(p^{\rm p}_0 + v_{\rm p}^{j} p^{\rm p}_{j} \right)$ $\textstyle \,=\,$ 0 (A.5)

and inserting (9) this explicitly becomes
$\displaystyle \mu^{\rm n} + m^{\rm n} \Phi - {1\over2}m^{\rm n} v_{\rm n}^{\,2}$ $\textstyle \,=\,$ $\displaystyle {\cal C}^{\rm n} ,$ (A.6)
$\displaystyle \mu^{\rm p} + m^{\rm p} \Phi - {1\over2}m^{\rm p} v_{\rm p}^{\,2}$ $\textstyle \,=\,$ $\displaystyle {\cal C}^{\rm p}.$ (A.7)


Copyright ESO 2001