The f-mode instability in neutron stars

# Nonlinear viscous damping and gravitational wave detectability of the f-mode instability in neutron stars

Andrea Passamonti Kostas Glampedakis

Theoretical Astrophysics, University of Tübingen, Auf der Morgenstelle 10, Tübingen D-72076, Germany
Departamento de Física, Universidad de Murcia, E-30100 Murcia, Spain
September 15, 2019
###### Abstract

We study the damping of the gravitational radiation-driven -mode instability in rotating neutron stars by nonlinear bulk viscosity in the so-called supra-thermal regime. In this regime the dissipative action of bulk viscosity is known to be enhanced as a result of nonlinear contributions with respect to the oscillation amplitude. Our analysis of the -mode instability is based on a time-domain code that evolves linear perturbations of rapidly rotating polytropic neutron star models. The extracted mode frequency and eigenfunctions are subsequently used in standard energy integrals for the gravitational wave growth and viscous damping. We find that nonlinear bulk viscosity has a moderate impact on the size of the -mode instability window, becoming an important factor and saturating the mode’s growth at a relatively large oscillation amplitude. We show similarly that nonlinear bulk viscosity leads to a rather high saturation amplitude even for the -mode instability. In addition, we show that the action of bulk viscosity can be significantly mitigated by the presence of superfluidity in neutron star matter. Apart from revising the -mode’s instability window we provide results on the mode’s gravitational wave detectability. Considering an -mode-unstable neutron star located in the Virgo cluster and assuming a mode amplitude at the level allowed by bulk viscosity, we find that the emitted gravitational wave signal could be detectable by advanced ground-based detectors such as Advanced LIGO/Virgo and the Einstein Telescope.

###### keywords:
methods: numerical – stars: neutron – stars: oscillation – star: rotation.
pagerange: Nonlinear viscous damping and gravitational wave detectability of the -mode instability in neutron starsNonlinear viscous damping and gravitational wave detectability of the -mode instability in neutron starspubyear:

## 1 Introduction

Neutron stars are exciting cosmic laboratories of matter at extreme conditions. However, unveiling the physics of these ultra-dense and relativistic systems is a real challenge even after more than four decades of astronomical observations. The situation is likely to improve significantly with the advent of gravitational wave astronomy. The forthcoming generation of improved gravitational-wave detectors, such as the Advanced LIGO and Virgo interferometers and the planned Einstein Telescope (ET), should be sensitive to various aspects of neutron star dynamics. The observation of neutron star oscillations (normal modes), usually dubbed neutron star “asteroseismology”, could become an excellent tool for probing the interior properties of these objects (Andersson & Kokkotas, 1998).

Of particular importance for neutron star asteroseismology are the oscillation modes that can become unstable (i.e. grow in amplitude) under the emission of gravitational radiation. This is the well-known CFS instability which develops when a retrogade mode (with respect to the stellar rotation) becomes prograde, as a result of rotational dragging, and then grows while emitting gravitational radiation (Chandrasekhar, 1970; Friedman & Schutz, 1978). Exhaustive work on the subject has shown that among the various stellar modes the ones that could become CFS-unstable in realistic neutrons stars, and likely to be strong sources of gravitational waves, are the inertial -mode and the fundamental -mode (see Andersson (2003) for a review).

The study of dissipative mechanisms that could limit or entirely suppress the gravitational wave instability of the - and -modes has become a small industry since the early days of this field (Andersson & Kokkotas, 2001). Given that astrophysical neutron stars are viscous systems, a clear understanding of the various dissipative mechanisms is a necessary ingredient of any realistic theoretical model. In fact, this is another facet of neutron star asteroseismology as the absence of unstable modes can be related to the strength of dissipation, thus providing indirect information about the stellar interior. For instance, this possibility has been recently illustrated by Ho, Andersson & Haskell (2011), who compared the latest theoretical parameter space of the -mode instability against the observed population of accreting neutron stars in low mass X-ray binaries.

The present paper makes another contribution to our understanding of the -mode instability and its damping. We focus on the action of bulk viscosity, which accounts for the energy drained from the oscillating fluid as a result of the induced departure from local beta equilibrium (Sawyer, 1989; Shapiro & Teukolsky, 1983). The novelty of our work is the study of bulk viscosity in the nonlinear regime (with respect to the oscillation amplitude), thus extending previous models based on amplitude-independent viscosity (e.g. Ipser & Lindblom, 1991; Lindblom, 1995; Gaertig et al., 2011).

The bulk viscosity formalism used in our analysis is based on the work by Alford and collaborators (Alford et al., 2010, 2011, 2012). These authors have recently studied the role of nonlinear bulk viscosity in the context of the -mode instability (we also briefly discuss the -mode below in Sec 4.3). A key result in these papers is an analytical expression for the bulk viscosity coefficient in the nonlinear regime which can be easily incorporated in the standard machinery for calculating viscous mode damping. This latter procedure consists of solving the inviscid hydrodynamics equations and determining the mode frequency and eigenfunctions, and subsequently using these results in the available expressions for the mode’s volume-integrated energy and damping rate (see section 2).

Our -mode calculation is numerical and is based on a 2D time-domain code developed by Passamonti et al. (2009a). The code evolves the linearised hydrodynamical equations of rapidly rotating stars in Newtonian theory. The frequency and eigenfunctions of specific modes can be then extracted from the time-evolved perturbations with a method developed by Stergioulas et al. (2004) and Dimmelmeier et al. (2006). An overview and a sample test of the code is provided in section 3. In this work we focus on the sextupole () -mode, which is known to be the dominant unstable multipole both in Newtonian and Relativistic stars (Friedman, 1983; Ipser & Lindblom, 1991; Gaertig et al., 2011). We calculate the resulting instability “window” (that is, the unstable portion of the parameter space described by the stellar rotational frequency and temperature) after accounting for the action of shear and bulk viscosity (section 4.1). Our calculation provides results for the -mode saturation amplitude as determined by the balance between nonlinear viscous damping and gravitational radiation-driven growth (section 4.1). We also calculate the gravitational wave strain associated with an unstable -mode and with its amplitude limited by nonlinear bulk viscosity (section 4.2). Our conclusions and a discussion of remaining issues can be found in section 5.

## 2 Formalism

The strategy for calculating the viscous damping and the gravitational wave-driven growth of oscillation modes is a well established two-step procedure (Ipser & Lindblom, 1991) and therefore it will only be outlined here. The first step consists of calculating the mode properties (frequency, eigenfunctions) of the inviscid system. In the second step, these results are inserted in the volume integral expressions for the damping/growth rates. This procedure relies on the implicit assumption that the mode evolves quasi-adiabatically and the dissipation and growth rates are much lower than the oscillation frequency.

The relevant linearised hydrodynamics equations for non-superfluid matter are the familiar Euler, mass conservation and Poisson equations. In a frame rotating with the stellar angular frequency these are:

 ∂tδvi+2ϵijkΩjδvk+∇i(δh+δΦ)=1ρDi, (1) ∂tδρ+∇j(ρδvj)=0, (2) ∇2δΦ=4πGδρ, (3)

where , and are, respectively, the perturbed velocity, mass density and gravitational potential. We have also introduced the perturbation of the specific enthalpy, , where is the fluid pressure perturbation. The term in the Euler equation (1) is the collective viscous force and will be discussed below.

The system of equations (1)-(3) is closed once an equation of state (EoS) is provided. In this work we consider a polytropic EoS

 P=Kργ, (4)

where is a constant and is related to the adiabatic index by .

Within the present single-fluid hydrodynamics model the viscosity force consists of shear and bulk viscosity contributions (Ipser & Lindblom, 1991),

 Di=∇j(νδσij)+∇i(ζδσ), (5)

where and are, respectively, the shear and bulk viscosity coefficient and

 δσij =12(∇iδvj+∇jδvi−23gij∇δσ), (6) δσ =∇jδvj, (7)

where is the metric tensor.

The viscosity force leads to the following dissipation rates (Ipser & Lindblom, 1991),

 ∂tEs=2∫dVηδσijδσ∗ij,∂tEb=∫dVζδσδσ∗, (8)

where the integrals are taken over the stellar volume and an asterisk denotes a complex conjugate. The gravitational radiation power is given by the standard multipole formula (Thorne, 1980),

 ∂tEgw=ω∑l≥2Nl(ω−mΩ)2l+1(|δDlm|2+|δJlm|2), (9)

where is the mode frequency as measured in the rotating frame. This expression features mass and current multipoles, denoted as and respectively, and the coupling constant

 Nl=4πGc2l+1(l+1)(l+2)l(l−1)[(2l+1)!!]2. (10)

The -mode oscillation radiates mainly through the mass multipole moments. These are given by

 δDlm=∫dVrlδρY∗lm, (11)

where is the standard spherical harmonic function with indices .

The above damping/growth rates can be readily calculated once the mode solution has been obtained from the inviscid Euler equation, i.e. equation (1) with . The inviscid equations also lead to the conserved mode energy

 E=12∫dV[ρδviδv∗i+12(δρδU∗+δρ∗δU)], (12)

where .

With viscosity included, the mode acquires a complex-valued frequency. For example, the perturbed velocity has a time profile

 δvi∼eiωt−t/τ, (13)

and similarly for the other parameters. As a result the mode energy evolves according to

 ∂tE=−2Eτ. (14)

The net combined effect of damping and growth can then be expressed as (Ipser & Lindblom, 1991)

 1τ=1τgw+1τs+1τb, (15)

where

 1τgw=∂tEgw2E,1τs=∂tEs2E,1τb=∂tEb2E, (16)

are the timescales associated with, respectively, gravitational wave growth, shear and bulk viscosity damping. A mode goes CFS-unstable when , which requires, as a necessary condition, a negative flux . This can happen only when , which is the mathematical statement of the mode changing from retrograde to prograde.

Up to this point our discussion has been essentially a repetition of well known theory. The following sections, however, offer a new analysis of the effect of bulk viscosity on the -mode instability. On the other hand there is nothing new here about shear viscosity, which anyway plays a minor role in our calculation. This is because shear viscosity has a significant impact on the instability in the low temperature regime () where the instability is anyway likely to be suppressed by superfluid vortex mutual friction (Lindblom & Mendell, 1995). Our modelling does not account for superfluidity. Without superfluidity the shear viscosity is dominated by neutron collisions and the resulting coefficient is (Flowers & Itoh, 1979; Ipser & Lindblom, 1991)

 η=347ρ9/4T−2g cm−1s−1. (17)

We use this coefficient for the shear viscosity in equations (8) and (16).

### 2.1 Nonlinear bulk viscosity

Considering a neutron star core composed of matter, the beta equilibrium is established through the modified Urca (mUrca) reactions

 ~N+n→~N+p+e−+¯νe,~N+p+e−→~N+n+νe, (18)

where is a spectator nucleon that guarantees overall energy-momentum conservation. In more massive neutron stars where the proton fraction is sufficiently large, , the more efficient direct Urca (dUrca) process (i.e. the reaction (18) without ) becomes energetically favourable (Lattimer et al., 1991). Although both processes are discussed in the following, most attention is given to the mUrca case which is the relevant one in a ‘canonical’ neutron star.

The coefficient of the bulk viscosity associated with beta reactions takes the form (Haensel et al., 2001)

 ζ0=C2ΓT2qω2, (19)

where is the oscillation frequency in the rotating frame, and represents the particle reaction rate (see Table 1). The temperature power-law is reaction-dependent and takes the value () for mUrca (dUrca). The quantity is the so-called strong interaction susceptibility

 C≡ρ∂μΔ∂ρ∣∣∣xp, (20)

where is the difference between the particle chemical potentials on the two sides of (18). For the simple case of neutron matter modelled as a gas of free (non-interacting) hadrons we have in natural units (Alford et al., 2010)

 C=(3π2ρ)2/36m5/3n, (21)

where is the neutron mass. The resulting viscosity coefficient for this type of matter has been calculated by Sawyer (1989) and is the most widely used coefficient in the literature:

 ζ0=6×10−59ρ2ω−2T6g cm−1s−1. (22)

Strictly speaking, the bulk viscosity formulae (19) and (22) are valid for a ‘small’ oscillation amplitude in the so-called sub-thermal regime, where . However, if the oscillation amplitude is sufficiently large the system may enter the supra-thermal regime . As discussed in Alford et al. (2010, 2011), this condition may be actually reached at relatively small amplitudes . The Lagrangian density perturbation is defined in terms of the displacement vector as .

In this nonlinear bulk viscosity regime, the formula (19) is replaced by the amplitude-dependent expression (Alford et al., 2010, 2011),

 ζ=ζ0[1+q∑j=1(2j+1)!!χj2j(j+1)!(CTΔρρ)2j], (23)

where are matter/reaction-dependent parameters given in Table 1. This expression is approximate, but remains accurate as long as the amplitude is not too large. For mUrca (dUrca) bulk viscosity the relevant amplitude threshold is ((Alford et al., 2010, 2011). For an oscillation amplitude above these thresholds equation (23) becomes inaccurate, failing to capture the nonlinear saturation of the supra-thermal bulk viscosity and the subsequent decreasing behaviour (see figure 1 in Alford et al., 2011). For the purpose of this paper, where the fluid perturbation remains in the linear regime (i.e. ), the approximation (23) remains accurate and, in fact, includes the maximum viscosity attained in the supra-thermal regime.

The nonlinear coefficient (23) provides the main microphysical input for our -mode instability calculation. A large amplitude -mode is likely to ‘activate’ the nonlinear terms in (23) and suffer an enhanced viscous dissipation. If the resulting damping rate is sufficiently high then the mode may actually saturate and stop growing at some maximum amplitude. In the following sections we explore to what extent this scenario could be relevant for unstable -modes of rapidly rotating neutron stars.

In the following, is calculated assuming the simple case of a free hadron gas, see Table 1 for the relevant parameters. For the proton fraction (which is assumed uniform throughout the core) we take a typical value when mUrca is active and a higher value in the case of dUrca viscosity.

Given the functional form of it is also necessary to define an appropriate mode-amplitude. We introduce the amplitude parameter defined by the following expression:

 δvreq=αωReqYlm(π/2,0), (24)

where is the equatorial radius of the star, and is the radial component of the velocity perturbation determined at and . Note that our definition of the mode amplitude is similar to the -mode definition used by Owen et al. (1998), provided that is replaced by .

## 3 Neutron star model

Our analysis is based on a Newtonian polytropic stellar model with uniform rotation. More specifically, we study two polytropes with indices and , which provide a reasonable approximation to realistic equations of state. The axisymmetric equilibrium configurations are determined by solving the stationary equations with the self-consistent method developed by Hachisu (1986). A detailed discussion of the numerical techniques can be found in Passamonti et al. (2009a, b). In the following calculations we fix the stellar mass at and the constant (in cgs units) at () for the () model. The resulting stellar radius of the nonrotating model is then  km ( km) for (). As discussed by Ipser & Lindblom (1991) and Lindblom (1995), the -mode damping/growth timescale in neutron stars with different mass and radius can be obtained by a simple rescaling.

Although our -mode calculation assumes a fluid star we implicitly account for the presence of a solid crust by limiting the action of bulk viscosity in the liquid core, where the reactions (18) take place. This means that the volume integral in equation (8) is taken only over the core. This truncation should be accurate below the crust melting temperature  (Shapiro & Teukolsky, 1983), a region which overlaps with the temperature range where the -mode instability is likely to operate. In our model the “crust” occupies the outer layer of the star with its boundary located at the transition density .

A step of the calculation that deserves some special discussion concerns the constraints imposed on the -mode amplitude by the chosen stellar models. A first requirement is that the mode velocity in the inertial frame should not exceed the speed of light. In the rotation range relevant for the instability of the -mode (which is the most unstable multipole), we find that this condition is satisfied when () for the () model.

Another physically motivated (albeit less rigorous) condition is imposed by the linear perturbation formalism used here. We can quantify this requirement in terms of the Lagrangian density perturbation , as our perturbation formalism (section 2) should break down when . This condition also limits the use of the approximate bulk viscosity coefficient given in equation (23). Indeed, as detailed in Alford et al. (2011), this approximation for the mUrca (dUrca) bulk viscosity is accurate only when (). In other words, we need to make sure that remains within the limits of linearised theory.

The quantity is shown in figure 1 for the -mode with a fiducial amplitude . The density eigenfunction is calculated at the equatorial plane (), and in two radial locations, namely, the surface and the crust-core interface. All the rapidly rotating models shown in figure 1 have at the surface. On the other hand, the perturbation remains ‘well-behaved’ at the crust-core boundary, i.e. . The same is true for every point inside the stellar core. Therefore our analysis can be safely used for studying the bulk viscosity in the core.

Another way of gauging the size of the mode amplitude is by comparing the mode energy defined in equation (12) against the bulk rotational energy of the star . It is reasonable to expect that in the instability regime we have , as the latter energy is the one tapped by the instability. In figure 1 we show these two quantities for an -mode with . Near the maximum rotation rate, the mode energy is much smaller than the rotational energy, , but becomes a significant fraction of it at relatively lower rotations. For instance, the condition is violated in the model when . However this matters little for our calculation since, as we will see in section 4.1, these stars lie outside the instability window of the -mode. In summary, our -mode calculations preserve self-consistency with respect to the mode amplitude (24) and the approximate supra-thermal bulk viscosity coefficient (23).

### 3.1 Numerical Code

The evolution of the perturbation equations (1)-(3) is a three-dimensional problem in space. In order to have an easier numerical implementation, we can exploit the background symmetry and expand the non-axisymmetric perturbations in a Fourier series with respect to the azimuthal harmonic index . As a result, the problem becomes two-dimensional (Papaloizou & Pringle, 1980; Passamonti et al., 2009a, b). For the mass density perturbation the expansion is given by

 δρ(t,r,θ,ϕ)=m=∞∑m=0[δρ+m(t,r,θ)cosmϕ+δρ−m(t,r,θ)sinmϕ], (25)

and similarly for the other variables. For any , we evolve a system of ten partial differential equations for the ten variables , while the enthalpy perturbation, , is determined by an EoS. For a polytropic star it is given by

 δh=c2sδρρ, (26)

where is the speed of sound,

 c2s=∂P∂ρ=Kγργ−1. (27)

The linearised equations are evolved numerically with a code developed by Jones et al. (2002) and Passamonti et al. (2009a, b), where the reader can find all the details. From the time evolution of each perturbation variable, we determine the mode frequency with a Fast Fourier Transformation (FFT) and then extract the eigenfunctions with the method developed by Stergioulas et al. (2004) and Dimmelmeier et al. (2006). Both the frequencies and eigenfunctions are used in section 4.1 to determine the damping time of the bulk and shear viscosity and the growth time of the gravitational-wave-driven instability.

### 3.2 Boundary Conditions

The time evolution of non-axisymmetric oscillations requires the prescription of boundary conditions. We study rotating stars with equatorial and rotation axis symmetry, where the two-dimensional numerical grid extends over the region and . The quantity is the star’s radius at a given latitude . At each time step, we must then specify the variables at the surface, origin, rotational axis and equator.

At the origin () and the rotational axis (), the perturbation variables and equations must be regular. For the non-axisymmetric modes the regularity condition leads to the following equations:

 δρ=δh=δΦ=0,andδvi=0, (28)

both at and .

At the equator (), the reflection symmetry divides the perturbations into two sets with opposite parity (Passamonti et al., 2009a). For the -mode the parity conditions read

 ∂δρ∂θ=∂δh∂θ=∂δΦ∂θ=0, (29) ∂δvr∂θ=δvθ=∂δvϕ∂θ=0. (30)

We impose at the surface, , the vanishing of Lagrangian enthalpy perturbation,

 Δh=δh+ξi∇ih=0, (31)

which is the standard free surface condition used in stellar oscillations. The vector field is the Lagrangian displacement, and the value of the perturbed enthalpy at the surface is determined at each time step from equation (31).

### 3.3 Numerical tests

We test our numerical framework using the results of Ipser & Lindblom (1990, 1991) and Lindblom (1995). These papers study the -mode instability for various polytropic Newtonian models, and use an eigenvalue method to extract the frequencies and eigenfunctions of the -mode.

First of all, we determine the frequency of the onset of the CFS instability for inviscid stars. This corresponds to models with vanishing -mode frequency in the inertial frame, i.e. . We find for and for . These values agree to better than 0.5% with the results reported in Ipser & Lindblom (1990). Note that is the frequency unit used by Ipser & Lindblom (1990), where denotes the average mass density of the non-rotating model. In figure 2 (left panel) we show the -mode frequencies determined with our numerical code.

During the tests we have noticed a slight difference between the maximum rotation rate determined with our numerical code and that reported in Ipser & Lindblom (1991). Although the difference is very small (less than 0.6%) it has a visible impact on the instability window. For the and polytropes, we find and , respectively, while Ipser & Lindblom (1991) report and . We have no reason to question the reliability of our code, given that it reproduces to high accuracy the axisymmetric background configurations of Hachisu (1986). Hence, we think that this small discrepancy is due to the different numerical techniques used in Ipser & Lindblom (1991).

As a final test, we compare the damping/growth times of the gravitational radiation and the viscous damping times with the data extracted from Ipser & Lindblom (1991) and Lindblom (1995). The results are shown in figure 2 for the rotating models with temperature  K. We have rescaled the data reported in Ipser & Lindblom (1991) with the of our equilibrium configuration. By using the standard bulk viscosity given in equation (22), we find a very good agreement with previous results.

## 4 Results

We are now ready to present our results on the impact of nonlinear bulk viscosity on the gravitational wave-driven -mode instability. We first discuss how the mode’s instability window is modified as a result of viscous damping (section 4.1). We then calculate the gravitational wave strain associated with the -mode and provide upper limits for its detectability by present day and future detectors (section 4.2). In the final part of our calculation (section 4.3) we make a brief digression and discuss the role of nonlinear bulk viscosity in the context of the -mode instability.

### 4.1 Instability window

The boundary of the -mode instability window on the plane is determined by

 1τ=0. (32)

For a given this condition provides a critical rotational frequency above which the mode is CFS-unstable.

The previous studies of the -mode instability included bulk viscosity in the sub-thermal regime, typically using Sawyer’s coefficient (22(Sawyer, 1989). In that case equation (32), and consequently , are independent of the mode amplitude . However, once the full nonlinear bulk viscosity expression (23) is used, the instability window depends explicitly on .

We first consider bulk viscosity due to the mUrca process. In figure 3 we show the -mode instability window for our and polytropic models and for different mode amplitudes. The effect of nonlinear bulk viscosity is clearly more pronounced in the model, and for amplitudes . For an amplitude as large as the window area is decreased to a significant degree, resulting in a minimum . The stiffer neutron star model is less sensitive to nonlinear viscous damping. This is not surprising given that the instability’s growth time is shorter than the one in the polytrope. To summarise, we can conclude that nonlinear bulk viscosity due to mUrca reactions has a moderate effect on the -mode-instability even for the largest mode amplitudes consistent within our framework.

The instability window in models with dUrca reactions is shown in figure 4. As expected, the higher efficiency of the dUrca reactions leads to a stronger bulk viscosity which suppresses the instability above . This temperature is very close to the expected onset of neutron superfluidity in the core, and once neutron superfluidity is present the instability is suppressed as a result of vortex mutual friction (Lindblom & Mendell, 1995). The resulting -mode instability window is much smaller than the (already small) window of the mUrca models. Focusing on the damping of nonlinear bulk-viscosity we find that the situation is similar to the previous mUrca models. The critical frequency is moderately increased, more for the polytrope and less for the one.

Our results suggest that bulk viscosity, despite its increased efficiency due to the nonlinear supra-thermal contribution, can only cause a moderate degree of additional damping. How robust is this conclusion? Based on recent results on the -mode instability in relativistic stars (Gaertig et al., 2011) which generally predict a shorter growth timescale than the Newtonian one (), one would suspect that the role of nonlinear bulk viscosity is even less important. At the same time, however, we also need to account for the fact that bulk viscosity itself becomes stronger, roughly by a factor , as we move to a more realistic model with interacting hadrons (Alford et al., 2010, 2012). With these two effects combined and the ensuing balance between them (at least approximately) it is likely that our Newtonian analysis leads, after all, to robust results. Obviously, a definitive answer can only be reached by means of a more detailed calculation (e.g relativistic stellar models with a more realistic equation of state).

As we have already pointed out, our analysis neglects the presence of neutron and proton superfluidity. At the same time, however, we know that superfluidity plays a key role by suppressing the -mode instability in sufficiently cold neutron stars through vortex mutual friction (Lindblom & Mendell, 1995). Bulk viscosity is also directly affected by superfluidity as a consequence of the reduced rates of the mUrca/dUrca reactions in superfluid neutron stars. Thus we have good reasons for revising the bulk viscosity damping in superfluid neutron stars. Within our framework we can do that by introducing the appropriate superfluid correction factor in the viscosity coefficient (see below). Other important aspects of superfluid hydrodynamics, like damping due to vortex mutual friction and other viscous degrees of freedom (Andersson & Comer, 2006), can not be accessed without the full machinery of multi-fluid hydrodynamics (e.g. Andersson, Glampedakis & Haskell, 2009).

As a newborn neutron star cools down the transition to proton superconductivity is likely to occur first, with neutrons in the core becoming superfluid at a later stage. This order is the one suggested by pairing calculations (e.g. Kaminker et al., 2002). The recent models for the observed thermal evolution of the Cassiopeia A neutron star provide constraints for the critical temperatures for the onset of proton and neutron superfluidity in the core, and  (Page et al., 2011; Shternin et al., 2011). There is therefore a temperature range where the stellar core contains normal neutrons and superconducting protons. In this regime, where vortex mutual friction cannot yet operate, it makes sense to study how bulk viscosity is modified by proton superconductivity.

This can be achieved by using the modified viscosity coefficient

 ζsup=ζR, (33)

where is the bulk viscosity coefficient for normal matter defined by equation (23). The factor represents the superfluidity-induced reaction rate suppression and is a function of the critical temperature for the onset of proton superconductivity. The explicit form of is given by equation (34) in Haensel et al. (2001). It should be pointed out that the modification (33) is the appropriate one for the bulk viscosity coefficient in the sub-thermal regime. As far as we know there is no similar result available for the supra-thermal , hence the superfluid modification (33) should be viewed as an approximation.

The impact of proton superconductivity on the -mode window is shown in figure 5. We see that proton superconductivity can mitigate to a significant degree the damping action of nonlinear bulk viscosity. This effect becomes more prominent with an increasing critical temperature .

Nonlinear bulk viscosity, when efficient, is a mechanism for saturating a growing mode due to its intrinsic dependence on the mode-amplitude. Of course this does not exclude the simultaneous presence of other competing (and perhaps more efficient) saturation mechanisms. A strong candidate alternative saturation mechanism is nonlinear mode coupling. For instance, mode couplings with other inertial modes is known to saturate the -mode instability at small amplitudes, irrespectively of rotation/temperature (Arras et al., 2003; Bondarescu et al., 2007). Thus it is also possible that a similar effect could be at work in the case of the -mode. Up to now there is limited work on the nonlinear hydrodynamics of the -mode (Ou et al., 2004; Shibata & Karino, 2004; Kastaun et al., 2010), and none for the multipole which is of interest to us.

Still, we can obtain some hints from the available work on the quadrupole -mode. The nonlinear numerical simulations of Kastaun et al. (2010) suggest that the -mode amplitude in relativistic polytropic stars is generically limited by wave breaking at the surface and mode couplings. Considering a rotating polytropic model with axis ratio , where is the polar radius, Kastaun et al. (2010) found that an -mode with an initial energy is efficiently damped by nonlinear effects on a timescale of about 10 ms. Assuming that this result remains valid also for the -mode we can deduce (with the help of figure 1 and noting that the above axis ratio corresponds to ) that nonlinear mode couplings could saturate the mode amplitude at . This exercise suggests that nonlinear hydrodynamics may provide the dominant saturation mechanism for the -mode instability. However, without a detailed analysis of the -mode dynamics we cannot yet draw any secure conclusion.

### 4.2 Gravitational wave signal

The final part of our -mode instability analysis is devoted to the calculation of the gravitational wave strain and its detectability. Unlike previous studies (Lai & Shapiro, 1995; Ou et al., 2004; Shibata & Karino, 2004) we concentrate on the -mode which is, as we have already pointed out, the most instability-prone multipole in both Newtonian and relativistic stars (Gaertig et al., 2011).

The gravitational wave strain can be computed by means of the standard multipole formula (Thorne, 1980)

 h44ij=Gc61rd4Idt444T44ij, (34)

where is proportional to the mass moment (11),

 I44=8π945√5δD44, (35)

and is the corresponding spin tensor harmonic angular function of the “electric-type” (Thorne, 1980).

It turns out that the two independent wave polarizations can be expressed in terms of the spin-weighted spherical harmonic as

 h44θθ−ih44θϕ=h44−2Y44, (36)

where the wave amplitude is defined as

 h44≡√109454πrGc6d4dt4δD44. (37)

For a monochromatic wave equation (37) reduces to

 h44=√109454πrGc6ω4I∫dVδρ(r)r4Y∗44eiωIt=^heiωIt, (38)

where is the mode frequency in the inertial frame.

As a representative gravitational wave amplitude we numerically calculate using the relevant -mode frequency and eigenfunctions. In figure 6 (left panel) we show our results for assuming a source located at a distance of 20 Mpc and a mode amplitude . As expected, the gravitational-wave signal is stronger for models with , where the star is more unstable and is higher.

For assessing the actual detectability of a neutron star undergoing oscillations due to an unstable -mode we need to take into account the cumulative effect of the observed number of oscillation cycles in a given frequency bandwidth. This information is encoded in the characteristic wave strain

 hc≡⟨∣∣h44−2Y44∣∣⟩√ν2∣∣∣dtdν∣∣∣ (39)

where and denotes the average over the angles .

The time profile of the mode frequency can be determined following the two-stage evolution of the instability (Owen et al., 1998). The onset of the instability leads to the first (linear) phase of exponential growth where . Eventually this growth arrives to a halt as the mode is saturated by nonlinear physics (for instance, the nonlinear bulk viscosity studied in this paper). During this second (nonlinear) phase the mode amplitude remains constant and the emitted gravitational radiation removes angular momentum from the star’s bulk rotation causing its spin-down. This evolution leads to a monotonically decreasing mode frequency . Here we consider only the second nonlinear saturation phase and neglect the thermal evolution of the star and other mechanisms which can accelerate the spin-down. Our results provide therefore an approximate upper limit for the gravitational wave strain. A more detailed analysis of the detectability of unstable -modes in neutron stars will be presented in future work.

In order to estimate the number of oscillations near the given mode frequency we need to consider the angular momentum balance for the system,

 dJdt=−dJdt∣∣∣gw, (40)

where is the stellar angular momentum 111To be more precise we should also consider the variation of the mode’s canonical angular momentum, but its contribution is generally smaller than the stellar angular momentum (Owen et al., 1998) and hence can be neglected for simplicity.. Writing , and solving with respect to the derivative of the mode frequency, we obtain

 dνdt=−12πdωdΩdJdt∣∣∣gw(dJdΩ)−1. (41)

The radiated angular momentum is given by (Thorne, 1980)

 dJidt∣∣∣gw=2N4ω9I|δD44|2δiz. (42)

where is the Kronecker delta. From this expression it is obvious the scaling , which means that . Then from (39) we can deduce that the characteristic strain is independent of the mode amplitude .

The above calculation is idealised in the sense that it implicitly assumes a gravitational wave observation lasting as long as the -mode instability evolution timescale. This timescale is given by (Owen et al., 1998)

 τevol=∣∣∣νdtdν∣∣∣. (43)

For a fiducial mode amplitude and a model, we find that for the evolution timescale is  weeks, and increases to  yr for slower rotating models with . The evolution timescale becomes even longer for modes with given the scaling . Thus in the most relevant cases the mode evolution could be much longer than the typical observation time of LIGO and Virgo. What that means is that from the detectors’ point of view the signal can be considered as essentially monochromatic. Besides the temporal limits set by the detectors’ capability, there might be additional ones from the physical system itself. For instance, the star may efficiently spin-down as a result of a strong magnetic field, thus quenching the -mode instability and the emission of gravitational radiation.

Therefore, in order to provide a ‘realistic’ estimate for the gravitational wave strain associated with the -mode instability we assume a monochromatic signal of duration , for a source at a 20 Mpc distance. For the purpose of comparison we also consider the ideal case of a signal lasting a time period . The results are shown in figure 6 together with the sensitivity curves of the current and future generation Earth-based gravitational wave detectors. The detectors’ noise curves are computed in the standard way using , where is the detector’s power spectral density. In our simplified spin-down model we have neglected the effects of the neutron star’s cooling. The curves shown in figure 6 (right panel) may represent therefore the signal emitted by an unstable star which moves at constant temperature from the Kepler frequency down to the minimum of the instability window (see figure 3). However the thermal evolution may significantly affect the trajectory through the instability window and the star may therefore leave the instability region at larger rotation rates and thus at higher mode-frequencies than those indicated in figure 6.

Our results suggest that the -mode instability may be detectable from the future ET for a mode amplitude . Detection by Advanced LIGO/Virgo requires a larger amplitude and stars with a stiff equation of state (in the figure the polytrope is clearly favoured over the one).

### 4.3 Digression: the r-mode instability

So far we have discussed the role of (nonlinear) bulk viscosity in the damping of the -mode instability. However, as mentioned in the introduction, the physics of nonlinear bulk viscosity has already been discussed in the context of the -mode instability (Alford et al., 2010, 2011, 2012). The main conclusion of that recent work is that bulk viscosity operating in the supra-thermal regime (see section 2.1) could suppress the -mode instability (regardless of the stellar temperature), and prevent the growth of the mode amplitude above some maximum saturation level. In this section we revisit the topic of -mode nonlinear viscous damping with the aim of providing a better understanding of the recent results. Why this is necessary will become obvious in the following discussion.

To begin with, we introduce the ‘standard’ dimensionless -mode amplitude defined as (Owen et al., 1998)

 δvj=αΩR(rR)lYBjlmeiωt, (44)

where is the magnetic vector harmonic. Obviously must not be confused with the -mode amplitude defined in equation (24). This expression describes the -mode velocity at leading order with respect to the stellar spin in the familiar slow-rotation approximation (see Andersson & Kokkotas, 2001, for further discussion). In this approximation the star maintains its spherical shape and the mode’s energy is dominated by the kinetic part, i.e.

 E=12α2Ω2R2−2m∫R0ρr2m+2dr, (45)

where the -mode is considered. The calculation of the bulk viscosity damping timescale requires the input of the density fluctuation associated with the perturbation (44). After adding a missing factor in Owen et al. (1998) (see also a comment in Lindblom et al., 1999)

 δρρ ≈α(ΩR)2c2s2m+1m(m+1)√2m+3 ×[2m2m+1√mm+1(rR)m+1+δΦ0]Ymm+1eiωt, (46)

where is the radial eigenfunction of the perturbed gravitational potential. At this point it is important to make a clarification. In their study of the -mode viscous damping Alford et al. (2011, 2012) have used a different definition for the -mode amplitude. More specifically, Alford et al. use the amplitude defined in equation (4.1) of Lindblom et al. (1999).222Shortly after the submission of this paper, Alford et al. have revised their r-mode calculation and used the more suitable amplitude defined by equation (44). The relation between the two different definitions of the -mode amplitude is (see also the footnote in Lindblom et al., 1999)

 α=~α√πm(m+1)3(2m+1)!. (47)

For the most unstable -mode we find .

The fact that Alford et al. (2011, 2012) use the notation for the non-standard -mode amplitude may be the cause for some confusion. According to their results, e.g. fig. 2 in Alford et al. (2011), bulk viscosity in the supra-thermal regime is capable of suppressing the -mode instability for a wide range of stellar temperatures provided . However, in terms of the standard -mode amplitude (44) this result corresponds to a saturation amplitude . Obviously this is a rather large amplitude333This can be demonstrated if we take a polytropic star with mass and radius  km, rotating at the Kepler limit (the maximal rotation rate can be approximated with , where is the average density at ). We then obtain and . Hence, the -mode saturation amplitude obtained in (Alford et al., 2011, 2012) would imply a mode energy comparable to, or even larger, than the stellar rotational energy., much higher than the saturation amplitude due to the -mode’s nonlinear coupling with other inertial modes (Arras et al., 2003; Bondarescu et al., 2007). Thus, we are led to conclude that for an amplitude the -mode instability is affected very little by the nonlinear contribution of bulk viscosity.

## 5 Concluding remarks

We have studied the gravitational wave-driven instability of the -mode in Newtonian polytropic neutron star models. The novelty of our work is the inclusion of bulk viscosity in the supra-thermal regime. We have found that the ensuing nonlinear enhancement of bulk viscosity leads to a moderate change in the -mode instability window. For the particular case of the -mode (the most unstable multipole) we have shown that nonlinear bulk viscosity associated with the modified Urca process has a notable effect provided the mode amplitude, as defined in equation (24), is (see figure 3). This amplitude translates to a mode energy for , increasing to for . However, even for the extreme case of a mode amplitude as large as bulk viscosity cannot entirely suppress the instability.

If the more powerful direct Urca reactions can operate (which would be the case in the most massive neutron stars) then the resulting -mode instability window is significantly reduced already in the sub-thermal bulk viscosity regime. The transition to the supra-thermal regime has again a moderate effect to the shape of the instability window. It turns out that in the temperature range where the change is important the instability is already likely to be suppressed by superfluid mutual friction (see figure 4).

The presence of superfluid components in neutron star matter could diminish the impact of bulk viscosity by suppressing the rates of the chemical reactions involved. By considering a temperature regime where the protons are superfluid and the neutrons normal and by treating in an approximate manner the superfluid suppression, we have shown that nonlinear bulk viscosity is significantly weakened (see figure 5).

Apart from calculating the bulk viscosity-modified -mode instability window we have also provided upper limits for the gravitational wave signal from a mode with , i.e. the amplitude range where nonlinear bulk viscosity can partially saturate the mode (figure 6). Our results suggest that the -mode instability (operating in a neutron star at the fiducial distance of 20 Mpc) could be detectable by a third generation instrument like the ET for an amplitude and an observation time of week. Detectability by Advanced LIGO/Virgo requires a rather high amplitude () and/or a longer observation time (e.g. the evolution timescale (43)). Moreover, the -mode detectability is moderately sensitive to the stellar polytropic index, with lower indices (stiffer equation of state) being favoured. It should be pointed out that our work provides only ‘static’ upper limits, in the sense that the combined rotational/temperature evolution of an -mode unstable neutron star has not been properly incorporated. A more detailed analysis of this issue will the subject of future work.

Our brief excursion to the study of the damping of the -mode instability by nonlinear bulk viscosity has revealed that a large mode amplitude (, see equation (47)) is required for the effect to be relevant. Thus, most likely, the saturation of the -mode is controlled by nonlinear mode couplings rather than bulk viscosity. A similar statement could also be true for the -mode but this remains to be (dis)proved, given the lack of any concrete results on the -mode nonlinear couplings. Some light on this issue can be shed by a recent study of the nonlinear dynamics of the -mode (Kastaun et al., 2010). Extrapolation of those results to the mode suggests that nonlinear hydrodynamics could outperform bulk viscosity in saturating a growing mode.

Besides the obvious need for a better understanding of the -mode nonlinear dynamics, several other improvements are required in this research area. For instance, our own analysis can be clearly improved by moving from Newtonian to General Relativistic gravity and from polytropic matter to more realistic equations of state. Indeed, as suggested by recent work on the relativistic -mode, the instability’s growth time may shrink as much as a factor compared to the Newtonian result. However, as we pointed out in section 4.2, this effect is likely to be counter-balanced by a similar increase in the bulk viscosity strength as a result of using more realistic equations of state. The presence of a solid crust adds one more level of complexity to the problem. Despite its central role in limiting the -mode instability (Andersson & Kokkotas, 2001), the dynamics of the crust has been totally ignored in the context of the -mode instability. Another interesting possibility arises when exotic matter (for example, deconfined quarks) is present in the stellar interior. Previous work on the -mode instability in neutron stars with exotic cores has unveiled several important differences with respect to normal hadronic neutron stars (e.g. Madsen, 1992; Andersson et al., 2002). This is likely to be the case also for the -mode instability. Clearly, all the above issues are of great interest and we plan to address some of them in the near future.

## Acknowledgements

We thank L. Lindblom for providing us with the numerical data shown in figure 2. We also would like to thank N. Andersson and K. Kokkotas for useful feedback on our paper, and W. Kastaun for fruitful discussions. AP acknowledges support from the German Science Foundation (DFG) via SFB/TR7. KG is supported by the Ramón y Cajal Research Programme of the Spanish Ministerio de Ciencia e Innovación

## References

• Alford et al. (2010) Alford M. G., Mahmoodifar S., Schwenzer K., 2010, Journal of Physics G Nuclear Physics, 37, 125202
• Alford et al. (2011) Alford M., Mahmoodifar S., Schwenzer K., 2011, in F. J. Llanes-Estrada & J. R. Peláez ed., American Institute of Physics Conference Series Vol. 1343 of American Institute of Physics Conference Series, Non-linear viscous saturation of r-modes. pp 580–582
• Alford et al. (2012) Alford M. G., Mahmoodifar S., Schwenzer K., 2012, Phys. Rev. D, 85, 044051
• Andersson & Comer (2006) Andersson N., Comer G. L., 2006, Class. Quantum Gravity, 23, 5505
• Andersson et al. (2009) Andersson N., Glampedakis K., Haskell B., 2009, Phys. Rev. D, 79, 103009
• Andersson et al. (2002) Andersson N., Jones D. I., Kokkotas K. D., 2002, MNRAS, 337, 1224
• Andersson & Kokkotas (1998) Andersson N., Kokkotas K. D., 1998, MNRAS, 299, 1059
• Andersson & Kokkotas (2001) Andersson N., Kokkotas K. D., 2001, International Journal of Modern Physics D, 10, 381
• Arras et al. (2003) Arras P., Flanagan É. É., Morsink S. M., Schenk A. K., Teukolsky S. A., Wasserman I., 2003, ApJ, 591, 1129
• Bondarescu et al. (2007) Bondarescu R., Teukolsky S. A., Wasserman I., 2007, Phys. Rev. D, 76, 064019
• Chandrasekhar (1970) Chandrasekhar S., 1970, Phys. Rev. Lett., 24, 611
• Dimmelmeier et al. (2006) Dimmelmeier H., Stergioulas N., Font J. A., 2006, MNRAS, 368, 1609
• Flowers & Itoh (1979) Flowers E., Itoh N., 1979, ApJ, 230, 847
• Friedman (1983) Friedman J. L., 1983, Phys. Rev. Lett., 51, 11
• Friedman & Schutz (1978) Friedman J. L., Schutz B. F., 1978, ApJ, 222, 281
• Gaertig et al. (2011) Gaertig E., Glampedakis K., Kokkotas K. D., Zink B., 2011, Phys. Rev. Lett., 107, 101102
• Hachisu (1986) Hachisu I., 1986, ApJSS, 61, 479
• Haensel et al. (2001) Haensel P., Levenfish K. P., Yakovlev D. G., 2001, A&A, 372, 130
• Ho et al. (2011) Ho W. C. G., Andersson N., Haskell B., 2011, Phys. Rev. Lett., 107, 101101
• Ipser & Lindblom (1990) Ipser J. R., Lindblom L., 1990, ApJ, 355, 226
• Ipser & Lindblom (1991) Ipser J. R., Lindblom L., 1991, ApJ, 373, 213
• Kaminker et al. (2002) Kaminker A. D., Yakovlev D. G., Gnedin O. Y., 2002, A&A, 383, 1076
• Jones et al. (2002) Jones D. I., Andersson N., Stergioulas N., 2002, MNRAS, 334, 933
• Kastaun et al. (2010) Kastaun W., Willburger B., Kokkotas K. D., 2010, Phys. Rev. D, 82, 104036
• Lai & Shapiro (1995) Lai D., Shapiro S. L., 1995, ApJ, 442, 259
• Lattimer et al. (1991) Lattimer J. M., Prakash M., Pethick C. J., Haensel P., 1991, Phys. Rev. Lett., 66, 2701
• Lindblom (1995) Lindblom L., 1995, ApJ, 438, 265
• Lindblom & Mendell (1995) Lindblom L., Mendell G., 1995, ApJ, 444, 804
• Lindblom et al. (1999) Lindblom L., Mendell G., Owen B. J., 1999, Phys. Rev. D, 60, 064006
• Ou et al. (2004) Ou S., Tohline J. E., Lindblom L., 2004, ApJ, 617, 490
• Owen et al. (1998) Owen B. J., Lindblom L., Cutler C., Schutz B. F., Vecchio A., Andersson N., 1998, Phys. Rev. D, 58, 084020
• Page et al. (2011) Page D., Prakash M., Lattimer J. M., Steiner A. W., 2011, Phys. Rev. Lett., 106, 081101
• Papaloizou & Pringle (1980) Papaloizou J. C., Pringle J. E., 1980, MNRAS, 190, 43
• Passamonti et al. (2009a) Passamonti A., Haskell B., Andersson N., Jones D. I., Hawke I., 2009a, MNRAS, 394, 730
• Passamonti et al. (2009b) Passamonti A., Haskell B., Andersson N., 2009b, MNRAS, 396, 951
• Sawyer (1989) Sawyer R. F., 1989, Phys. Rev. D, 39, 3804
• Shapiro & Teukolsky (1983) Shapiro S. L., Teukolsky S. A., 1983, Black holes, white dwarfs, and neutron stars.  New York, Wiley-Interscience, 1983
• Shibata & Karino (2004) Shibata M., Karino S., 2004, Phys. Rev. D, 70, 084022
• Shternin et al. (2011) Shternin P. S., Yakovlev D. G., Heinke C. O., Ho W. C. G., Patnaude D. J., 2011, MNRAS, 412, L108
• Stergioulas et al. (2004) Stergioulas N., Apostolatos T. A., Font J. A., 2004, MNRAS, 352, 1089
• Thorne (1980) Thorne K. S., 1980, Reviews of Modern Physics, 52, 299
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters