A&A 381, 178-196 (2002)
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 . 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
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 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.
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 "'' and "'', representing the
neutrons and protons, respectively. The fundamental variables of our
description are the respective particle number densities
and the corresponding particle currents
We will only consider situations where the particle
numbers are conserved individually, so that we have
The internal energy density or "equation of state''
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
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
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,
say, of particles with momentum
to the non-vacuum background. The velocity
of the particles as
a function of their momentum is then given by
and the acceleration
dot represents the total time derivative) is linked to the force
via a relation of the form
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 (where represents nuclear density, i.e. g/cm3), that ; Sjöberg (1976) finds, for , that ; and more recently, Baldo et al. (1992) find, for the same range as Chao et al., that . Thus, we see that the proton effective mass apparently can range over the values . 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.
Stationary two-fluid configurations are
and therefore the equations of
motion (12) take the form
It is interesting to note that the stationary two-fluid model
with different rotation rates
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
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
for densities, and
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
For realistic equations of state the mass shedding limit is well
approximated by the empirical expression
is the mean density.
Therefore the dimensionless rotation rates
generally smaller than 1, and for typical rotation rates
most observed pulsars we will have
|Figure 1: Comparison of one-fluid slow-rotation configurations with numerical results obtained using the LORENE code. and are the star's equatorial and polar radius repectively. The stellar model is a N=1 polytrope with mass 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|
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 and being small compared to the natural rotation rate . In the dimensionless form of the equations, the small parameters of the expansion are simply the rotation rates, as and . We can therefore write the solution as a Taylor expansion in and . 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
We know that the static solution
is purely spherical and
We now expand all physical quantities of (31) and
(34) up to second-order in
We can find an algebraic relation between the chemical potential
and the density coefficients nAXY,
which will allow us to considerably simplify the problem. This relation
is found by expanding the chemical potential
defined in (5) in its arguments up to
Using (44) we find
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 .
The total mass density coefficient
can now be written
We note that a necessary condition for this method to work, i.e. for
(67) to be well-defined, is that the structure function
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.
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
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
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
We now separate the variables r and
in the orthogonal basis of Legendre polynomials
These are the eigenfunctions of the angular part of the
Using (73) together with the identity
Once the solutions 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).
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
First, we require that the
solution must be regular at the centre of the star. At r=0 the
is singular, which leads to the first
We write the slow-rotation expansion of the exterior potential
the usual way as
Up to second-order in
the surface of the rotating star is
As stated earlier, the l=0 equation requires two additional conditions in order to fix the constants of integration and . These conditions are crucial as they determine which type of rotating star sequence (as a function of the rotation rates and ) 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
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
to second-order in the rotation rates. Starting from the obvious
The respective moments of inertia
IA to second-order can be written as
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 and (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.
Within the present slow-rotation approximation any equation of state can
(quite generally) be expressed as
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 . 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
which leads to the term
It will turn out that
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
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
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
The Fixed Central Density solution for the l=0 component is determined by
(94) together with (96). In this case one finds that
|+ EXYAr2 .||(137)|
The Fixed Mass solution for the l = 0 component is determined by
(94) together with (98). One finds
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
has three degrees of freedom. One of
these is already determined by the static radius constraint (129) and another is associated with
the parameter discussed earlier. A
further constraint comes from the proton fraction ,
which we define as
As we have seen earlier, cf. (19), the entrainment
can be completely described by the effective proton mass
for our EOS (119) is a constant, and so we can choose
without loss of generality
With the "structure functions''
completely in terms of the proton fraction ,
the "symmetry energy''
and the entrainment coefficient
the explicit expressions for EAXY are
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)
For the moments of inertia,
integrating (110) for the analytic solution leads to the
For the analytical solution we can express the correction to
the Kepler rotation rate (114) as
If we take the simple case of the two fluids co-rotating at the maximal
rotation rate, i.e.
this simply leads
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
and radius 10 km.
The proton fraction is taken to be
choose the neutron rotation
rate to be (except in Figs. 4 and 7)
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 ,
and the relative rotation rate
For the first two parameters
we will consider only the
relative rotation rate will be assumed to lie in the range
(which allows the protons to be
counter-rotating with respect to the neutrons). The values chosen for
are in accordance with the discussion in
Sect. 8.1. We also recall that the best current
estimates for entrainment
imply a range of
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).
|Figure 2: Plots of the radial profiles of the proton density corrections, normalized by the static number density, in the equatorial plane ( ) and along the polar axis ( ) for (left graph), (right graph), , and and .|
|Open with DEXTER|
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
In Fig. 3 we show the corresponding corrections to the
neutron number densities.
For the neutrons we show plots for
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.
|Figure 3: Plots of the radial profiles of the neutron density corrections, normalized by the static number density, in the equatorial plane ( ) and along the polar axis ( ) for , , and and . In this, and the following figures, the dotted lines correspond to , dashed have , whereas the solid lines are for .|
|Open with DEXTER|
From Fig. 2 we notice that the effect of changing and is more pronounced for the protons. That this should be the case is easy to understand since the proton fraction is small and only of the mass of the star is in the protons. The relative effect of a changing is more apparent in the proton plots, then, because the absolute magnitude of the proton number density corrections are always 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 and affect at leading order ( ), they only lead to higher order corrections to (since there are also terms of O(1) in (146)). We can also see from the right-hand panel of Fig. 2 that changing 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.
|Figure 4: Neutron (thick line), proton (thin line) and static configuration (dashed line) isodensity curves in the meridional plane, for , , and and .|
|Open with DEXTER|
|Figure 5: The neutron and proton ellipticities as functions of the relative rotation , for , , and .|
|Open with DEXTER|
|Figure 6: Neutron () and proton () moments of inertia for , , and .|
|Open with DEXTER|
Figure 4 provides a different perspective on the rotationally-induced redistribution of the particles. Here we show isodensity surfaces in the -plane, for , (in the left panel) and , (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.
Now we focus our attention on the roles of , , and the relative rotation rate 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 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 . 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 and "symmetry energy'' parameter in determining the minima of the curves. For both the protons and the neutrons, we see that an increase in for a fixed value of leads to a deeper value for the minimum. Decreasing the value of 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 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.
|Figure 7: Plots of the neutron (n) and proton (p) Kepler limits as functions of the relative rotation , for and .|
|Open with DEXTER|
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 and and solve the resulting quadratic for as a function of the ratio , 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 , and the proton curves when . Of course, the various curves always intersect when , the case that corresponds to corotation of the two fluids. For the case of , 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 . As explained by Andersson & Comer (2001a), this is due to the fact that the neutrons make up of the mass of the star, and the star is behaving like a single-fluid star with a small proton component. When the Kepler limit increases as 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.
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 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.
A more elegant way of finding the first integrals of motions
(31) consists of using