# r-modes and mutual friction in rapidly rotating superfluid neutron stars

## Abstract

We develop a new perturbative framework for studying the r-modes of rotating superfluid neutron stars. Our analysis accounts for the centrifugal deformation of the star, and considers the two-fluid dynamics at linear order in the perturbed velocities. Our main focus is on a simple model system where the total density profile is that of an polytrope. We derive a partially analytic solution for the superfluid analogue of the classical r-mode. This solution is used to analyse the relevance of the vortex mediated mutual friction damping, confirming that this dissipation mechanism is unlikely to suppress the gravitational-wave driven instability in rapidly spinning superfluid neutron stars. Our calculation of the superfluid r-modes is significantly simpler than previous approaches, because it decouples the r-mode from all other inertial modes of the system. This leads to the results being clearer, but it also means that we cannot comment on the relevance of potential avoided crossings (and associated “resonances”) that may occur for particular parameter values. Our analysis of the mutual friction damping differs from previous studies in two important ways. Firstly, we incorporate realistic pairing gaps which means that the regions of superfluidity in the star’s core vary with temperature. Secondly, we allow the mutual friction parameters to take the whole range of permissible values rather than focussing on a particular mechanism. Thus, we consider not only the weak drag regime, but also the strong drag regime where the fluid dynamics is significantly different.

## 1 Introduction

Oscillations of rapidly rotating neutron stars have attracted interest for a considerable time. During the last decade significant effort was aimed at understanding whether gravitational-wave emission sets the upper speed limit for pulsars, e.g. via the r-mode instability (Andersson, 1998). This possibility is of particular interest since such unstable systems may radiate detectable gravitational waves. In particular, it was suggested that r-modes in rapidly rotating neutron stars in low-mass X-ray binaries (LMXBs) may lead to a persistent source of gravitational radiation (Bildsten, 1998; Andersson et al., 1999b) that may be detected by advanced LIGO [see Watts et al. (2008) for the most recent analysis of this problem]. It is of great importance to understand whether internal fluid dissipation allows the instability to develop in such systems or whether it suppresses the r-modes completely.

There are, in fact, many mechanisms at work in a real neutron star which compete with the gravitational-wave driving of the r-mode. One can express the relative strength of these mechanisms in terms of timescales which will, in general, depend on a large number of parameters. For a basic instability analysis, the most important parameters are the rotation rate and the temperature of the star. The instability can only develop when the gravitational radiation growth time-scale is shorter than the damping timescales due to the various viscosity mechanisms. This defines a region in the spin-temperature parameter space where the r-mode instability is active. Several studies have been devoted to calculating the shear- and bulk viscosity timescales for more or less realistic neutron star models. These studies tend to agree that the notion of persistent gravitational-wave emission accreting neutron stars in LMXBs is quite robust [see Nayyar & Owen (2006) and references therein].

The neutron stars in LMXBs are believed to be old cold recycled pulsars that have been spun up by accretion. These stars are expected to have superfluid cores and are thus expected to have rather different fluid dynamics and dissipation mechanisms. The simplest description of these systems is provided by a two-fluid model, where the superfluid neutrons are distinguished from the proton-electron plasma (usually treated as a charge neutral fluid). The oscillations of neutron stars with such superfluid cores have been studied by Lindblom & Mendell (1994) and Andersson & Comer (2001). The latter showed that there are no g-modes in such stars, but on the other hand one has a new class of “superfluid” modes. In particular, there are two families of r-modes (Prix et al., 2004). One of these resembles the ordinary r-modes of a barotropic star in the sense that the neutrons and protons “move together”. For the other set of modes the neutrons and protons are “counter-moving”. The relative motion of neutrons and protons in a superfluid is, however, rapidly damped by mutual friction, a mechanism that is directly associated with the rotational neutron vortices. The most commonly considered mutual friction mechanism is the scattering of electrons off the magnetic fields associated with the superfluid neutron vortices (Alpar et al., 1984; Andersson et al., 2006). Mutual friction is known to have a decisive impact on neutron star dynamics. In fact, Lindblom & Mendell (1995) have argued that mutual friction completely suppresses the gravitational radiation driven f-mode instability [see Andersson et al. (2008) for a recent discussion]. It is thus important to understand if the r-mode instability suffers a similar fate. Previous work on the subject by Lindblom & Mendell (2000) and Lee & Yoshida (2003) presents a somewhat mixed picture. While the ordinary r-mode is very weakly damped for most values of the entrainment parameter, the damping becomes very strong for a small range of parameter values. Hence, there may be situations where the mutual friction is strong enough to stop the instability from growing.

The purpose of this paper is to clarify the role of the mutual friction damping and understand whether it can stabilize the r-modes against the gravitational-wave instability. In order to do this we extend the formulation of Andersson et al. (2008) to second order in rotation to calculate the mutual damping timescale for r-modes. Although we do not expect our analysis to lead to results that change previous conclusions, there are good reasons to return to this problem. First of all, previous studies have focussed entirely on the weak drag regime for the mutual friction. Meanwhile, recent studies of superfluid turbulence (Andersson et al., 2007; Peralta et al., 2005, 2006), neutron star free precession (Glampedakis et al., 2008a, b) and pulsar glitches (Glampedakis & Andersson, 2008) demonstrate that systems in the strong drag regime may have very dynamics. Since one can make convincing arguments for the strong drag regime being relevant (Ruderman et al., 1998; Link, 2003, 2006; Sedrakian & Sedrakian, 1995), we clearly need to understand the r-mode problem for this parameter range as well. Secondly, we want to improve the treatment of the critical temperature/density at which superfluidity comes into play. The density dependence of the various superfluid pairing gaps translates into distinct regions of normal- and super-fluids that vary with the core temperature. So far, the best analysis in this respect is that of Lindblom & Mendell (2000) who consider a two-fluid core surrounded by a single fluid envelope. We will build on this by considering superfluid layers, the size of which are determined by a (qualitatively) realistic pairing gap. The size of the superfluid region affects not only the mutual friction but also the shear- and bulk viscosities, and it is relevant to demonstrate how this alters the r-mode instability window. Finally, we want to lay the foundation for more detailed work on exotic neutron star cores. It is well established that exotic cores, dominated by either hyperons or deconfined quarks, may be associated with a very strong bulk viscosity [see, for example, Nayyar & Owen (2006)]. Although it has been suggested that this would lead to the suppression of the r-modes, it is clear that such a conclusion is premature. In reality one would expect superfluidity to play a role. If the hyperons are superfluid the reactions that produce the bulk viscosity are suppressed and the effect on the r-mode instability may not be that severe (Nayyar & Owen, 2006). A similar argument applies to colour superconducting quarks (Alford et al., 1999). However, it may be important to keep in mind that a superfluid system has extra dynamical degrees of freedom. In order to truly understand the dynamics of these exotic phases of matter we need to explore the multi-fluid nature of these systems. The present study paves the way for future work in this direction.

## 2 The two-fluid equations

### 2.1 The unperturbed problem

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

(1) |

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

(2) |

where (), and represents the chemical potential (in the following we assume that ). Moreover, represents the gravitational potential, and the parameter encodes the entrainment effect. The force on the right-hand side of (2) can be used to represent other interactions, including dissipative terms (Andersson & Comer, 2006). We will focus on the vortex-mediated mutual friction force for a system that, in equilibrium, rotates uniformly. This means that we consider a force of form (Andersson et al., 2006)

(3) |

Here, is the angular frequency of the neutron fluid (a hat represents a unit vector).

One can express the mutual friction force in terms of a dimensionless ”drag” parameter such that (Andersson et al., 2006)

(4) |

In the standard picture, the mutual friction is due to the scattering of electrons off of an array of neutron vortices (Alpar et al., 1984) which leads to , thus placing the problem in the weak drag regime. This means that and . Hence, all terms can be neglected. However, it is commonly thought (Ruderman et al., 1998; Link, 2003, 2006) that the strong drag limit, , may apply. Intermediate values for the drag, are also interesting. In particular since one would then have both and (note that in the case one still has but ) . This leads to the presence of terms that may have significant effect on the dynamics of the system (Glampedakis et al., 2008a, b; Glampedakis & Andersson, 2008). At this point, our understanding of neutron star core physics is sufficiently rudimentary that we should avoid ruling out the various possibilities. Hence, we will consider the entire permissible range of values for the drag parameter.

Let us consider a frame rotating with the star at fixed angular velocity . The equations of motion then take the form:

(5) |

where we have included the centrifugal term in the potential

(6) |

The continuity equations maintain the form (1) and the Poisson equation is

(7) |

### 2.2 Perturbations

To keep the problem tractable we will assume that the background configuration is such that the two fluids rotate together [see Prix et al. (2004); Glampedakis & Andersson (2008) for discussions of oscillations in systems that are not in co-rotation]. Perturbing the equations of motion and working in a frame rotating with we then have

(8) |

and

(9) |

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. We shall thus assume that the superfluid is only present in a distinct region determined by the core temperature and the superfluid pairing gap. We will discuss the implementation of this idea once we have set up the relevant system of equations.

From previous work on superfluid neutron star oscillations, e.g. Lindblom & Mendell (1994); Andersson & Comer (2001); Prix & Rieutord (2002); Andersson et al. (2008), we know that the problem has two ”natural” degrees of freedom. One of them represents the total mass flux. Introducing

(10) |

and combining the two Euler equations accordingly, we get

(11) |

where and the pressure is obtained from

(12) |

We also have

(13) |

At this point we have two equations which are identical to the perturbation equations for a single fluid system. Of course, we are considering a two-fluid problem and there is a second dynamical degree of freedom. To describe this, it is natural to consider the difference in velocity. Thus, we introduce

(14) |

Combining the two Euler equations in the relevant way we have

(15) |

We have defined

(16) |

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

(17) |

and we remind the reader that

(18) |

Finally, we need a second ”continuity” equation. To close the system, it seems natural to consider an equation for the proton fraction . We then find that

(19) |

### 2.3 Model equation of state

Let us now consider the set of equations that we have written down. We see that the two degrees of freedom and only couple “directly” through (19). For compressible models, the two degrees of freedom also couple “indirectly”, since we can use the equation of state to relate to .

Deciding to work with and [cf. Andersson et al. (2008)], we use

(20) |

and

(21) |

As a simple model, we shall consider an equation of state such that

(22) |

where is the background sound speed. We combine this with a simple linear expression for the proton fraction;

(23) |

This leads to

(24) |

where we have used the relation (Andersson & Comer, 2001)

(25) |

to get

(26) |

In our model calculations we shall take where g/cm is the nuclear saturation density.

In a superfluid system the momentum of each component may not be parallel to it’s velocity. Rather, it acquires a component along the relative velocity. This is evident from equation (2). This non-dissipative coupling is usually parametrised in terms of the entrainment parameter . One can show (Prix et al., 2004) that this parameter is linked to the effective proton mass, , according to

(27) |

For simplicity, we shall take to be constant throughout the superfluid region of the star. Recent work suggests that, while this may not be a good approximation for an entire neutron star, it is approximately true if we consider a shell. From Chamel (2008) we see that a reasonable range for the entrainment parameter is

In order to compare our results to previous work, it is worth recalling that the entrainment parameter used by Lindblom & Mendell (2000) is related to by

(28) |

Hence, their range of corresponds, if we take a volume average of , to .

To summarise; When we solve the r-mode problem, we will consider an equation of state represented by i) the overall density profile, represented by the sound speed , ii) the proton fraction and iii) the entrainment parameter . These quantities allow us to specify the background needed for the perturbation problem. This model is quite simplistic, and it would not be difficult to make it more realistic. However, our main interest is to explore how the r-mode results depend on the different parameters. This question is easier to address with a simple parametrised model.

## 3 Slow rotation analysis

In order to determine the rotational corrections to the superfluid r-modes, we will apply the formalism of Saio (1982) to the problem. Thus, we consider corrections up to order in the slow rotation approximation. We will also make the Cowling approximation, i.e. neglect perturbations of the gravitational potential, .

We start by considering a slowly, and uniformly, rotating star. It is well known that such a star will not be spherical and that a distorted potential surface can be written in the form

(29) |

Here is a function of and which represents the deformation of the equilibrium structure from the background spherical state. The equilibrium physical quantities are functions only of . Following Saio (1982); Smeyers & Martens (1983), we write the equations of motion in a frame corotating with the star, denoted by , starting from a static Cartesian frame, denoted by . Our coordinates in the rotating frame will be a spherical polar system, explicitly:

(30) | |||||

(31) | |||||

(32) |

The metric in the new coordinates is then given by

(33) |

which, after linearizing with respect to , leads to

(34) |

The connection coefficients, in the coordinate basis

(35) |

can be obtained from

(36) |

We shall, however, use the vector basis

(37) |

for which the covector basis is

(38) |

Note that the basis vectors are not orthogonal.

Before we proceed, it is worth recalling that Smeyers & Martens (1983) have argued that the derivation of Saio’s results is flawed, even though the final results for the rotational frequency correction are correct. Since this is an important point, we will outline how the second order slow-rotation perturbation equations should be derived following Saio’s strategy.

As we are considering linear perturbations on a rotationally distorted background, the first step is to calculate such a background, using the coordinates defined above. The equations that govern the background equilibrium are particularly simple (if we assume that neutrons and protons do not move relative to each other)

(39) |

where , defined in equation (6), includes the centrifugal terms. This is due to the fact that the coordinate has been defined in such a way as to label the equipotential surfaces of , which are also isobaric (and isopycnic as the background equation of state is barotropic) surfaces of the star, thus leading to all the background quantities being a function of only [for a detailed description see Tassoul (1978)]. The Euler equations are thus identical to those one would obtain for a spherical background and the solution can be formally obtained by replacing the radial coordinate with the variable in the spherically symmetric solution. It is important to remember, however, that we are calculating physical quantities at a point in the deformed star labeled by the coordinates and not at a point in the spherical star labeled by coordinates . Note, in fact, that as the geometry is not spherical the measures of distance and volume change. For example, the mass of the star is now obtained from

(40) |

where is the mass of the spherically symmetric star and is the difference in mass due to rotation.

### 3.1 Perturbations

Let us now consider linear perturbations on our rotationally deformed background. We shall express the perturbations in terms of a total displacement vector , such that , and a difference displacement vector , such that . We can then write equation (11) as:

(41) |

where the that denotes the Eulerian perturbations should not be confused with the Kronecker delta . We are using to represent the Coriolis force, and the vector has the following components:

(42) | |||||

(43) | |||||

(44) |

The differential operator is the same as for spherical polar coordinates, i.e.

(45) |

The continuity equation (13) becomes

(46) |

We shall also need the linearized equation of state. In our model case, we have

(47) |

That is, in comparing our equations to Saio’s results one should only consider the barotropic limit.

## 4 The r-mode problem

We are mainly interested in understanding the effect that mutual friction has on the r-mode instability. It is well-known that the r-modes have zero frequency and are purely axial in the non-rotating limit. Rotation breaks this degeneracy and leads to modes that are purely axial to leading order but which have a poloidal component of order . These modes have frequency

(53) |

The r-mode frequency can be calculated, to second order in rotation, for a single fluid star using a number of different formalisms (Andersson et al., 1999a; Lindblom et al., 1999; Yoshida & Lee, 2000). The superfluid problem has been solved both to linear (Andersson & Comer, 2001; Prix et al., 2004; Andersson et al., 2008) and second (Lindblom & Mendell, 2000; Lee & Yoshida, 2003) order in rotation. To leading order one finds that the ordinary r-mode has frequency

(54) |

where, in fact, only the modes exist. In addition, a constant density model supports a set of purely axial counter-moving modes with frequency (Andersson et al., 2008)

(55) |

Since the frequency in (55) contains the term , it is easy to see that it can only be the solution for a global mode if is constant. In the present analysis we will assume that is constant. Hence, the second class of r-modes can only exist if we also assume that the proton fraction is constant. If is constant, the co-rotating and counter-rotating degrees of freedom are, in fact, completely decoupled also at second order in rotation. The r-mode is then exactly the same as for a barotropic single fluid star. In particular, r-mode solutions exist only for . However, we are not generally considering a constant (in fact we are assuming that the proton fraction scales linearly with the total mass density, which is not constant). Hence, we see there will be no counter-moving r-modes in our model. In fact, if is not a constant, then the coupling between the degrees of freedom leads to the counter-moving r-mode becoming a general inertial mode, with a mixed toroidal/poloidal velocity field to leading order. These modes have been determined numerically by Lee & Yoshida (2003). If we consider a generic profile for and work at second order in rotation the leading order r-mode will drive the counter-moving degrees of freedom, leading to mutual friction damping. This is the effect that we are interested in.

### 4.1 Perturbation equations

So far, we have written down the equations for a general perturbation problem on the rotationally deformed background. We will now focus on modes that are purely toroidal to leading order. That is, we write the total displacement vector as

(56) |

Analogously, the difference displacement vector can be written

(57) |

Note that we use uppercase variables for the co-moving degree of freedom and the corresponding lowercase variable for the counter-moving degree of freedom.

As noted by Smeyers & Martens (1983), the above form for the displacement vectors does not correspond to the usual decomposition in spheroidal and toroidal components. At second order in rotation, the above basis differs from the standard basis on a spherical background as the vector is no longer orthogonal to . This is, however, not a problem as we do not need to explicitly identify the poloidal and toroidal parts of the mode at second order in rotation. Such an identification would not be very useful anyway since the components are coupled at .

We want to study the classical r-mode, i.e. a mode that to first order in rotation involves only the co-moving degree of freedom and which is purely toroidal in the slow-rotation limit. This identification is unique since our decomposition coincides with the standard one into poloidal and toroidal components at . The requirement is that the leading term is of order unity while the amplitudes of the ”spheroidal” components are of order for the total displacement. Meanwhile, following the first order results of Andersson et al. (2008), all components of the difference displacement must be of order . As discussed in Appendix A, this ordering for the displacement vector, together with the frequency from equation (53), leads to the following equations for the r-modes;

(58) | |||||

(59) | |||||

(60) | |||||

(61) | |||||

(62) | |||||

(63) | |||||

where

(65) | |||||

(66) |

In the last two expressions we have used

(67) | |||||

(68) |

It is easy to verify that these equations coincide with the results of Saio (1982) in the barotropic limit. In fact, in order to facilitate this comparison we have used essentially the same notation as Saio. Thus, we have defined the normalised mode frequency

(69) |

and the following background quantities;

(70) | |||||

(71) | |||||

(72) | |||||

(73) | |||||

(74) | |||||

(75) | |||||

(76) | |||||

(77) |

where is the Legendre polynomial:

(78) |

We are also using the definition

(79) |

For later convenience, it is worth noting that this means that .

### 4.2 Rotating polytropes

We will now restrict ourselves to the case where the density profile of the equilibrium configuration is that of an polytrope. This greatly simplifies the analysis, as the background quantities we are interested in can be obtained in closed form. Explicitly, we have

(80) |

with

(81) |

Introducing the dimensionless variable we then find that

(82) |

from which we can define

(83) |

Furthermore, we can calculate, following the classical work of Chandrasekhar (1933), the rotationally induced deformation

(84) |

and

(85) |

where

(86) | |||||

(87) |

The mass of the star is thus, from equation (40), given by

(88) |

In the following, when we quote the mass of the star we shall in fact be refering to the mass of the spherical star, . Formally what we are doing is thus considering a sequence of stellar models with the same central density, but with a mass that varies with the rotation rate. The difference in mass does not, however, enter the r-mode frequency correction to order and is thus irrelevant for the present discussion.

Finally, we will also need the sound speed, which simply follows from

(89) |

### 4.3 Decoupling the degrees of freedom

Let us now examine the consequences of the ordering we assumed at the beginning of the analysis. From equation (62) we see that if is of order unity then (i.e. the pressure perturbation) must be of order . Similarly, from equation (63) we see that, as all components of the difference displacement are of order , the variable (i.e. ) must be of order . It follows that in equation (59) the term involving is of higher order than the others and can be neglected. The consequence of this is that the equations for the co-moving degree of freedom are completely decoupled. Hence, they can be solved independently and then used as source terms for the counter-moving degrees of freedom. This approach has obvious advantages compared to previous fully numerical calculations (Lindblom & Mendell, 2000; Lee & Yoshida, 2003). After all, we can now solve for the r-mode throughout the star, imposing regularity at the centre of the star and the vanishing of the Lagrangian perturbation of the pressure, , at the surface. As a second independent step, we solve the equations for the counter-moving degree of freedom. Adopting this strategy it becomes straightforward to account for, for example, the temperature dependence of the superfluid gap and the associated variation of the superfluid region in the star.

The solution for the comoving degree of freedom is particularly simple. Not only is it decoupled from the counter-moving degree of freedom, the equation for also decouples and takes the form:

(90) |

leading to the solution

(91) |

One can use this solution to determine from equation (58). This leads to

(92) |

where and (for the modes we are considering)

(93) |

with defined in equation (154). To solve equation (92) let us first of all consider the solution to the homogeneous problem. This is readily obtained and takes the form

(94) |

Note that this solution diverges both at the centre and the surface of the star. We next need to determine a particular solution to the problem. To do this let us write as

(95) |

where

(96) | |||||

(97) |

Inspired by the solution to the homogeneous problem we make the following ansatz for the particular solution to (92);

(98) |

This leads to

(99) |

The first term is thus of the form

(100) |

which is easily integrated. For the second part we need to be able to integrate, for the term proportional to ,

(101) |

so that, for an polytrope we require

(102) |

For the remaining part, involving , we require the integral

(103) |

which, using the explicit form for in equation (85), gives

(104) |

where

(105) | |||||

(106) |

This analysis is perhaps a little bit messy, but the result is very useful. Collecting the various results we can write the final solution to equation (99) as

(107) |

That is, we have an analytic solution to the problem.

As we now have the full solution to the problem we can impose boundary conditions. First of all we require the solution to be regular at the centre of the star. This determines the constant in the homogeneous solution. The remaining condition is that the Lagrangian variation of the pressure vanish at the surface of the star, which in terms of our variables means that