Hydrodynamics of superfluid neutron stars

Hydrodynamics of rapidly rotating superfluid neutron stars with mutual friction


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 spin-up phase of pulsar glitches. We linearise the dynamical equations of a Newtonian two-fluid 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 non-stratified 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 non-axisymmetric 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 non-dissipative 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 pre-glitch state, and derive the equations that are used to extract the gravitational-wave signal.

methods: numerical – stars: neutron – stars: oscillation – star:rotation – gravitational waves

1 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 superfluids-superconductors (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 spin-up events seen in a number of young pulsars (and also some magnetars) during their magnetic slow-down 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 neutron-star seismology and a renewed assault on the problem of magnetic star oscillations, a seriously challenging problem from the theory point-of-view (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 high-density 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 counter-moving 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 multi-fluid dynamics.

The present work presents recent progress towards this goal. We study the oscillations of superfluid neutron stars by evolving in time the linearized two-fluid 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 gravitational-wave 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 counter-moving degrees of freedom are coupled in stratified models. From time-evolutions of the relevant perturbation equations, with the gravitational potential perturbation included, we determine the axi- and non-axisymmetric 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 non-axisymmetric 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 two-fluid 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 Euler-type equations and the Poisson equation for the gravitational potential (Prix, 2004). These take the form;


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 non-dissipative 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


where represents the bulk rotation (later we will assume that the two fluids co-rotate 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


that ensures Galilean invariance. The chemical potential and the entrainment parameter are then defined by


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:


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


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 counter-moving degrees of freedom are coupled in stratified neutron stars, which means that the gravitational-wave spectrum may contain the imprints of “superfluid” modes (see Sec. 6). This would not be the case in a non-stratified model.

The first EoS is determined by the following expression (Prix et al., 2002; Yoshida & Eriguchi, 2004; Passamonti et al., 2009a):


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 co-rotating 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):


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 non-rotating 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 non-rotating 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


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:


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 co-rotate 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 co-rotating 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 non-corotating models as small deviations from a co-rotating equilibrium. Subsequently, we use this approach to determine the initial conditions for hydrodynamical glitch evolutions.

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
Table 1: This table provides the main parameters for the two sequences of rotating models. The first column labels each model. In the second and third columns we give, respectively, the ratio of polar to equatorial axes and the angular velocity of the star. In the fourth column, the rotation rate is compared to the Kepler velocity that represents the mass shedding limit. The ratio between the rotational kinetic energy and gravitational potential energy and the stellar mass are given in the fifth and sixth columns, respectively. Finally, the seventh column gives the value of the chemical potential at the centre of the star. All quantities are given in dimensionless units, where is the gravitational constant, represents the central mass density and is the equatorial radius.

3.1 Corotating background

The equations that describe rapidly and uniformly rotating background models can be derived by imposing the conditions of stationarity and axi-symmetry on the Euler-type equations (2) and the Poisson equation (3(Prix et al., 2002; Yoshida & Eriguchi, 2004). This leads to


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


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 self-consistent 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 non-rotating model up to the mass shedding limit. In the numerical code, we re-write 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:


where the asterisk denotes the dimensionless quantities and . Combining equations (17) with the dimensionless mass , we can derive the equatorial radius of the star,


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:


where .

Figure 1: In this figure, we compare our numerical results to the analytical solution of Prix et al. (2002) for two slowly rotating stellar models with axis ratio (left panel) and (right panel). These two background stars are described by the EoS (10) and have the same proton fraction and symmetry energy term . The non-corotating corrections are determined by choosing the relative angular velocity  and imposing the constant central chemical potential condition (28). In the two panels, we show the radial profile of the perturbed neutron mass density  for three different angular directions, i.e. and . Our numerical results (solid line) agree very well with the analytical solution (empty circle) for the slowest rotating model (left panel). For faster rotating models the slow-rotation solution is expected to be less accurate. This is already evident for the case in the right panel, where the numerical and analytical solutions start to disagree.

3.2 Non-corotating solutions

In a multi-fluid 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 spin-up 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 spin-down 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 non-corotating quantities are considered as small deviations from a stationary, rapidly corotating configuration. We have already used these non-corotating corrections as initial data for studying the post-glitch 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


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 non-corotating corrections, Equations (14)–(15) become


The system of equations is closed by


that relates the mass density and the chemical potential perturbations for co-rotating backgrounds.

Non-corotating 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):


For solutions with constant mass, we impose a constraint on the mass of each fluid component, i.e.


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:


For the EoS (10), the adiabatic index is and the boundary condition (30) therefore reduces to:

Figure 2: This figure displays, for model C2, the non-corotating solutions of the proton and neutron mass density , in the left and right panel, respectively. The results correspond to constant mass solutions with parameters  and . These solutions were used as initial conditions by Sidery et al. (2010) for studying glitch hydrodynamics.

The system of equations (25)–(26) can be solved iteratively. First of all, for a given EoS we determine the co-rotating background with the self-consistent 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 non-corotating 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 slow-rotation 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 non-corotating 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 slow-rotation solution becomes less accurate as the star’s rotation increases. The same behaviour is found for non-corotating 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 non-corotating 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).

Figure 3: We show, for the C2 model and the non-corotating configuration from Fig. 2, the result for the dimensionless gravitational potential .

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


where is the azimuthal angle associated with the rotational motion, and the perturbed mutual friction force is (in the case of a co-rotating background)


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


where we recall that . By inverting these relations we can determine the velocity fields at any time step,


where .

The time evolution of the non-axisymmetric perturbation equations is a three-dimensional 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)


With this Fourier expansion the perturbation equations decouple with respect to and the problem becomes two-dimensional. In particular, for the axisymmetric case () only the component survives.

4.1 Boundary Conditions

In this work, we study axisymmetric () and non-axisymmetric 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 axi-symmetric and non-axisymmetric oscillations, we have to impose the following conditions, respectively :


For the velocity fields , we impose that there must be no mass flux across the origin () for both axisymmetric and non-axisymmetric perturbations:


At the rotational axis (), we impose the following conditions:


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:


Meanwhile, the Type II class is such that:


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.


The vector field is the Lagrangian displacement of the x-fluid 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 Gravitational-wave Extraction

In order to study the gravitational-wave signal emitted by pulsating superfluid neutron stars, we have implemented the quadrupole formula for both axisymmetric and non-axisymmetric oscillations. We will now discuss this implementation, in particular, the momentum and stress formula that we use to improve the numerical gravitational-wave extraction.

The gravitational-wave strain can be determined using the quadrupole formula (Thorne, 1980):


where is the pure spin tensor harmonic which has “electric-type” 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




The quantity is the quadrupole moment, in the case of a two-fluid star defined by;


where the spherical harmonics for the and cases are given by


where is the Legendre polynomial.

It is well-known 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 gravitational-wave 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:


where the quantity is defined by


We can reduce the order of the time derivative by using the method developed by Finn & Evans (1990), and obtain the perturbed momentum formula:


and the perturbed stress formula:


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).

Figure 4: We compare the gravitational-wave extraction results for axisymmetric and non-axisymmetric oscillations. The signal is generated by perturbing the stellar model C2, and the illustrated quantities are dimensionless. The left panel shows the dimensionless code quantity determined from three equivalent equations, respectively, the second time derivative of the quadrupole moment (dot-dashed line), the momentum-formula (solid-line) and the stress formula (dashed-line). In the right panel, we show the waveform of the non-axisymmetric oscillations for the perturbed C2 model. The upper and lower right panels displays respectively the real part and the imaginary part of the dimensionless quantity determined by the code. We compare the signal extraction to the momentum-formula (solid-line) and the stress formula (dashed-line).

For non-axisymmetric oscillations with , the two independent polarizations of the strain can be written as follows:


where is the spin-weighted spherical harmonics,


and we have defined the quantity


We can then re-write equation (62) as follows:




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:


For linear perturbations on a corotating background, we can further transform equation (65) into the following expression:


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):


By using Parseval’s Theorem we can write equation (67) for the and components as follows:


where , and is its Fourier transformation.

The characteristic strain of the gravitational-wave signal is then given by (Flanagan & Hughes, 1998):


where is the source distance. The strains and are related to the dimensionless quantities and used in the numerical code by the following expressions:


Similar relations provide the characteristic strain


As a first test of the numerical implementation, we compare the gravitational-wave extraction formulae for the axisymmetric and non-axisymmetric 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 gravitational-wave extraction routine. In the relativistic case, the linear perturbations of non-rotating 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 gravitational-wave 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 gravitational-wave strain is almost monochromatic and for a source at the maximal amplitude is .

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
Table 2: Comparison of the first three ordinary and superfluid mode frequencies and the Prix & Rieutord (2002) results. The star is the non-rotating C0 model, which corresponds to model III of Prix & Rieutord (2002), where the entrainment parameter is zero. Frequencies are given in units of and have been determined with an FFT of the gravitational strain time evolution. For this specific numerical simulation the frequency error bar is . The Prix and Rieutord values are denoted by PR. In the final column we show the relative error between our and PR results.

With our 2D Newtonian code, we then evolve in time non-radial oscillations of a superfluid non-rotating star for both the EoS (10) and (11). The kinetic energy of oscillating superfluid stars can be determined by the following expression:


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 gravitational-wave strain.

6 Results

Having formulated the time-evolution problem and described our implementation of the gravitational-wave 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 non-axisymmetric oscillations. We also provide a more detailed analysis of the gravitational-wave 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 two-dimensional 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 Mac-Cormack algorithm. The numerical simulations are stabilised from high frequency noise with the implementation of a fourth order Kreiss-Oliger numerical dissipation. More technical details have been discussed in Passamonti et al. (2009a, b).

Figure 5: This figure displays the effect of rotation on the quasi-radial and axisymmetric quadrupole modes. We use the sequence of models A, with entrainment parameter , proton fraction and vanishing symmetry energy term. On the horizontal axis the angular velocity is rescaled with the Kepler angular velocity , while the mode frequencies are given in dimensionless units and for a rotating frame. In the left panel, we show some “ordinary modes”, which are due to the co-moving degrees of freedom. We identify the and fundamental modes and the first three quasi-radial overtones and pressure modes. In the right panel, we show instead the “superfluid modes”, which correspond to the counter-moving degrees of freedom. In this case we show the modes up to the second overtones. For these non-stratified models, the ordinary and superfluid modes are decoupled.

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 slow-down 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 two-fluids components instead of the “mass flux” perturbations of the co-moving and counter-moving 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:


where neutrons and protons are initially counter-moving. For Type II perturbations, we excite mainly normal and superfluid r-modes with the following initial data: