On the oscillations of dissipative superfluid neutron stars

On the oscillations of dissipative superfluid neutron stars

N. Andersson, K. Glampedakis and B. Haskell School of Mathematics, University of Southampton, Southampton SO17 1BJ, United Kingdom SISSA, via Beirut 2-4, 34014 Trieste, Italy Theoretical Astrophysics, Auf der Morgenstelle 10, University of Tuebingen, Tuebingen D-72076, Germany

We investigate the oscillations of slowly rotating superfluid stars, taking into account the vortex mediated mutual friction force that is expected to be the main damping mechanism in mature neutron star cores. Working to linear order in the rotation of the star, we consider both the fundamental f-modes and the inertial r-modes. In the case of the (polar) f-modes, we work out an analytic approximation of the mode which allows us to write down a closed expression for the mutual friction damping timescale. The analytic result is in good agreement with previous numerical results obtained using an energy integral argument. We extend previous work by considering the full range of permissible values for the vortex drag, e.g. the friction between each individual vortex and the electron fluid. This leads to the first ever results for the f-mode in the strong drag regime. Our estimates provide useful insight into the dependence on, and relevance of, various equation of state parameters. In the case of the (axial) r-modes, we confirm the existence of two classes of modes. However, we demonstrate that only one of these sets remains purely axial in more realistic neutron star models. Our analysis lays the foundation for companion studies of the mutual friction damping of the r-modes at second order in the slow-rotation approximation, the first time evolutions for superfluid neutron star perturbations and also the first detailed attempt at studying the dynamics of superfluid neutron stars with both a relative rotation between the components and mutual friction.


I Introduction

Neutron stars have a complex interior structure. With core densities reaching several times the nuclear saturation density, these objects require an understanding of physics that cannot be gained from laboratory experiments. This makes the modelling of neutron star dynamics an interesting challenge. On the one hand, one has to consider exotic physics that is, at best, poorly constrained. On the other hand, one may ask to what extent observations can distinguish between different possible models. An excellent example of this interplay concerns the possibility that the quarks may deconfine in the high density region. If this is the case, it will have a considerable effect on transport properties associated with viscosity and heat conductivity. In fact, such a quark core is expected to be a colour superconductor alford . The dynamics of this exotic state of matter, and the relevance of its different possible phases, is not yet certain. In order to improve our understanding of this problem, we need to build more precise stellar models and study, for example, their oscillation properties in detail. In this context, considerable attention has been focused on the inertial r-modes of a rotating star. The r-modes are interesting because they can be driven unstable by the emission of gravitational radiation, see nk ; nareview for literature reviews. The r-mode instability window is, however, sensitive to the physics of the neutron star interior. Since the bulk and shear viscosities are quite different in a quark core, compared to “normal” npe matter, one may hope to use observations to constrain the theory, see alford for a discussion of the relevant literature. In absence of a direct detection of an r-mode gravitational-wave signal, this analysis would have to be based on the nature of the instability window. The idea would be that, if an observed neutron star spins at a rate that would place it inside a predicted instability region, one may be able to rule out this particular theoretical model. Of course, this argument comes with a number of caveats. It could, for example, be that additional physics places a stronger constraint on the r-modes than the considered mechanisms. Inevitably, this becomes a “work in progress” where improved theoretical models are tested against better observational data.

In order to consider “realistic” neutron stars, it is important to appreciate the relevance of superfluidity. A neutron star is expected to contain a number of superfluid/superconducting components haensel , and it is crucial to understand to what extent this affects the stars oscillation properties. It is well established that the behaviour of a superfluid system can differ significantly from standard hydrodynamics. The most familiar low-temperature system is, perhaps, He, which exhibits superfluidity below a critical temperature near 2 K. Experimentally, it has been demonstrated that this system is very well described by the Navier-Stokes equations above the critical temperature. Below the critical temperature the behaviour is very different, and a “two-fluid” model is generally required (see helium for a very recent discussion). Superfluid neutron stars are, to some extent, similar. The second sound in Helium is analogous to a set of, more or less distinct, “superfluid” oscillation modes epstein ; mendell1 ; ac01 in a neutron star. These additional modes arise because the different components of a superfluid system are allowed to move “through” each other. The dissipation channels in a superfluid star are also quite different. Basically, the superfluid flows without friction. In the outer core of a neutron star, which is dominated by npe matter, one expects the neutrons to be superfluid while the protons form a superconductor. As a result, the shear viscosity is dominated by e-e scattering cutler ; npa . The bulk viscosity, which is due to the fluid motion driving the system away from chemical equilibrium and the resultant energy loss due to nuclear reactions, is also expected to be (exponentially) suppressed in a superfluid haensel . These effects have direct implications for the damping of neutron star oscillations, and play a key role in determining the r-mode instability window for a mature neutron star. This is, however, not the end of the story. A superfluid exhibits an additional dissipation mechanism, usually referred to as “mutual friction”. The mutual friction is due the presence of vortices in a rotating superfluid. In a neutron star core, the electrons can scatter dissipatively off of the (local) magnetic field of each vortex (see als ; mendell ; trev1 for discussions and references). This effect may dominate the damping of realistic neutron star oscillation modes.

The basic requirements of a rudimentary model for superfluid neutron star oscillations should now be clear. One must account for the additional dynamical degree(s) of freedom, and also account for the mutual friction damping. This is obviously only a starting point, but the problem is sufficiently complicated that one may want to proceed with care. There has already been a number of studies of dissipative superfluid oscillations. The area was pioneered by Lindblom and Mendell who, in particular, demonstrated that the gravitational-wave instability of the fundamental f-modes would be suppressed in a superfluid star lm95 . Following the discovery of the r-mode instability, they also provided the first accurate estimates of the relevance of the mutual friction for these modes lm00 . Similar results were subsequently obtained by Lee and Yoshida yl1 . These studies provide important assessments of the relevance of the mutual friction damping. There are, however, a number of reasons why we need to return to this problem. Most importantly, we want to consider more realistic neutron star models, including finite temperature effects, magnetic fields and the possible presence of exotic (hyperon and/or quark) cores. The additional physics brings additional complications, like additional fluid degrees of freedom, boundary layers at phase-transition interfaces and fundamental issues concerning dissipative multifluid systems helium ; monster . We also need to move away from the assumption that the vortex drag, which leads to the mutual friction, is weak. Strong arguments suggest that this is not going to be the case when the protons form a type II superconductor and there are magnetic fluxtubes present in the system sauls ; ruderman ; link03 . The neutron vortices may be “pinned” to the fluxtubes leading to a strong drag regime. The strong drag problem has only been considered recently precess ; kostas , and the first results demonstrate the presence of a new instability in systems where the two components rotate at different rates. This instability, which may be relevant for the understanding of pulsar glitches (kostas, ), provides a direct demonstration that the dynamics in the strong drag regime may be both complicated and interesting. The present investigation lays the foundation for future work in this direction by allowing for strong drag. In particular, we retain the dynamic contribution to the mutual “friction” that has previously been neglected as a matter of course.

Ii The two-fluid equations

Our discussion is based on the standard two-fluid model for neutron star cores prix ; monster . That is, we consider two dynamical degrees of freedom loosely speaking representing the superfluid neutrons (labeled ) and a charge-neutral conglomerate of protons and electrons (labeled ). Assuming that the individual species are conserved, we have the usual conservation laws for the mass densities ,


where the constituent index x may be either p or n. Meanwhile, the equations of momentum balance can be written


where the velocities are , the relative velocity is defined as and represents the chemical potential (we will assume that throughout this paper). represents the gravitational potential, and the parameter encodes the non-dissipative entrainment coupling between the fluids prix ; monster . The force on the right-hand side of (2) can be used to represent various other interactions, including dissipative terms monster .

In the following we will focus on the vortex-mediated mutual friction. Assuming that the two fluids exhibit solid body rotation we have a force of form trev1 (see also als ; mendell )


Here, is the angular frequency of the neutron fluid (a hat represents a unit vector). The mutual friction parameters are intimately related to the induced friction on the vortex. The latter is often described in terms of a dimensionless “drag” parameter such that


In the standard picture, the mutual friction is due to the scattering of electrons off of the array of neutron vortices. This leads to , i.e., , and hence the first term in the mutual friction force can be ignored. There are, however, good arguments for why the problem may be in the opposite regime. In particular if one considers the interaction between the fluxtubes in a type II proton superconductor and the neutron vortices sauls ; ruderman ; link03 . Then one would expect to be in the strong drag regime where , i.e., while remains small. Superfluid oscillations in this regime have not previously (with the exception of kostas ) been considered.

Anyway, from (3) we see that the mutual friction will not be present in a non-rotating star. This is obvious since there would then be no vortices in the first place. Of course, any non-trivial motion of the superfluid neutrons leads to vortex generation. This means that a generic perturbation of a non-rotating star will be associated with a local vorticity which could lead to mutual friction. However, in this context the resulting mutual friction interaction would require a perturbative calculation to be carried out to second order. As far as we are aware, such calculations have not yet been attempted. It may be an interesting problem for the future.

Iii The perturbation equations

iii.1 Decoupling the degrees of freedom

If we want to consider the effects of mutual friction we need to consider rotating stars. To keep the problem tractable (at least initially) we assume that the background configuration is such that the two fluids rotate together. Perturbing the equations of motion and working in a frame rotating with we then have




where represents an Eulerian variation.

From previous work on superfluid neutron star perturbations (and indeed the large body of work on superfluid Helium) we know that the problem has two ”natural” degrees of freedom, see for example epstein ; mendell1 ; ac01 ; prix1 ; prix2 ; trev2 . One of the degrees of freedom represents the total mass flux. Introducing


and combining the two Euler equations we find that


where and the pressure is obtained from 111Caution: In the general case when the background fluids are not co-rotating, there will be additional contributions associated with the entrainment in this relation.


In deriving this relation we have used


where it has been assumed that the two fluids are in chemical equilibrium in the background. That is, we have . We also have the usual continuity equation


At this point we have two equations which are identical to the perturbation equations for a single fluid system. It is particularly notable that (8) does not have a force term. This follows immediately from the fact that we are only considering the mutual friction interaction. In other situations, say including shear viscosity, we would no longer have a homogeneous equation.

Of course, we are considering a two-fluid problem and there is a second degree of freedom to take into account. To describe this, it is natural to consider the difference in velocity. Thus, we introduce


Combining the two Euler equations in the relevant way we have


Here we have defined


which represents the (local) deviation from chemical equilibrium induced by the perturbations. We have also introduced the simplifying notation


where is the proton fraction. Again, equation (13) does not couple the different degrees of freedom.

The coupling is entirely due to the second continuity equation. It is natural to use the proton fraction to complement the total density . Then we find that


This equation shows that the two dynamical degrees of freedom are explicitly coupled unless the proton fraction is constant. This fact has already been pointed out by Prix and Rieutord prix1 .

Before moving on, it is worth asking to what extent it is possible to find solutions that are purely co-moving, i.e. for which . From the above equations it is easy to see that such a solution would have to satisfy


This condition is trivially satisfied if the proton fraction is uniform. In addition, it will be satisfied for fluid motion that has (for a spherical background configuration) no radial component and also do not lead to variations in . Are there oscillation modes with this character? Indeed, to leading order in the slow-rotation approximation the standard r-mode satisfies these criteria. It is purely axial and the associated density perturbations appear at order . However, in general we do not expect to find any oscillations of a “realistic” neutron star model to be purely co-moving. This means that a generic neutron star oscillation mode will be affected by mutual friction.

iii.2 Boundary conditions

To completely specify the perturbation problem, we need boundary conditions. At the centre of the star we simply require that all variables are regular. The surface of the star is somewhat more complex. In reality one does not expect the superfluid region to extend all the way to the surface. A real neutron star will always be covered by a single fluid envelope (eg. the outer parts of the elastic crust). However, for simplicity we do not want to deal with the various interfaces in the present analysis lm00 ; layers ; lin . Instead we will consider stars with a two-fluid surface, which is obviously somewhat artificial.

A reasonable approach is to assume that the perturbed star has a unique surface. That is, let the two perturbed fluids move together (in the radial direction) at the surface. In the two-fluid problem we have two distinct Lagrangian displacements grosart . These follow from


We have assumed that the two fluids corotate in the background, i.e. we have . If we also impose that there is a common surface, then we have at and it follows that we should require


From (13) we see that this implies that, for a non-rotating configuration we must also have


When we determine the rotational corrections to the f-mode we will impose this condition also at first slow-rotation order. This is not entirely consistent, cf. (13), but it is straightforward to relax this condition should it be required.

If the two fluids move together at the surface, it also follows that


where at the surface. This is the usual single fluid condition of a vanishing Lagrangian pressure variation .

iii.3 A bit of chemistry

Consider the various equations that we have written down. At this point the two degrees of freedom and only couple explicitly through (16). In fact, if we assume that the two fluids are incompressible, then there is no coupling at all. Since the mutual friction only enters the problem via (13), it is thus the case that any incompressible dynamics in the sector will be unaffected by mutual friction. This shows that, if we are interested in the effect of mutual friction on (say) the f-mode oscillations of a star it is not meaningful to consider an incompressible model. We know already from the outset that we would only find the mutual friction effects on the counter-moving “superfluid” modes. This may be an interesting problem, but it is not our main motivation here.

For compressible models, the two degrees of freedom also couple indirectly. Basically, we need to use the equation of state to relate to . For models where the two fluids co-rotate in the background we can use 222Note that there will be entrainment contributions here in the general case when the two background fluids are not co-rotating.




Using these relations, or their “inverse”, we see that the two degrees of freedom couple in a more subtle way. If we choose to reduce the problem by eliminating and then the coupling arises through the boundary conditions and the Euler equations. If, on the other hand, we eliminate and then the coupling enters through the continuity equations.

Iv Dissipation integrals

In order to estimate the damping associated with various dissipation mechanisms one can either (i) account for the dissipative terms in the equations of motion and solve for the damped modes directly, or (ii) solve the non-dissipative problem and use an energy integral argument to estimate the damping rate. In the typical situation when the damping is very slow the second strategy should be reliable. Indeed, all studies of damped neutron star oscillations have used this approach (see nareview for a discussion). Given this, it is natural to pause and consider the energy integral approach to the mutual friction problem.

iv.1 The conserved energy

To work out a suitable energy associated with a given perturbation, we first multiply (5) with (where the bar represents complex conjugation). Then we add the result to its complex conjugate. Combining the individual contributions from the neutron and proton fluids and integrating over the star we find that, when , the result is a total time derivative of two terms. The first term is the “kinetic energy”, which follows from


where . Alternatively, expressing this in the co- and countermoving variables, we have


The “potential” energy requires a bit more work. Using the divergence theorem and the continuity equation one can show that we need


(c.c. represents the complex conjugate). The surface term vanishes if at the surface of the star. It also vanishes for modes that have no radial component. Adding the contributions for the neutron and proton fluids we see that we need


Again, one can argue that the surface term vanishes.

Finally, we have terms of form

It is natural to express these in terms of the perturbed densities using (this is valid for co-rotating background models only)


Adding all the terms together we find that the “potential energy” follows from




With these definitions it follows that the total “energy” is conserved, i.e.


when . These energy expressions are equivalent to those used by Lindblom and Mendell lm95 .

iv.2 Mutual friction

Even though we will include the mutual friction terms in the equations of motion, it is useful to work out the corresponding dissipation integrals. After all, this is the way that the mutual friction damping has traditionally been evaluated lm95 ; lm00 ; yl1 and we want to be able to compare the two approaches.

First consider the terms. It is easy to show that these terms are not dissipative. We find that


by symmetry. This result is not surprising. In fact, we see from (13) that the terms enter the equations of motion in the same way as the Coriolis force. Since the Coriolis terms vanish identically when we multiply each Euler equation with this should be true also for the non-dissipative part of the mutual friction.

Finally, it is straightforward to show that the dissipative terms lead to


Let us now ask how we can use these results to estimate the mutual friction damping timescale. Let us assume that we have a mode solution to the full dissipative problem. That is, we have a solution with time dependence where such that is the damping timescale. From the fact that the energy is quadratic in the perturbations it follows that il91


Moreover, since the solution satisfies the dissipative equations of motion we also know that


Hence, we can equally well use


As long as we are using the complete solution to evaluate this expression, it is an identity. However, in many cases we do not have access to the solution to the dissipative problem. (If we did, we would not need the energy integrals in the first place.) In these cases we can still estimate the damping timescale by evaluating the right-hand side of (36) using the non-dissipative mode solution. When the damping is sufficiently slow, in the sense that the dissipative terms have a small effect on the eigenfunctions, this estimate should be reliable. Of course, one should not expect it to yield exactly the same result as the solution to the full dissipative problem.

iv.3 Gravitational-wave emission

Finally, let us work out the multipole formulas for gravitational-wave emission from a two-fluid star. This exercise is particularly relevant if we are interested in oscillations that may be driven unstable by gravitational-wave emission nareview . The main motivation for including it here is that it demonstrates the intuitive result that gravitational waves are only generated by the co-moving degree of freedom.

Following thorne we need the mass multipoles


and the current multipoles


In these expressions are the standard spherical harmonics and are the magnetic multipoles thorne .

To work these out we start with the usual expression for the two-fluid stress-energy tensor in general relativity lin


In the relativistic formulation, see livrev for a review and a survey of the literature, the central variables are the fluxes . The associated momenta follows from


Hence, the coefficients encode the entrainment effect. Let us now work in the frame of an observer moving with four-velocity such that


where is the associated three-velocity. Then it follows that


We want the Newtonian (low-velocity) limit of this expression. Thus we let


and we get


Finally use the definition


and to arrive at


as one would have expected.

For the current multipoles we need


In the low-velocity limit, this leads to


Finally, we need to write this expression in terms of the Newtonian variables. This can be done by comparing the momenta,


This suggests that we should identify


Using analogous expressions for the protons we see that


which leads to (for a co-rotating background)


These results show that it is only the co-moving degree of freedom that radiates gravitationally.

V Slow rotation perturbation equations

Let us now return to the problem of oscillating superfluid neutron stars. We will first derive the general perturbation equations for a slowly rotating superfluid star. To do this we expand all variables in spherical harmonics. Since we expect rotation to couple the various multipoles, we represent the velocity perturbations by the general expressions




Note that we represent the “co-moving” degree of freedom by the uppercase amplitudes while the “counter-moving” degree of freedom corresponds to the lowercase quantities . All scalar perturbations are expanded in spherical harmonics, i.e. we have etcetera. From now on the sum over will be implied.

One can use a number of different strategies in writing down the perturbation equations. To some extent this is a matter of taste. However, in the slow-rotation problem it can be advantageous to work with a set of equations where the coupling between different multipoles is minimal. The set of equations that we use was chosen using this criterion. We also decided to use the velocity perturbations as our main variables. This approach is analogous to that used by Lockitch and Friedman lockitch in their analysis of inertial modes of single fluid stars. It is notably different from the two-potential formalism pioneered by Ipser and Lindblom ipser , which was extended to superfluid stars by Lindblom and Mendell lm00 .

We replace each of the perturbed Euler equations with three equations. The first is the radial component of the vorticity equation that follows if we take the curl of (8) or (13). Assuming that the perturbations have a harmonic dependence on time, , we get




In deriving these equation we have made use of the standard recurrence relations






For future reference, note that and .

Next we could use also the (or ) components of the vorticity equation. However, as discussed in ekman there is a slightly simpler alternative. We first of all use a pair of equations analogous to the “divergence” equation in ekman . These can be written




Meanwhile the radial components of the Euler equations lead to




Finally, the continuity equations become




This completes the description of the general first order slow-rotation problem.

In order to deduce the relevant recurrence relations from the above equations we need to recall that we have been implying summation over . That is, we are considering relations of form


Using orthogonality of the spherical harmonics, i.e. multiplying by and integrating over the sphere, we obtain the recurrence relation


Given this result, it is straightforward to write down recurrence relations for the various classes of oscillation modes of a rotating superfluid star. However, since the level of rotational coupling is different for different kinds of modes, it is not particularly useful to write down the general relations. Instead, we focus on two specific examples.

Vi The f-modes

Let us begin by considering modes that are non-trivial already in a non-rotating star. Then we first need to solve the non-rotating (and non-dissipative since the mutual friction damping requires rotation) problem. Simply setting in our perturbation equations we see that the polar and axial degrees of freedom decouple (as they should). It is also clear, cf. (55) and (56), that there will not exist any purely axial modes in the non-rotating case. This means that we can make the Ansatz


together with




and similarly for the various scalar perturbation quantities. For example, in the case of the proton fraction we have with


vi.1 The non-rotating problem

At the non-rotating level the equations in Section IV provide the following relations


Meanwhile the continuity equations lead to




Before we proceed, we will simplify the problem. Our aim is to determine analytic approximations for the fundamental modes of the system, including the mutual friction damping. Solving the problem numerically is, of course, straightforward but does not lead to the same level of insight into the dependence on the various parameters. To facilitate an analytic solution, we will combine an incompressible background model with compressible perturbations. This simplifies the calculations considerably. In addition, since this is the same model that was considered by Lindblom and Mendell lm95 we can compare our final results directly to the available literature. We thus assume that and are both constant, while and are not.

It is also useful to introduce a new variable for the co-moving degree of freedom. Let us define


For a single barotropic fluid, corresponds to the perturbed enthalpy. For a compressible background model we would have


However, for the uniform density model the gradient of the proton fraction vanishes so we simply have


It is also worth noting that has the same dimension as .

We now find that (72) and (73) can be written




Before we proceed, we need to decide what variables we want to work with. We can either remove or (or some other combination of these variables) from the problem using thermodynamic identities. Opting for the latter possibility, we use




In these relations we have defined, first of all, the speed of sound as


We have also introduced




and made use of the identity 333Ultimately, this relation follows from the fact that the partial derivatives with respect to the number densities commute, so the mixed second derivatives of the energy functional (the “equation of state”) must be equal.