Hydrodynamics of rapidly rotating superfluid neutron stars with mutual friction
Abstract
We study time evolutions of superfluid neutron stars, focussing on the nature of the oscillation spectrum, the effect of mutual friction force on the oscillations and the hydrodynamical spinup phase of pulsar glitches. We linearise the dynamical equations of a Newtonian twofluid model for rapidly rotating backgrounds. In the axisymmetric equilibrium configurations, the two fluid components corotate and are in equilibrium. We use analytical equations of state that generate stratified and nonstratified stellar models, which enable us to study the coupling between the dynamical degrees of freedom of the system. By means of time evolutions of the linearised dynamical equations, we determine the spectrum of axisymmetric and nonaxisymmetric oscillation modes, accounting for the contribution of the gravitational potential perturbations, i.e. without adopting the Cowling approximation. We study the mutual friction damping of the superfluid oscillations and consider the effects of the nondissipative part of the mutual friction force on the mode frequencies. We also provide technical details and relevant tests for the hydrodynamical model of pulsar glitches discussed by ?2010MNRAS.tmp..554S. In particular, we describe the method used to generate the initial data that mimic the preglitch state, and derive the equations that are used to extract the gravitationalwave signal.
keywords:
methods: numerical – stars: neutron – stars: oscillation – star:rotation – gravitational waves1 Introduction
Mature neutron stars are expected to have superfluid and superconducting components in their interior. Shortly after a neutron star’s birth the temperature decreases below , at which point superfluid neutrons should be present both in the inner crust and the outer core, while the core protons should form a superconductor. At all relevant temperatures, the electrons form a “normal” fluid that is tightly locked to the protons due to the electromagnetic interaction. This suggests that the dynamics of mature neutron stars depends on the detailed interaction between coupled superfluidssuperconductors (Glampedakis, Andersson & Samuelsson, 2010), i.e. represents a complex physics problem. The situation is not expected to simplify if one also accounts for the inner neutron star core, at several times the nuclear saturation density, where exotic states like hyperon superfluid mixtures or deconfined quark condensates may be present.
Although it is generally appreciated that neutron stars have this very complicated structure, the evidence for the presence of the different superfluid phases remain indirect. The strongest support comes from observed pulsar glitches, rapid spinup events seen in a number of young pulsars (and also some magnetars) during their magnetic slowdown phase. The typical glitch size is very small, representing a relative change () in the observed rotation rate () in the range . The currently accepted model for these events relies on the transfer of angular momentum between a (faster spinning) superfluid neutron component and the star’s (slower spinning) elastic crust (to which the magnetic field is anchored). The exchange is thought to be mediated by neutron vortices (by means of which the superfluid mimics bulk rotation) and the associated mutual friction (Alpar et al., 1984).
A challenge for future observations is to probe the detailed physics of a neutron star’s interior. In this context, asteroseismology associated with either gravitational or electromagnetic signals seems particularly promising. In fact, the quasiperiodic oscillations seen in the tails of giant magnetar flares may have provided us with the first opportunity to test our theoretical models against observational data (see for instance Watts & Strohmayer, 2007, and references therein). The observed variability likely originates from crustal oscillations and depends on the detailed crust dynamics and the interaction with the neutron star’s magnetic field. These observations have led to a resurgence of interest in neutronstar seismology and a renewed assault on the problem of magnetic star oscillations, a seriously challenging problem from the theory pointofview (see Colaiuda, Beyer & Kokkotas, 2009, for a discussion of the literature). In the context of the present paper, the potential relevance of the neutron superfluid that penetrates the neutron star crust is particularly relevant (Andersson, Glampedakis & Samuelsson, 2009; Samuelsson & Andersson, 2009). The prospect of detecting gravitational waves from oscillating neutron stars is also exciting, especially since the associated signals will allow us to probe the highdensity region and hence the supranuclear equation of state (EoS) (Andersson & Kokkotas, 1998; Benhar, Ferrari & Gualtieri, 2004; Samuelsson & Andersson, 2007; Andersson et al., 2009).
In order to faciliate future observations and the decoding of collected data, we need to improve our models considerably. The superfluid aspects are particularly interesting in this respect, since the oscillation spectrum of a superfluid star is more complex than that of a single fluid model. In superfluid regions fluid elements can execute both co and countermoving motion, leading to the existence of unique “superfluid” oscillation modes. Our understanding of the nature of the additional degree(s) of freedom and the effect on observables must be improved by detailed modelling, ultimately in the context of general relativistic multifluid dynamics.
The present work presents recent progress towards this goal. We study the oscillations of superfluid neutron stars by evolving in time the linearized twofluid equations in Newtonian gravity. We improve on the analysis of Passamonti et al. (2009a) by including the perturbations of the gravitational potential. We also account for the mutual friction force associated with vortices, and implement quadrupole extraction of the gravitationalwave signal associated with the fluid motion. We provide the detailed analysis (and relevant code tests) for the configurations that we recently used to study the hydrodynamics of pulsar glitches (Sidery et al., 2010). We consider two simple analytical EoS and construct two distinct sequences of rapidly rotating stars, the main difference being the presence or absence of composition gradients. Such gradients impact on the superfluid dynamics, as the co and countermoving degrees of freedom are coupled in stratified models. From timeevolutions of the relevant perturbation equations, with the gravitational potential perturbation included, we determine the axi and nonaxisymmetric oscillation modes for models that rotate up to the mass shedding limit. Finally, we account for the (standard form of the) mutual friction force. This adds two coupling terms to the equations of motion. One component is dissipative and damps an oscillation mode, while the other modifies the frequencies of the superfluid modes. We study both these effects and infer an analytical relation for the associated frequency change of the nonaxisymmetric superfluid fundamental and inertial modes.
2 Equations of Motion
In a basic model for superfluid neutron stars, the matter constituents are superfluid neutrons, superconducting protons and normal electrons. Given the typical dynamical timescale of stellar oscillations, one would expect the charged particles to be efficiently locked together by the electromagnetic interaction. Therefore, the dynamics of superfluid stars depends on two components, a neutron superfluid and a neutral conglomerate of protons and electrons. For simplicity, we will refer to the latter mixture as the “protons” in the following. More detailed discussion and justification for the twofluid model is provided by Mendell (1991a, b), Prix (2004) and Andersson & Comer (2006).
When the mass of each fluid component is conserved, i.e. when we neglect the various particle reactions, the dynamics of a superfluid star is described by two mass conservation laws, two Eulertype equations and the Poisson equation for the gravitational potential (Prix, 2004). These take the form;
(1) 
(2) 
(3) 
These equations are given in a coordinate basis, which means that the indices and denote spatial components of the various vectors. Meanwhile the indices and label the two fluid components. In the present case these constituent indices will be for the neutrons and for the protons. Throughout this work, the summation rule for repeated indices applies only for spatial indices. In equations (1)–(3), the total mass density is , is the chemical potential for each fluid component (scaled with the particle mass ), is the gravitational potential, while the relative velocity between the two fluids is . The parameter accounts for the nondissipative entrainment effect. In a neutron star core the entrainment is due to the strong interaction between the nucleons. From equation (2), it is clear that it leads to a momentum that is not longer aligned with the individual component velocity. The vector field represents the force density acting on the fluid component. In this paper we consider only the vortex mediated mutual friction force. The general form of this force is
(4) 
where represents the bulk rotation (later we will assume that the two fluids corotate in the unperturbed background), and and are the mutual friction parameters.
2.1 Equation of State
The equation of state (EoS), that is needed to close the system of equations, can be described by an energy functional
(5) 
that ensures Galilean invariance. The chemical potential and the entrainment parameter are then defined by
(6)  
(7) 
When the relative velocity between the two fluids is small, as is the case in most systems of practical relevance, equation (5) can be expanded in a series:
(8) 
This has the advantage that the bulk EoS and the entrainment parameter can be independently specified at . From equation (7) it follows that the entrainment parameter is related to the function by
(9) 
Despite recent developments (Chamel, 2008), we do not yet have a realistic EoS that consistently describes the superfluid properties of a neutron star core. Therefore, we consider two analytical EoS, based on generalisation of the familiar polytrope. These models are particularly useful if we want to explore the role of entrainment, composition stratification and symmetry energy. Moreover, since the two EoS have been used elsewhere we have “independent” tests of our numerical results. The main difference between our two sets of models is the presence, or absence, of composition gradients. This is important since the co and countermoving degrees of freedom are coupled in stratified neutron stars, which means that the gravitationalwave spectrum may contain the imprints of “superfluid” modes (see Sec. 6). This would not be the case in a nonstratified model.
The first EoS is determined by the following expression (Prix et al., 2002; Yoshida & Eriguchi, 2004; Passamonti et al., 2009a):
(10) 
where is a polytropic constant, is the proton fraction and is a parameter that can be related to the symmetry energy (Prix et al., 2002). In this EoS, both and are taken to be constant (Passamonti et al., 2009a). Using equation (10), we construct a sequence of corotating axisymmetric configurations without composition gradients. These correspond to the A models used by Passamonti et al. (2009a).
In order to study the effects of stratification on the oscillation spectrum we consider a second EoS, defined by (Prix & Rieutord, 2002; Andersson et al., 2002; Passamonti et al., 2009a):
(11) 
Here, the coefficients and are constants. We consider , for all the rotating models. In the numerical code, the coefficients are given in units of , where is the gravitational constant and is the equatorial radius of the stellar model. We take them to have the values and for the nonrotating model, which corresponds to model III used by Prix & Rieutord (2002). Note that for rotating models, the dimensionless can assume a different value with respect to the nonrotating star. For instance, when we impose that the central proton fraction is constant for all the sequence of rotating models (see Section 3.1). From equations (6) and (11) it follows that the chemical potential and mass density are related by
(12) 
where the polytropic index is given by . From this result we can determine the proton fraction for a given stellar model by imposing equilibrium. After some calculations, we obtain:
(13) 
This EoS will be used to construct a sequence of stratified rotating models, where the central proton fraction is fixed, . These equilibrium configurations have already been used by Sidery et al. (2010), and will be refered to as models C in the following.
3 Equilibrium configurations
We study the oscillations of rotating axisymmetric background models where neutrons and protons are in equilibirum and corotate with constant angular velocity, i.e. we have . In this work, we also assume that two fluid components coexist throughout the stellar volume. This is obviously artifical; the outer region of a real neutron star will not be superfluid. However, at this stage our main interest is in the bulk core dynamics. In the future we plan to extend our model to account appropriately for the expected superfluid regions. At that point we will also consider the role of the elastic crust.
In Sec. 3.1, we introduce the equations that govern stationary corotating equilibrium configurations and solve them numerically for the EoS (10) and (11). In Sec. 3.2, we then describe the perturbative approach developed by Yoshida & Eriguchi (2004) for determining stationary configurations in which the two fluids rotate with a small velocity lag. With this method we can obtain noncorotating models as small deviations from a corotating equilibrium. Subsequently, we use this approach to determine the initial conditions for hydrodynamical glitch evolutions.
Model  

A0  1.00000  0.00000  0.00000  0.00000  1.2732  1.2732 
A1  0.99792  0.05913  0.08153  0.05802  1.2701  1.2701 
A2  0.98333  0.16675  0.22992  0.38482  1.2479  1.2477 
A3  0.95000  0.28729  0.39613  1.16918  1.1967  1.1962 
A4  0.90000  0.40268  0.55524  2.38295  1.1186  1.1179 
A5  0.80000  0.55626  0.76700  4.93320  0.9557  0.9576 
A6  0.70000  0.65789  0.90713  7.56798  0.7801  0.7917 
A7  0.60000  0.71733  0.98909  9.86465  0.5794  0.6176 
A8  0.55625  0.72524  1.00000  10.2760  0.4749  0.5361 
C0  1.00000  0.00000  0.00000  0.00000  1.0826  1.1755 
C1  0.99792  0.55856  0.08403  0.04561  1.0798  1.1725 
C2  0.98333  0.15764  0.23716  0.36682  1.0601  1.1516 
C3  0.95000  0.27145  0.40837  1.11303  1.0146  1.1034 
C4  0.90000  0.38006  0.57177  2.26285  0.9447  1.0301 
C5  0.80000  0.52334  0.78732  4.64841  0.7984  0.8792 
C6  0.70000  0.61536  0.92576  7.01965  0.6392  0.7223 
C7  0.60000  0.66249  0.99667  8.77258  0.4562  0.5559 
C8  0.57656  0.66471  1.00000  8.87526  0.4077  0.5147 
3.1 Corotating background
The equations that describe rapidly and uniformly rotating background models can be derived by imposing the conditions of stationarity and axisymmetry on the Eulertype equations (2) and the Poisson equation (3) (Prix et al., 2002; Yoshida & Eriguchi, 2004). This leads to
(14)  
(15) 
where and are, respectively, the angular velocities and the integration constants for the neutron and proton fluids. For corotating background models, i.e. when , in which the two fluids are in equilibrium and share a common surface, the hydrostatic equilibrium equation (14) becomes
(16) 
where is the background chemical potential and . When the system of equations (15)–(16) is closed by an EoS, we can numerically determine a corotating stationary axisymmetric background via the selfconsistent field method (Hachisu, 1986; Passamonti et al., 2009a). The solution is such that the surface of the star corresponds to the zero chemical potential surface, (Yoshida & Eriguchi, 2004).
For the EoS (10) and (11) we construct the sequences of rotating models A and C, respectively. The set of models extends from a nonrotating model up to the mass shedding limit. In the numerical code, we rewrite the background equation in dimensionless form by using the gravitational constant , the central mass density and the equatorial radius . All stellar models of the sequence have the central proton fraction set to . By specifying the axis ratio between the polar and equatorial radius , the iterative numerical routine determines all other quantities of an axisymmetric configuration. The main properties of the rotating models are given in Table 1. From these quantities we can easily construct stellar models in physical units. For instance for models C, we can evaluate equation (12) at the centre and obtain:
(17) 
where the asterisk denotes the dimensionless quantities and . Combining equations (17) with the dimensionless mass , we can derive the equatorial radius of the star,
(18) 
where the physical mass and the EoS parameters can be arbitrarily chosen. The central mass density, , and the rotational period, , are determined by the following equations:
(19)  
(20) 
where .
3.2 Noncorotating solutions
In a multifluid system, like an astrophysical neutron star, the various fluid components can have different velocities. This is, in fact, an essential element in the favoured model for pulsar glitches where the sudden observed spinup is explained as a transfer of angular momentum between an interior superfluid neutrons and the charged component. In this model, the momentum transfer is due to the interaction between the crust and an array of quantised neutron vortices that are generated by the stellar rotation. During the magnetically driven spindown of a neutron star, these vortices are pinned to the crust and corotate with the charged components. Therefore, a velocity lag develops between superfluid neutrons and the crust and an increasing Magnus force acts on the vortices. When this force becomes stronger than the pinning force, the vortices should unpin. At this point they are free to move and can accelerate the crust, generating a glitch.
Typically, the spin variation observed in a glitch is very small, . This means that the effects of a glitch on the stellar structure is expected to be tiny and can be studied perturbatively. The approach developed by Yoshida & Eriguchi (2004) is particularly appropriate for this kind of problem, as the noncorotating quantities are considered as small deviations from a stationary, rapidly corotating configuration. We have already used these noncorotating corrections as initial data for studying the postglitch dynamics and the associated stellar oscillations (Sidery et al., 2010). We will now provide further details about the method.
Adopting the Yoshida & Eriguchi (2004) approach, we expand equations (14)(15) up to the first order in
Thus, we have
(21)  
(22)  
(23)  
(24) 
where the subscript “c” denotes the corotating values. Note that by definition represents the relative deviation of the x fluid angular velocity with respect to the corotating background, i.e. . For the noncorotating corrections, Equations (14)–(15) become
(25)  
(26) 
The system of equations is closed by
(27) 
that relates the mass density and the chemical potential perturbations for corotating backgrounds.
Noncorotating solutions can be constructed with either fixed central chemical potential or total mass. For the first class of models, we can impose the condition at the star’s centre, and determine the integration constant from equation (25) (Yoshida & Eriguchi, 2004):
(28) 
For solutions with constant mass, we impose a constraint on the mass of each fluid component, i.e.
(29) 
In equation (25), we can replace the chemical potential by the mass density perturbation using equation (27), and integrate over the star’s volume . The integration constant is then given by the following expression:
(30) 
For the EoS (10), the adiabatic index is and the boundary condition (30) therefore reduces to:
(31) 
The system of equations (25)–(26) can be solved iteratively. First of all, for a given EoS we determine the corotating background with the selfconsistent field method of Hachisu (1986), where we specify the axis ratio of the star. Secondly, we choose the relative angular velocity of each fluid component. The iteration algorithm then proceeds as follows: i) we solve the perturbed Poisson equation (26) for an initial guess of the perturbed mass density , ii) we get the integration constant imposing either the condition (28) or (30), iii) we determine the chemical potential from equation (25) and then the new mass density from the EoS. This procedure is iterated until the difference between the quantities is smaller than a prescribed error.
An important property of this linear perturbation approach is that we can construct two independent solutions to equations (25)–(26), respectively corresponding to and . Since the problem is linear, any noncorotating configuration can be obtained as a linear combination of these two solutions.
We have tested our code against the analytical solution for the EoS (10) determined by Prix et al. (2002) in the slowrotation approximation. We select two slowly rotating models with axis ratio and , respectively. The models have the same proton fraction and symmetry energy term, i.e. and . The noncorotating corrections correspond to a relative angular velocity with constant central chemical potential, c.f. (28). In Fig. 1, we show the radial profile of the perturbed neutron mass density for the three angles and , respectively. In the slowest rotating model, the agreement between the numerical and the analytical solutions is evident. In the second model, with axis ratio the two solutions begin to differ, as expected. The slowrotation solution becomes less accurate as the star’s rotation increases. The same behaviour is found for noncorotating solutions with constant mass, i.e. when . This comparison gives us confidence in our numerically generated background models.
For the sequence of constant mass models, we show in Figs. 2 and 3 the noncorotating mass density pertubations and the gravitational potential perturbation for the C2 model. These are solutions to equations (25)–(26) with and , which were used as initial conditions for the glitch simulations discussed by Sidery et al. (2010).
4 Perturbation Equations
The dynamics of a superfluid neutron star can be studied by linearizing the system of differential equations (1)–(3). In the inertial frame, the Eulerian perturbation equations are given by
(32)  
(33)  
(34) 
where is the azimuthal angle associated with the rotational motion, and the perturbed mutual friction force is (in the case of a corotating background)
(35) 
The chemical potential perturbations can be expressed in terms of the mass density perturbations using equation (27).
In order to solve numerically equations (32)–(33) we use the conjugate momentum perturbations as dynamical variables. These are given by
(36)  
(37) 
where we recall that . By inverting these relations we can determine the velocity fields at any time step,
(38)  
(39) 
where .
The time evolution of the nonaxisymmetric perturbation equations is a threedimensional problem in space. However, linear perturbations on an axisymmetric background can be expanded in terms of a set of basis functions , where is the azimuthal harmonic index (Papaloizou & Pringle, 1980). The mass density perturbations as well as the other perturbation quantities then take the following form (Jones et al., 2002; Passamonti et al., 2009a)
(40) 
With this Fourier expansion the perturbation equations decouple with respect to and the problem becomes twodimensional. In particular, for the axisymmetric case () only the component survives.
4.1 Boundary Conditions
In this work, we study axisymmetric () and nonaxisymmetric oscillations () of a superfluid neutron star with equatorial and rotational axis symmetry. The numerical domain extends over the region and , and we need to impose boundary conditions at the surface, origin, rotational axis and equator.
We first discuss the boundary conditions at the origin () and the rotational axis (), where the perturbation equations must be regular. Let us denote by a general scalar perturbation, such as the mass density , the chemical potential and the gravitational potential . For axisymmetric and nonaxisymmetric oscillations, we have to impose the following conditions, respectively :
(41)  
(42) 
For the velocity fields , we impose that there must be no mass flux across the origin () for both axisymmetric and nonaxisymmetric perturbations:
(43) 
At the rotational axis (), we impose the following conditions:
(44)  
(45) 
At the equator (), the reflection symmetry divides the perturbations into two sets with opposite parity (Passamonti et al., 2009b). In the Type I parity class, the scalar perturbations and the velocity satisfy the following conditions:
(46) 
Meanwhile, the Type II class is such that:
(47) 
The outer layers of a mature neutron star form an elastic crust made up of nuclei. The crust is an important aspect that is yet to be implemented in our numerical model (although we are making progress on it). Our current model is simplified, in the sense that we assume that superfluid neutrons and protons are present throughout the stellar volume. We then impose the standard boundary condition of a free surface, i.e. require that the Lagrangian perturbation of the individual chemical potentials vanish at the surface, i.e.
(48) 
The vector field is the Lagrangian displacement of the xfluid component (Andersson, Comer & Grosart, 2004). The value of the perturbed chemical potential at the surface is determined from equation (48) at each time step.
5 Gravitationalwave Extraction
In order to study the gravitationalwave signal emitted by pulsating superfluid neutron stars, we have implemented the quadrupole formula for both axisymmetric and nonaxisymmetric oscillations. We will now discuss this implementation, in particular, the momentum and stress formula that we use to improve the numerical gravitationalwave extraction.
The gravitationalwave strain can be determined using the quadrupole formula (Thorne, 1980):
(49) 
where is the pure spin tensor harmonic which has “electrictype” parity, i.e. (Thorne, 1980). In this work, we focus only on the and pulsations. In the orthonormal basis of spherical coordinates, the components of the and spin tensor harmonics are, respectively, given by
(50) 
and
(51)  
(52) 
The quantity is the quadrupole moment, in the case of a twofluid star defined by;
(53) 
where the spherical harmonics for the and cases are given by
(54)  
(55) 
where is the Legendre polynomial.
It is wellknown that, the numerical calculation of the second order time derivative of the quadrupole moment in equation (49) could lead to inaccurate results (Finn & Evans, 1990). However, the accuracy of the gravitationalwave extraction can be improved by transforming equation (49) into either the perturbed momentum formula, with a first order time derivative, or the perturbed stress formula, where the time derivatives are absent (Finn & Evans, 1990). In this work, we use both these prescriptions in order to check the wave extraction accuracy.
For axisymmetric oscillations, , the gravitational strain can be written as follows:
(56) 
where the quantity is defined by
(57) 
We can reduce the order of the time derivative by using the method developed by Finn & Evans (1990), and obtain the perturbed momentum formula:
(58) 
and the perturbed stress formula:
(59)  
where the gradient components in eqaution (59) are determined in the orthonormal spherical basis, i.e. . At the end of the day, the quantity in the strain equation (56) can be determined from either of the three equations (57)–(59).
For nonaxisymmetric oscillations with , the two independent polarizations of the strain can be written as follows:
(60) 
where is the spinweighted spherical harmonics,
(61) 
and we have defined the quantity
(62) 
We can then rewrite equation (62) as follows:
(63) 
where
(64) 
In equation (64), the order of the time derivatives can be reduced by using the equations of motion (see Appendix A for more details). This leads to the following expression:
(65) 
For linear perturbations on a corotating background, we can further transform equation (65) into the following expression:
(66)  
where the time derivatives are absent. In equations (65) and (66), the perturbations are determined in the inertial frame.
The energy radiated as gravitational waves is determined by the following equation (Thorne, 1980):
(67) 
By using Parseval’s Theorem we can write equation (67) for the and components as follows:
(68)  
(69) 
where , and is its Fourier transformation.
The characteristic strain of the gravitationalwave signal is then given by (Flanagan & Hughes, 1998):
(70) 
where is the source distance. The strains and are related to the dimensionless quantities and used in the numerical code by the following expressions:
(71)  
(72) 
Similar relations provide the characteristic strain
(73)  
(74) 
As a first test of the numerical implementation, we compare the gravitationalwave extraction formulae for the axisymmetric and nonaxisymmetric oscillations. We evolve the C2 model with a density perturbation and extract the signal using equations (57)–(59) for the pulsations, and (65)–(66) for the oscillations. Typical results are shown in Figure 4. We generally find good agreement between the different numerical results, although we note that (as expected) the momentum and stress formulae produce a smoother signal than the “raw” quadrupole formula (57).
As an additional test, we have used the relativistic numerical code developed by Nagar & Diaz (2004) and Passamonti et al. (2007) to test the results of the gravitationalwave extraction routine. In the relativistic case, the linear perturbations of nonrotating relativistic stars were evolved and the signal was extracted using the Zerilli function (Zerilli, 1970). From the Newtonian approach used in the current work, it is evident that we cannot accurately reproduce the relativistic results. However, we can establish that our calculations provide a good estimate of the amplitude of the gravitationalwave strain. To this end, we consider a star with mass and radius , and evolve the relativistic code with an initial enthalpy perturbation, which produces an averaged pulsational kinetic energy of , where is the speed of ligth. The related gravitationalwave strain is almost monochromatic and for a source at the maximal amplitude is .
Mode  

PR  [ % ]  
1.91361  1.92743  0.7  
2.52823  2.53376  0.2  
3.94917  3.94911  0.1  
4.20552  4.20420  0.1  
5.61069  5.52870  1.5  
5.93799  5.92165  0.3  

1.33511  1.33178  0.2 
1.83142  1.82281  0.5  
3.47686  3.48786  0.3  
3.68465  3.69878  0.4  
5.24187  5.25802  0.3  
5.51946  5.52876  0.2 
With our 2D Newtonian code, we then evolve in time nonradial oscillations of a superfluid nonrotating star for both the EoS (10) and (11). The kinetic energy of oscillating superfluid stars can be determined by the following expression:
(75) 
If we evolve oscillations that have the same pulsational kinetic energy as in the case studied with the relativistic code, we obtain for model A0 and for model C0. In this calculation, we have used equation (71) with the parameters of the relativistic stellar model. This test shows that we can be confident that the implementation of the quadruple formula in our code provides reasonable results, in accordance with the expected relation between the pulsational kinetic energy and the gravitationalwave strain.
6 Results
Having formulated the timeevolution problem and described our implementation of the gravitationalwave extraction, we will now discuss our results. In this section, we focus on the effects of the gravitational potential perturbation and the mutual friction force on axisymmetric and nonaxisymmetric oscillations. We also provide a more detailed analysis of the gravitationalwave signal generated by the basic glitch model that we discussed in a previous work (Sidery et al., 2010).
The pulsation dynamics is studied with a numerical code that evolves in time the system of hyperbolic perturbation equations (32)–(33), solving at each time step the perturbed Poisson equation (34). The part of the code that evolves the hyperbolic equations uses the same technology as in previous work (Passamonti et al., 2009a, b), whereas the elliptic equation (34) is solved using a pseudo spectral method. The numerical grid is twodimensional and covers the volume of the star, i.e. the region and . The implementation uses a new radial coordinate , which is fitted to surfaces of constant chemical potential. This allows us to consider stars that are highly deformed by rotation. The perturbation variables are discretized on this grid and updated in time with a MacCormack algorithm. The numerical simulations are stabilised from high frequency noise with the implementation of a fourth order KreissOliger numerical dissipation. More technical details have been discussed in Passamonti et al. (2009a, b).
In order to solve the elliptic equation (34) with a spectral method, and save computational time, we set up a second numerical grid with lower resolution. This is important, since the spectral solver must be used at each time step, leading to a significant slowdown of the simulations. However, the lower resolution on the spectral grid does not affect the results, as spectral elliptic solvers provide highly accurate and rapidly convergent solutions already for relatively coarse grids (Grandclément & Novak, 2009). Therefore, at each time step we first fit the mass density perturbation on the spectral grid and then use the spectral routines to determine the gravitational potential perturbation . Subsequently, we fit the new value of to the original grid for the hyperbolic equations and carry on the evolution. The numerical code provides stable simulations for all rotating stellar models considered in this paper.
In this work, our choice of variables differs from that of Passamonti et al. (2009a). We evolve the velocity perturbations of the twofluids components instead of the “mass flux” perturbations of the comoving and countermoving degrees of freedom. The two formulations are obviously mathematically equivalent, but we wanted to develop a code based on the new set of variables in order to explore which formulation is best suited for future extensions. This is important, as we plan to add more realistic physics to our models by implementing an elastic crust region. As a first test, we compare the results of the new code to those obtained in Cowling approximation by Passamonti et al. (2009a). Neglecting the perturbation of the gravitational potential, i.e. setting , we find a complete agreement between the two numerical codes.
In order to study the spectral properties discussed below, in Sec. 6.1 and 6.2, we consider “generic” initial conditions that excite a large set of oscillation modes. For Type I perturbations we provide the following expression for the mass density:
(76) 
where neutrons and protons are initially countermoving. For Type II perturbations, we excite mainly normal and superfluid rmodes with the following initial data:
(77) 
where