# Neutrino cooling and spin-down of rapidly rotating compact stars

###### Abstract

The gravitational-wave instability of r-modes in rapidly rotating compact stars is believed to spin them down to angular frequencies \Omega\sim 0.1\Omega_{\rm Kepler} soon after their birth in a Supernova. We point out that the r-mode perturbation also impacts the neutrino cooling and viscosity in hot compact stars via processes that restore weak equilibrium. We illustrate this fact with a simple model of spin-down due to gravitational wave emission in compact stars composed entirely of three-flavor degenerate quark matter (a strange quark star). Non-equilibrium neutrino cooling of this oscillating fluid matter is quantified. Our results imply that a consistent treatment of thermal and spin-frequency evolution of a young and hot compact star is a requisite in estimating the persistence of gravitational waves from such a source.

###### pacs:

26.60.+c, 24.85.+p, 97.60.Jd## I Introduction

Neutron stars that undergo pulsations of quadrupole order or higher can emit gravitational waves. Such quasi-normal pulsation modes are typically expected in the aftermath of a supernova or a sudden rearrangement of the neutron star crust or core (egs., a starquarke or a hadron-quark phase transition). The possibility of detecting the associated gravitational radiation in current experiments such as LIGO and VIRGO is enhanced by the existence of the r-mode instability in rotating neutron stars. Discovered by Andersson Andersson:1997xt (), and detailed in subsequent works Friedman:1997uh (); Lindblom:1998wf (), r-modes are fluid pulsations whose dynamics are governed by the Coriolis force, and the l\geq 2 modes grow by emitting gravitational waves in all rotating neutron stars irrespective of their rotation speed. So long as viscous damping effects are small, which is the case for hot young rapidly rotating neutron stars, the r-mode can become unstable to the emission of gravitational waves while the star spins down by losing angular momentum (the Chandrasekhar-Friedman-Schutz or CFS instability Chandra:1970 (); FS:1978 ()). This is one of several reasons why neutron stars are regarded as promising targets of opportunity for gravitational wave detection Schutz:2010xm ().

For experimental detection of gravitational waves (see eg., Schutz:2010xm () for a pedagogical treatment), it is important to have theoretical predictions for the gravitational waveform (eg. \bar{h}(f) in the frequency-domain). According to the simple model discussed in Owen:1998xg () which we employ in this work, the angle-averaged signal in the time domain h(t) depends on the distance, instantaneous angular velocity \Omega(t) and average density of the star, as well as the r-mode amplitude \alpha(t). The time evolution of the mode amplitude and of the angular velocity constitute two coupled first-order differential equations which can be solved explicitly for a particular stellar configuration once viscous damping and gravitational wave timescales are known. Since the viscous damping timescale is very sensitive to the temperature of the fluid, thermal evolution of the star affects r-mode evolution. In previous works Owen:1998xg (); Bildsten:1999zn (); Bondarescu:2007jw (), while this fact has been recognized, the converse effect, viz., that the r-mode perturbation also affects the cooling rate of the neutron star has been missed. A typical r-mode perturbation induces chemically equilibrating weak reactions that alter the neutrino emissivity from the standard equilibrium expressions and hence impact the cooling rate quantitatively. Therefore, adopting a simple power law for the cooling rate based on the temperature-dependence of equilibrium emissivities may not accurately describe the spin-down of the star. Thermal evolution, r-mode evolution and spin-down (i.e, rotational evolution) of a young and hot compact star are all coupled and a consistent analysis of these must be performed in any study of persistent gravitational waveforms from a rapidly rotating compact star.

In this paper, we illustrate this generic interplay of thermochemical and spin-down effects and quantify it in a particular case, viz., that of r-modes in three-flavor (u,d,s) degenerate quark matter. Such a phase can arise in the ultra-dense interior of a neutron star or comprise almost the entire matter in a quark star ^{1}^{1}1Although cold and dense quark matter is believed to be in a color superconducting stateAlford:1998mk (), we will not consider additional facets introduced by Cooper pairing here, leaving such considerations to future work. For r-modes in color superconducting quark matter, see Jaikumar:2008kh (); Sa'd:2008gf (); Andersson:2010sh (); Rupak:2010qg (). We also take into account leptonic (Urca) contributions to the viscosity of strange quark matter, calculated recently Sa'd:2007ud (). Our main findings are: (i) the star can undergo two distinct epochs of gravitational wave emission and rapid spin-down, of very different duration, separated by a decades-long hiatus where its rotation rate is stable;
(ii) the leptonic viscosity can be important for the spin-down history of the star if the initial spin-frequency \nu\lesssim 600 Hz (\nu\approx 0.5\Omega_{\rm Kepler}/(2\pi)); (iii) the star takes about an order of magnitude shorter time to spin down to a stable frequency if non-equilibrium neutrino emissivities are considered and the evolution equations are evolved consistently.

In Section II, we provide a general overview of the link between r-modes, viscosity and gravitational waves, outlining the phenomenological model for mode and spin evolution introduced in Owen:1998xg () and adopted for this work. In Section III, we explain how r-mode perturbations affect the viscosity and neutrino emissivity of the constituent matter. In Section IV, we present results from a coupled analysis of thermal and spin-down evolution of a strange quark star. Section V gathers our conclusions and discusses the implications for persistent gravitational waves from neutron/quark stars. We include some mathematical details in the Appendices.

## II r-mode evolution and gravitational waves

Pulsation modes of compact stars are classified by the nature of the restoring force in the perturbed Euler equation. Analysis in the co-rotating frame shows that the Coriolis force provides the restoring force for small non-radial fluid perturbations Andersson:2000mf (). It is conventional to assume that purely toroidal mode solutions to the perturbed Euler equation exist even for rotating stars, in which case the velocity perturbation can be expressed to {\mathcal{O}}(\Omega) as Lindblom:1998wf ()

\delta\vec{v}=R\Omega f(r/R)\vec{Y}_{lm}^{B}{\rm e}^{i\omega_{\rm rot}^{(l)}t}% \,\,, | (1) |

with R and \Omega being the star’s radius and angular velocity, f and \vec{Y}_{lm}^{B} the radial and angular functional dependence of the mode, and \omega_{\rm rot}^{(l)} the mode frequency in the co-rotating frame. If we restrict ourselves to isentropic perturbations (no composition or temperature gradients), it can be shown that only the l=m modes survive Provost:1981 (), that f(r/R)=\alpha(r/R)^{l} with \alpha an arbitrary constant (dimensionless mode amplitude), and that PP:1978 ()

\omega_{\rm rot}^{(m)}=\frac{2\Omega}{(m+1)}+{\cal O}(\Omega^{3})\,\,. | (2) |

These are the r-modes and they are termed quasi-normal when they lose energy to gravitational waves (l\geq 2). The energy of the r-mode is dissipated according to Lindblom:1998wf ():

\displaystyle\frac{dE}{dt}=-(\omega_{\rm rot}-m\Omega)^{2m+1}\omega_{\rm rot}|% \delta J_{mm}|^{2}-\int d^{3}r(2\eta\delta\sigma^{ab}\delta\sigma^{*}_{ab}+% \zeta\delta\sigma\delta\sigma^{*})\,\,, | (3) |

where \delta\sigma=\nabla_{a}\delta v^{a} is the volume expansion due to the r-mode and \delta\sigma_{ab}=\frac{1}{2}(\nabla_{a}\delta v_{b}+\nabla_{b}\delta v_{a}-% \frac{2}{3}\delta_{ab}\nabla_{c}\delta v^{c}) is (the traceless part of) the shear tensor. The first term is the energy radiated in gravitational waves to lowest order in \Omega with \delta J_{mm} being the current multipole

\displaystyle\delta J_{mm}\propto R^{2}\Omega\int_{0}^{R}dr\rho\left(\frac{r^{% 2}}{R}\right)^{(m+1)}\,\,. | (4) |

For m\geq 2, \omega_{\rm rot}<m\Omega, so that the r-mode energy grows with gravitational wave emission, triggering the r-mode instability. Viscosity from a variety of mechanisms in the fluid can suppress r-mode growth, and counter the instability. The fate of the r-mode is determined by a competition between the timescale for viscous damping and gravitational wave emission. According to the simple model for mode growth in the linear regime described in Owen:1998xg (), the evolution equations for \alpha (mode amplitude) and \Omega (angular velocity) are:

\displaystyle\dot{\Omega} | \displaystyle= | \displaystyle-\frac{2\Omega}{\tau_{V}}\frac{\alpha^{2}Q}{1+\alpha^{2}Q}\,\,, | (5) | ||

\displaystyle\dot{\alpha} | \displaystyle= | \displaystyle-\frac{\alpha}{\tau_{GR}}-\frac{\alpha}{\tau_{V}}\frac{1-\alpha^{% 2}Q}{1+\alpha^{2}Q}\,\,, |

where the viscous damping timescale is \tau_{V}={(\tau_{\zeta}^{-1}+\tau_{\eta}^{-1})}^{-1} (\zeta: bulk viscosity, \eta: shear viscosity) and the gravitational wave growth timescale is \tau_{GR}. In general, \tau_{i}^{-1}=-(dE_{i}/dt)/2E;\,i=\zeta,\eta,GR. Reading off (dE_{i}/dt) from Eqn. (3), and using the lowest order expression for the r-mode energy E, one obtains the expressions for the damping timescales \tau provided in Appendix B. The quantity Q=3\bar{J}/2\bar{I} where Owen:1998xg ()

\bar{J}=\frac{1}{MR^{4}}\int_{0}^{R}\rho r^{6}dr;\quad\bar{I}=\frac{8\pi}{3MR^% {2}}\int_{0}^{R}\rho r^{4}dr\,\,, | (6) |

with M being the stellar gravitational mass and \rho the energy density. In this model Owen:1998xg (), exponential mode growth is followed by non-linear saturation characterized by \alpha^{2}\sim\alpha_{0}^{2}\sim 1 where \dot{\alpha}\approx 0. While such large saturation amplitudes are probably unrealistic due to multi-mode coupling and onset of non-linear effects Arras:2002dw (); Bondarescu:2007jw (); Bondarescu:2008qx (), we prefer to use this large value to demonstrate the qualitative effect induced by the coupling of temperature and mode evolution. Ultimately, our analysis should be applied to more realistic saturation amplitudes and include the coupling to daughter modes Bondarescu:2007jw (); Bondarescu:2008qx (). Within our working model Owen:1998xg (), once saturation is achieved, only the angular velocity evolves

\dot{\Omega}=\frac{2\Omega}{\tau_{GR}}\frac{\alpha_{0}^{2}Q}{1-\alpha_{0}^{2}Q% }\,\,. | (7) |

It is important to note that \tau_{V} is strongly temperature dependent while \tau_{GR} is temperature-independent so that solution of Eqn. (5) requires an estimate of the thermal evolution T(t), which is determined in the simplest case (ignoring thermal gradients and conduction) by

C_{v}\frac{dT}{dt}=-L_{\nu}+L_{h}\,\,, | (8) |

where C_{v} is the specific heat, L_{\nu} is the neutrino luminosity and L_{h} is the heating rate due to dissipation ^{2}^{2}2Viscous heating by shear viscosity, estimated as in Owen:1998xg (), has a small effect on the temperature evolution in case of strange quark matter, while local heating due to bulk viscosity is rapidly radiated away in neutrinos and is unimportant for even the smallest timescales required for evolving eqs (5). Therefore, we only include viscous heating from shear viscosity in L_{h}.. While the r-mode is growing, the local baryonic number density fluctuates according to

n_{B}(t)=n_{B}^{0}+\delta n_{B}(t)=n_{B}^{0}+\delta\bar{n}_{B}{\rm Re}({\rm e}% ^{i\omega t})\,\,, | (9) |

where \omega\equiv\omega_{\rm rot}^{(2)} and \delta\bar{n}_{B} (proportional to \alpha, see Eqn. (12)) is the magnitude of the fluctuation. These fluctuations, which first appear at {\cal O}(\Omega^{2}) rather than {\cal O}(\Omega), destroy weak equilibrium, considerably enhancing the neutrino emissivity above its equilibrium value. They also affect the bulk viscosity, so that \zeta and L_{\nu} must be evaluated in the fluctuating foreground created by the r-mode. In the next section, we describe this computation, highlighting the effect of the r-mode perturbation on \zeta and L_{\nu}.

## III r-modes, Viscosity & Non-equilibrium neutrino rates

The r-mode perturbation described by the density fluctuation Eqn. (12) causes a change in the fluid volume as well as the chemical potentials of the electron (e) and the up (u), down (d) and strange (s) quarks. Since weak equilibrating processes are relatively slow, the pressure fluctuation, expressed as \delta P=\delta P_{0}{\rm e}^{i\omega t}, lags the density fluctuation through the former’s dependence on the fluctuating chemical potentials (i.e., \delta P_{0} acquires a phase). The dependence on temperature fluctuations is ignored since, as mentioned before, dissipated heat can be carried away much faster than weak equilibration timescales. Both leptonic and non-leptonic processes in dense quark matter can contribute to the viscosity at temperatures T\geq 10^{10}K for typical r-mode frequencies \omega\sim 1kHz Sa'd:2007ud (). Following the approach of Sa'd:2007ud (), it is useful to choose the two independent fluctuations in the chemical potentials corresponding to the leptonic processes as \delta\mu_{1}=\delta(\mu_{d}-\mu_{u}-\mu_{e}) and \delta\mu_{2}=\delta(\mu_{s}-\mu_{u}-\mu_{e}) which can be expressed in terms of variations in the following quantities: baryon number (\delta n_{B}/n_{B}), charge (\delta X_{e}) and strangeness (\delta X_{s})

\delta\mu_{i}=A_{i}\delta X_{s}+B_{i}\delta X_{e}+C_{i}\left(\frac{\delta n_{B% }}{n_{B}}\right);\quad X_{s}=\frac{n_{s}}{n_{B}}~{}\,,\quad X_{e}=\frac{n_{e}}% {n_{B}}\,\,, | (10) |

where n_{s} and n_{e} are the strange quark and electron number densities. The non-leptonic process d+u\leftrightarrow u+s is characterized by \delta\mu_{3}=\delta\mu_{1}-\delta\mu_{2}. The bulk viscosity, which is a measure of the power dissipated per unit volume over an oscillation cycle is given by

\zeta=-\frac{1}{\omega}{\rm Im}(\delta P_{0})\frac{n_{B}^{0}}{\delta\bar{n}_{B% }};\quad{\rm Im}(\delta P_{0})=-n_{B}\left[C_{1}{\rm Im}(\delta X_{e,0})+(C_{2% }-C_{1}){\rm Im}(\delta X_{s,0})\right]\,\,. | (11) |

For the density perturbation, we adopt an expression derived in Lindblom:1998ka () (also Eqn. (6.8) of Lindblom:1999yk ())

\displaystyle\frac{\delta\bar{n}_{B}}{n_{B}} | \displaystyle= | \displaystyle R^{2}\Omega^{2}\frac{d\rho}{dP}\left[\delta U_{0}+\delta\Phi_{0}% \right]\,\,, | (12) | ||

\displaystyle\delta U_{0} | \displaystyle= | \displaystyle\alpha\left(\frac{r}{R}\right)^{m+1}P_{m+1}^{m}({\rm cos}~{}% \theta){\rm e}^{im\phi};\quad\delta\Phi_{0}=\alpha\delta\phi_{0}(R)P_{m+1}^{m}% ({\rm cos}~{}\theta){\rm e}^{im\phi}\,\,, |

where \delta\Phi_{0} is the change in gravitational potential due to the r-mode oscillation (\delta\phi_{0}(R) is obtained from Eq.(4.4) of Lindblom:1999yk ()). The perturbation, and hence \delta\mu_{i}, has a radial and angular dependence, leading to different emissivities and viscosities in different directions and radial distances. We take a much simplified approach whereby, as is common in interpretation of neutron star thermal observations, the temperature T refers to or is related to a local surface temperature map (eg., a ”hot spot”) Page:1994dz (). Even if full angular dependence is retained, additional anisotropy in the bulk viscosity can be induced by strong magnetic fields in quark matter Huang:2009ue (). For now, we proceed with our simplifying assumptions of homogeneity and isotropy with magnetic fields set to zero. These assumptions should be relaxed in fully 3-dimensional calculations, which require extensive computations on a spatial grid. The computation of {\rm Im}(\delta X_{e,o}),{\rm Im}(\delta X_{s,0}) and the co-efficients A_{i},B_{i},C_{i} is straightforward once an equation of state is specified (see Appendix A for the equation of state for degenerate strange quark matter adopted in this work). We find

\displaystyle\zeta | \displaystyle= | \displaystyle\frac{n_{B}}{\omega}\frac{1}{g_{1}^{2}+g_{2}^{2}}\sum_{(i,j)=1;i<% j}^{3}\alpha_{i}\alpha_{j}\theta_{ij}\,\,, | (13) | ||

\displaystyle g_{1} | \displaystyle= | \displaystyle(A_{1}B_{2}-A_{2}B_{1})\Sigma-\Pi;\,\,\,g_{2}=\alpha_{1}\alpha_{2% }(A_{2}-A_{1})+\alpha_{1}\alpha_{3}(A_{2}-B_{2})-\alpha_{2}\alpha_{3}B_{1}\,\,, | |||

\displaystyle\Sigma | \displaystyle= | \displaystyle\alpha_{1}+\alpha_{2}+\alpha_{3};\,\Pi=\alpha_{1}\alpha_{2}\alpha% _{3}\,\,, | |||

\displaystyle\Theta_{12} | \displaystyle= | \displaystyle\Sigma(A_{1}C_{2}-A_{2}C_{1})^{2}+\Pi(C_{1}-C_{2})^{2}\,\,,\,% \Theta_{13}=\Sigma\left[C_{1}(A_{2}-B_{2})^{2}+C_{2}B_{2}\right]^{2}+\Pi C_{2}% ^{2}\,\,,\,\Theta_{23}=\Sigma(C_{1}B_{2}-C_{2}B_{1})^{2}+\Pi C_{2}^{2}\,\,, | (14) |

where A_{i},B_{i},C_{i} are given in Appendix A. In Eqn. (13), \alpha_{i}=n_{B}\omega/\lambda_{i} where the \lambda_{i}s are obtained from the decay rates (\Gamma) of the non-equilibrium processes at first order as

\displaystyle\Gamma_{\rm net}^{(1)} | \displaystyle= | \displaystyle\Gamma_{d\rightarrow u+e^{-}+\bar{\nu}_{e}}-\Gamma_{u+e^{-}% \rightarrow d+\nu_{e}}\approx\lambda_{1}\delta\mu_{1}\,\,, | (15) | ||

\displaystyle\Gamma_{\rm net}^{(2)} | \displaystyle= | \displaystyle\Gamma_{s\rightarrow u+e^{-}+\bar{\nu}_{e}}-\Gamma_{u+e^{-}% \rightarrow s+\nu_{e}}\approx\lambda_{2}\delta\mu_{2}\,\,, | |||

\displaystyle\Gamma_{\rm net}^{(3)} | \displaystyle= | \displaystyle\Gamma_{u+d\rightarrow s+u}-\Gamma_{s+u\rightarrow u+d}\approx% \lambda_{3}(\delta\mu_{1}-\delta\mu_{2})=\lambda_{3}\delta\mu_{3}\,\,, |

where the second relation in each line is a commonly made approximation to first order in \delta\mu. The net decay rates for the leptonic processes can be computed exactly from the corresponding non-equilibrium neutrino rates FloresTulian:2006fq ()

\frac{\partial\epsilon_{\rm tot}^{(i)}}{\partial\delta\mu_{i}}=3(\Gamma_{f}^{(% i)}-\Gamma_{b}^{(i)})=3\Gamma_{\rm net}^{(i)};\quad i=1,2 | (16) |

where \Gamma_{f} and \Gamma_{b} are the forward and backward reaction rates respectively, and \epsilon_{\rm tot}^{(i)}=\epsilon_{f}^{(i)}+\epsilon_{b}^{(i)}=\epsilon_{\bar{% \nu}_{e}}^{(i)}+\epsilon_{\nu_{e}}^{(i)}. As shown in Appendix B, the relevant non-equilibrium neutrino rates \epsilon_{\rm tot}^{(i)} are simply related to the equilibrium emissivities (\epsilon_{\rm tot}^{(i),{\rm eq}}) as Reise:1995 ()

\epsilon_{\rm tot}^{(i)}=\epsilon_{\rm tot}^{(i),{\rm eq}}\left[1+\frac{1071u_% {i}^{2}+315u_{i}^{4}+21u_{i}^{6}}{457}\right];\quad u_{i}=\frac{\delta\mu_{i}}% {\pi kT}\,\,. | (17) |

Combining Eqns.(15),(16) and (17) with Eqn. (21) generates a self-consistent solution for the \delta\mu_{i} even for large u_{i}. Since it turns out that \delta\mu_{i}\gg kT can occur for the r-mode perturbation, we retain the full expression in Eqn. (17). Inclusion of these ”large-amplitude” rates has been recently emphasized in the context of r-modes Alford:2010gw (). The decay rate for the non-leptonic process \Gamma_{\rm net}^{(3)} has been calculated in various works. We adopt the analytic result of Madsen Madsen:1992sx (); Madsen:1993xx (), which improves Sawyer’s result Sawyer:1989 () by including a \delta\mu_{3}^{3} term as well as m_{s}\neq 0 corrections (their inclusion increases the rate). This approximation matches quite well with the exact numerical result up to T\sim 10^{11}K, which is the upper temperature bound considered in this work. The equilibrium neutrino emissivities appearing in Eqn. (17) are Iwamoto ()

\displaystyle\epsilon_{\rm tot}^{(1),{\rm eq}} | \displaystyle= | \displaystyle\frac{457}{630}\alpha_{s}T^{6}G_{F}^{2}{\rm cos}^{2}~{}\theta_{c}% \mu_{e}\mu_{d}\mu_{u}\,\,, | (18) | ||

\displaystyle\epsilon_{\rm tot}^{(2),{\rm eq}} | \displaystyle= | \displaystyle\frac{457\pi}{1680}T^{6}G_{F}^{2}{\rm sin}^{2}~{}\theta_{c}\mu_{s% }m_{s}^{2}\,\,. | |||

## IV Neutron star cooling and spin-down

As mentioned in the previous section, \delta\mu_{i} can be much larger than kT over some time domain of the evolution. In our simulation, this happens in the initial few hundreds of seconds of the star’s spin down due to the r-mode instability. Consequently, non-equilibrium neutrino emissivities are much higher than equilibrium values, and we should anticipate a change in the cooling history of the star. In Fig 1 (bottom panel), we display results for the angular velocity versus temperature evolution from the combined analysis of Eqn.(5), for which we used the bulk and shear viscosity damping timescales displayed in the top panel of the same figure. Cooling with equilibrium and non-equilibrium neutrino emissivities (Eqn. (17)) is studied.

Firstly, we note that the double-valley feature of the damping timescale (top panel) is reflected in the double-peak feature of the critical frequency curve, which is the locus of points satisfying 1/\tau_{\rm f}(T)|_{\Omega_{c}}=[1/\tau_{\zeta}(T)+1/\tau_{\eta}(T)+1/\tau_{GW% }(T))]_{\Omega_{c}}=0. The critical frequency curve marks the boundary between rapid spin-down and unstable mode growth accompanied by gravitational waves in the region above the curve and stable rotation with a decaying r-mode amplitude below it. As expected, the critical frequency curve mirrors the behaviour of the bulk viscosity of strange quark matter (see eg., Fig. 3 of Sa'd:2007ud ()). The larger of the two peaks corresponds to the maximum in the bulk viscosity of the non-leptonic process Madsen:1992sx () while the smaller peak is from the maximum in the leptonic process Sa'd:2007ud (). The peak height is inversely proportional to the r-mode frequency (\omega=0.88kHz for this figure). For the simple model assumed here, these features determine the evolution of the angular velocity with temperature, as the star cools and spins down.

To determine the persistence of gravitational waves from r-modes in a cooling compact star, we track the angular momentum evolution against the critical frequency curve to see when the star enters/exits the stable region. In the case of equilibrium neutrino rates, which ignores the effect of the r-mode on neutrino emissivities, we find that the star rapidly spins down in about a few hours (\approx 10^{4} seconds) from its initial spin frequency, assumed to be the Kepler frequency (\Omega_{K}/(2\pi)\approx 1200 Hz) down to about a third of the starting value. This is the segment OA in Fig. 1 (bottom panel). The time to spin down to a stable rotation rate does not depend on the initial mode amplitude (typically, we varied this between 10^{-6} and 10^{-3}), nor does it depend on the mode saturation amplitude \alpha_{0}, which we varied between 0.3-1.0. However, the stable frequency at which the star settles down upon entering the stable region depends on \alpha_{0}, becoming smaller with increasing \alpha_{0}. Once the star enters this stable region, where bulk viscosity is large enough to damp the r-mode, the mode amplitude decreases to negligible values. Since the r-mode is damped, it does not affect the neutrino rates in any way and equilibrium rates apply. While a small time step of about 0.01 seconds is required in the initial rapid cooling phase, the time step in the stable region can be safely increased to hundreds of seconds, since temperature evolution and angular velocity evolution is extremely slow. The star takes of order 10^{9} seconds (nearly 100 years) to cool down to the point where it can once again enter the instability region. This hiatus between spin-down epochs is depicted by the segment AA^{\prime}. In the interim, the star has cooled by two orders of magnitude down to T\approx 0.7 keV, but hardly changed its spin-frequency (here, we do not take into account spin-down due to dipole braking or the possibility of glitches). Upon entering the unstable region at A^{\prime}, the mode amplitude begins a rapid growth from a negligibly small value and saturates within a year. At this point, the star spins down again due to the CFS instability and takes {\cal O}(10^{8}) seconds (10 years) to reach a stable value of \Omega/(2\pi)\sim 100 Hz. During this phase, the star is expected to once again radiate gravitational waves. Thus, a strange quark star might be expected to become unstable to r-modes over two distinct epochs, separated by about a 100 years, with the second one lasting about 10,000 times longer than the first.

For a consistent and more accurate estimate of the persistence of gravitational waves, non-equilibrium neutrino rates must be taken into account. A clear difference can be seen in the angular momentum history in Fig. 1. Although we begin with the same initial conditions, neutrino cooling is more rapid, since the emissivities are larger (Eqn. (17)). Consequently, the star enters the stable region at a higher angular velocity and in less time than the equilibrium case (segment OB). We find this time to be of order 10^{3} seconds (less than 1 hour). The cooling in the stable region is once again marked by slow cooling and a negligible r-mode amplitude, resulting in equilibrium-type cooling. However, the fact that the star is at a higher stable rotation rate now implies that it spends less time in the stable region (segment BB^{\prime}). This lasts approximately 5.10^{8} seconds (\approx 15 years), to be contrasted with the equilibrium case, where this era lasts for about 100 years. Upon exiting this region, the star spins down rapidly once more, taking \approx 10^{7} seconds (4 months) to reach its final stable value. This is the segment B^{\prime}O^{\prime}, shown in enlarged inset view in Fig. 1. The final rotation rate is exactly the same irrespective of the cooling history. Therefore, the importance of incorporating non-equilibrium emissivities is only evident in the first 100 years or so of the star’s spin-down history. Another important difference is the shorter time for which the star woud be active in gravitational waves if non-equilibrium rates are used. This indicates that the epoch of gravitational waves is shorter than expected, narrowing the window of detection. Similar results are expected for the case of neutron stars on the basis of the generality of the evolution equations, even though neutron matter has very different viscosities and cooling rates.

## V Conclusions

We have shown that the r-mode perturbation affects the neutrino cooling and viscosity in hot compact stars via processes that restore weak equilibrium. Using a simplified model of spin-down due to gravitational wave emission in compact stars composed entirely of three-flavor degenerate quark matter, we find that non-equilibrium neutrino rates associated to the r-mode fluctuation change considerably the cooling history of the star, as well as the angular velocity evolution, resulting in shorter epochs of gravitational wave emission for isolated quark stars. We have taken into account a recent calculation of leptonic contributions to the viscosity of quark matter, but find that it can only play a role if the initial spin frequency is chosen sufficiently small (\approx 600Hz, see Fig. 1). This is because of the smaller viscosity of the leptonic process, specially at r-mode frequencies characteristic of very rapidly rotating stars. However, the role of the non-leptonic process is very prominent. It leads to two distinct rapid spin-down events in the star’s evolution, caused by the CFS instability towards gravitational wave emission. Although we have focused on quark matter here as an exemplar, it should be clear that a similar exercise can be performed for any sufficiently simple phase of dense matter, including beta-stable neutron-rich matter. It is important for r-mode saturation and damping studies that one includes the large-amplitude bulk viscosity results, as argued in Alford:2010gw (), as well as realistic saturation amplitudes Bondarescu:2007jw (). While analytic approximations for the bulk viscosity exist for a simple quark matter equation of state like the one adopted in this work, a full numerical solution is required for hadronic matter, where non-linearities in the chemical equilibration rates can take more complex forms. Furthermore, non-equilibrium rates Haensel:1992 () and updated viscosities Benhar:2007yj () need to be taken into account in performing a similar calculation with hadronic equations of state. Despite these complications, one can predict that for isolated compact stars made up of neutron-rich matter, since there is no stable region between \approx 10^{7}-10^{10}K at rapid rotation rates Owen:1998xg (); Jaikumar:2008kh (), non-equilibrium cooling is probably important only over the first year or more of cooling.

We also mention some caveats that follow from the simplifications made in this work. The lowest order density fluctuation (Eqn. (12)) enters at order \Omega^{2}, while the expressions for the r-mode frequency (Eqn.(2)) and the current multipole (Eqn. (4)) are {\cal O}(\Omega). Higher-order effects change the mode frequency by less than 10% (see eg., Table 1 in Jaikumar:2008kh ()) while Eqn… (4) suffices for gravitational radiation from the current multipole. A more detailed discussion of higher-order rotation effects is given in Lindblom:1999yk (). Bulk viscosity timescales used in this work are approximate, and can be improved by replacing the Eulerian expression for the density perturbation, Eqn. (12) by a Lagrangian one Lindblom:1999yk (). This can change the bulk viscosity damping timescale by an order of magnitude, but the critical frequency curve is hardly modified because the bulk viscosity has a steep temperature dependence. Thus, the results presented in the bottom panel of Fig. 1 are accurate in this respect, even though the bulk viscosity timescale, as shown in the top panel is only approximate.

Non-equilibrium cooling can also be important for accreting neutron stars, which are believed to be persistent sources of gravitational waves Andersson:1998qs (), undergoing spin-cycles as they get spun-up by accretion torques and spun-down repeatedly by the r-mode instability. Larger duty-cycles are expected for exotic phases, such as quark matter or hyperonic matter Nayyar:2005th (), providing the motivation to extend our work to low-mass X-ray binaries. Even for isolated strange quark stars, Andersson has shown Andersson:2001ev () that if a strange star enters the instability window spinning near the break-up limit, it will emit bursts of gravitational waves as the spin evolution straddles the critical frequency curve. We plan to explore this interesting case as well as the spin-cycles of LMXBs including non-equilibrium cooling rates and improved quark matter viscosites. Finally, we note that LIGO data sets Abadie:2010hv () can already be adapted to search for imprints of the r-mode instability in gravitational waves Owen:2010ng (). In anticipation of the improved sensitivity of advanced LIGO, it is important to continue to examine various factors affecting the r-mode evolution.

## Acknowledgments

We thank C. D. Roberts, T. Klähn, B. El-Bennich and R. Young for useful discussions. This work was supported in part by the U.S. Department of Energy, Office of Nuclear Physics, contract no. DE-AC02-06CH11357 and by the United States ARMY High Performance Computing Research Center’s Research and Infrastructure Development Program.

## Appendix A

With reference to Eqn. (13), the definitions

\displaystyle e_{1} | \displaystyle= | \displaystyle\Sigma(B_{1}C_{2}-B_{2}C_{1});\quad e_{2}=(\alpha_{1}\alpha_{2}C_% {1}-\alpha_{1}(\alpha_{2}+\alpha_{3})C_{2})\,\,, | (19) | ||

\displaystyle f_{1} | \displaystyle= | \displaystyle\Sigma(C_{1}A_{2}-A_{1}C_{2});\quad f_{2}=(\alpha_{1}\alpha_{3}C_% {2}+\alpha_{2}\alpha_{3}C_{1})\,\,, | |||

\displaystyle g_{1} | \displaystyle= | \displaystyle\Sigma(A_{1}B_{2}-B_{1}A_{2})-\Pi;\quad g_{2}=\alpha_{1}\alpha_{2% }(A_{2}-A_{1})+\alpha_{1}\alpha_{3}(A_{2}-B_{2})C_{1}-\alpha_{2}\alpha_{3}B_{1% }\,\,, |

lead to the expressions for \delta X_{s,0} and \delta X_{e,0}

\delta X_{s,0}=\frac{e_{1}+ie_{2}}{g_{1}+ig_{2}}\left(\frac{\delta\bar{n}_{B}}% {n_{B}}\right);\quad\delta X_{e,0}=\frac{f_{1}+if_{2}}{g_{1}+ig_{2}}\left(% \frac{\delta\bar{n}_{B}}{n_{B}}\right)\,\,, | (20) |

from which it follows that

\delta\mu_{i}(t)=\left(\frac{\delta\bar{n}_{B}}{n_{B}}\right)\left[C_{i}{\rm cos% }~{}(\omega t)+\sqrt{\frac{e_{1}^{2}+e_{2}^{2}}{g_{1}^{2}+g_{2}^{2}}}A_{i}{\rm cos% }~{}(\omega t-\phi_{e})+\sqrt{\frac{f_{1}^{2}+f_{2}^{2}}{g_{1}^{2}+g_{2}^{2}}}% B_{i}{\rm cos}~{}(\omega t-\phi_{f})\right] | (21) |

with

{\rm tan}\phi_{e}=\frac{e_{1}g_{2}-e_{2}g_{1}}{e_{1}g_{1}+e_{2}g_{2}};\quad{% \rm tan}\phi_{f}=\frac{f_{1}g_{2}-f_{2}g_{1}}{f_{1}g_{1}+f_{2}g_{2}}\,\,. | (22) |

Given the following thermodynamic potential for 3-flavor degenerate quark matter (with electrons)

\displaystyle\Omega | \displaystyle= | \displaystyle\sum_{u,d,s,e}\Omega_{i};\quad\Omega_{e}=-\frac{\mu_{e}^{4}}{12% \pi^{2}}\,\,, | (23) | ||

\displaystyle\Omega_{f} | \displaystyle= | \displaystyle-\frac{\mu_{f}^{4}}{4\pi^{2}}\left\{\sqrt{1-z_{f}^{2}}\left(1-% \frac{5}{2}z_{f}^{2}\right)+\frac{3}{2}z_{f}^{4}{\rm ln}\left[F(z)\right]-% \frac{2\alpha_{s}}{\pi}\left[3\left(z_{f}^{2}{\rm ln}\left[F(z)\right]-\sqrt{1% -z_{f}^{2}}\right)^{2}-2(1-z_{f}^{2})^{2}\right]\right\}\,\,, | |||

\displaystyle F(z) | \displaystyle= | \displaystyle\frac{1+\sqrt{1-z_{f}^{2}}}{z_{f}};\quad z_{f}=m_{f}/\mu_{f};% \quad\alpha_{s}=g_{s}^{2}/4\pi\,\,, |

where f denotes flavor, the coefficients A_{i},B_{i},C_{i} follow from

\displaystyle A_{1} | \displaystyle= | \displaystyle-n_{B}\frac{\partial\mu_{d}}{\partial n_{d}};\quad B_{1}=-n_{B}% \left[\frac{\partial\mu_{d}}{\partial n_{d}}+\frac{\partial\mu_{u}}{\partial n% _{u}}+\frac{\partial\mu_{e}}{\partial n_{e}}\right];\quad C_{1}=\left[n_{s}% \frac{\partial\mu_{s}}{\partial n_{s}}-n_{d}\frac{\partial\mu_{d}}{\partial n_% {d}}\right]\,\,, | (24) | ||

\displaystyle A_{2} | \displaystyle= | \displaystyle n_{B}\frac{\partial\mu_{s}}{\partial n_{s}};\quad B_{2}=-n_{B}% \left[\frac{\partial\mu_{u}}{\partial n_{u}}+\frac{\partial\mu_{e}}{\partial n% _{e}}\right];\quad C_{2}=\left[n_{s}\frac{\partial\mu_{s}}{\partial n_{s}}-n_{% u}\frac{\partial\mu_{u}}{\partial n_{u}}-n_{e}\frac{\partial\mu_{e}}{\partial n% _{e}}\right]\,\,, |

upon using n_{i}=-\partial\Omega_{i}/\partial\mu_{i} and the expressions

n_{f}\frac{\partial\mu_{f}}{\partial n_{f}}=\frac{\mu_{f}}{3}\left\{1-z_{f}^{2% }+\frac{4\alpha_{s}}{\pi}z_{f}^{2}\left[{\rm ln}(2z_{f}^{-1})-\frac{2}{3}% \right]\right\};\quad n_{e}\frac{\partial\mu_{e}}{\partial n_{e}}=\frac{\mu_{e% }}{3}\,\,. | (25) |

## Appendix B

Neutrino emissivity: Omitting volume normalizations which ultimately cancel, the emissivity is given by

\epsilon_{\bar{\nu}_{e}}=\int\biggl{(}\,\prod_{i}\frac{d^{3}{\bm{p}}_{i}}{2E_{% i}(2\pi)^{3}}\biggr{)}\,\frac{d^{3}{\bm{k}}}{2\omega_{k}(2\pi)^{3}}\,\,(2\pi)^% {4}\,\delta(E_{f}^{\rm tot}-E_{i}^{\rm tot})\,\delta^{3}({\bm{P}}_{f}-{\bm{P}}% _{i})\,\biggl{(}\>\sum_{\text{spin}}\,\bigl{|}\,{\mathcal{M}}\bigr{|}^{2}% \biggr{)}\,\omega_{k}\,{\mathcal{F}}(E), | (26) |

where ({\bm{p}}_{i},E_{i}) denote the momenta and energy of the incoming and outgoing quarks and (\omega_{k},{\bm{k}}) label the neutrino energy and momenta. The delta functions account for total energy and momentum conservation, and the function {\cal{F}}(E)=\prod_{\rm in}f(E_{i})\prod_{\rm out}(1-f(E_{i})) is the product of occupation (incoming) and blocking (outgoing) statistical factors for the quarks, with Fermi-Dirac distribution function f(E_{i})=(\exp(\beta E_{i})+1)^{-1} (the neutrinos free stream). The matrix element {\mathcal{M}} is calculated from the corresponding tree level Feynman diagrams (see Sa'd:2007ud ()) in the low-energy limit. It can be shown that for the non-equilibrium case (E_{i}\approx\mu_{i}+\delta\mu_{i}), the neutrino energy integral in the total emissivity expression is modified from the equilibrium case Iwamoto () such that

\epsilon_{\bar{\nu}_{e}}+\epsilon_{\nu_{e}}\propto\int~{}d\epsilon_{\nu}% \epsilon_{\nu}^{3}\left[J\left(\frac{\epsilon_{\nu}+\delta\mu}{kT}\right)+J% \left(\frac{\epsilon_{\nu}-\delta\mu}{kT}\right)\right];\quad J(x)=\frac{T^{2}% }{2}\left[\frac{\pi^{2}+x^{2}}{1+{\rm e}^{x}}\right] | (27) |

yielding Eqn. (17).

Damping timescales: For viscous damping of the r-mode, the shear viscosity damping timescale is

\frac{1}{\tau_{\eta}}=\frac{(m-1)(2m+1)}{\int_{0}^{R}dr\rho r^{2m+2}}\int_{0}^% {R}dr\eta r^{2m}, | (28) |

where the shear viscosity is taken from Heiselberg:1993cr (), while the bulk viscosity damping timescale follows from

\displaystyle\frac{1}{\tau_{\zeta}} | \displaystyle= | \displaystyle-\frac{1}{2E}\left(\frac{dE}{dt}\right)_{\zeta}\,\,, | (29) | ||

\displaystyle E | \displaystyle= | \displaystyle\frac{\pi}{2m}\left(m+1\right)^{3}\left(2m+1\right)!R^{4}\Omega^{% 2}\int_{0}^{R}dr\rho\left(\frac{r}{R}\right)^{2m+2}\,\,. | |||

The gravitational radiation growth time scale is Friedman:1997uh ()

\frac{1}{\tau_{GW}}=-\frac{32\pi G\Omega^{2m+2}}{c^{2m+3}}\frac{\left(m-1% \right)^{2m}}{\left[\left(2m+1\right)!!\right]^{2}}\left(\frac{m+2}{m+1}\right% )^{2m+2}\int_{0}^{R}dr\rho r^{2m+2}\,\,. | (30) |

## References

- (1) N. Andersson, Astrophys. J. 502, 708 (1998) [arXiv:gr-qc/9706075].
- (2) J. L. Friedman and S. M. Morsink, Astrophys. J. 502, 714 (1998) [arXiv:gr-qc/9706073].
- (3) L. Lindblom, B. J. Owen and S. M. Morsink, Phys. Rev. Lett. 80, 4843 (1998) [arXiv:gr-qc/9803053].
- (4) S. Chandrasekhar, Phys. Rev. Lett. 24, 611 (1970).
- (5) J. L. Friedman and B. F. Schutz, Astrophys. J. 222, 881 (1978).
- (6) B. F. Schutz and F. Ricci, arXiv:1005.4735 [gr-qc].
- (7) B. J. Owen, L. Lindblom, C. Cutler, B. F. Schutz, A. Vecchio and N. Andersson,Phys. Rev. D 58, 084020 (1998) [arXiv:gr-qc/9804044].
- (8) L. Bildsten and G. Ushomirsky, Astrophys. J. 529, Issue 1, L33 (2000) [arXiv:astro-ph/9911155].
- (9) R. Bondarescu, S. A. Teukolsky and I. Wasserman, Phys. Rev. D 76, 064019 (2007) [arXiv:0704.0799 [astro-ph]].
- (10) M. G. Alford, K. Rajagopal and F. Wilczek, Nucl. Phys. B 537, 443 (1999) [arXiv:hep-ph/9804403].
- (11) P. Jaikumar, G. Rupak and A. W. Steiner, Phys. Rev. D 78, 123007 (2008) [arXiv:0806.1005 [nucl-th]].
- (12) N. Andersson, B. Haskell and G. L. Comer, Phys. Rev. D 82, 023007 (2010) [arXiv:1005.1163 [astro-ph.SR]].
- (13) G. Rupak and P. Jaikumar, arXiv:1005.4161 [nucl-th].
- (14) B. A. Sa’d, arXiv:0806.3359 [astro-ph].
- (15) B. A. Sa’d, I. A. Shovkovy and D. H. Rischke, Phys. Rev. D 75, 125004 (2007) [arXiv:astro-ph/0703016].
- (16) N. Andersson and K. D. Kokkotas, Int. J. Mod. Phys. D 10, 381 (2001) [arXiv:gr-qc/0010102].
- (17) J. Provost, G. Berthomieu, and A. Rocca, Astron. Astrophys. 94, 126, (1981).
- (18) J. Papaloizou and J. E. Pringle, Mon. Not. R. Astron. Soc. 182, 423 (1978).
- (19) P. Arras, E. E. Flanagan, S. M. Morsink, A. K. Schenk, S. A. Teukolsky and I. Wasserman, Astrophys. J. 591, 1129 (2003) [arXiv:astro-ph/0202345].
- (20) R. Bondarescu, S. A. Teukolsky and I. Wasserman, Phys. Rev. D 79, 104003 (2009) [arXiv:0809.3448 [astro-ph]].
- (21) L. Lindblom and J. R. Ipser, Phys. Rev. D 59, 044009 (1999) [arXiv:gr-qc/9807049].
- (22) L. Lindblom, G. Mendell and B. J. Owen, Phys. Rev. D 60, 064006 (1999) [arXiv:gr-qc/9902052].
- (23) D. Page, Astrophys. J. 442, 273 (1995) [arXiv:astro-ph/9407015].
- (24) X. G. Huang, M. Huang, D. H. Rischke and A. Sedrakian, Phys. Rev. D 81, 045015 (2010) [arXiv:0910.3633 [astro-ph.HE]].
- (25) S. Flores-Tulian and A. Reisenegger, Mon. Not. Roy. Astron. Soc. 372, 276 (2006) [arXiv:astro-ph/0606412].
- (26) A. Reisenegger, Astrophys. J. 442, 749 (1995).
- (27) M. G. Alford, S. Mahmoodifar and K. Schwenzer, arXiv:1005.3769 [nucl-th].
- (28) J. Madsen, Phys. Rev. D 46, 3290 (1992).
- (29) J. Madsen, Phys. Rev. D 47, 325 (1993).
- (30) R. F. Sawyer, Phys. Lett. B233, 412 (1989).
- (31) N. Iwamoto, Ann. Phys. 141, 1 (1982).
- (32) P. Haensel, Astron. Astrophys. 262, 131 (1992).
- (33) O. Benhar and M. Valli, Phys. Rev. Lett. 99, 232501 (2007) [arXiv:0707.2681 [nucl-th]].
- (34) N. Andersson, K. D. Kokkotas and N. Stergioulas, Astrophys. J. 516, 307 (1999) [arXiv:astro-ph/9806089].
- (35) M. Nayyar and B. J. Owen, Phys. Rev. D 73, 084001 (2006) [arXiv:astro-ph/0512041].
- (36) N. Andersson, D. I. Jones and K. D. Kokkotas, Mon. Not. Roy. Astron. Soc. 337, 1224 (2002) [arXiv:astro-ph/0111582].
- (37) J. Abadie et al. [LIGO Scientific Collaboration], arXiv:1006.2535 [Unknown].
- (38) B. J. Owen, arXiv:1006.1994 [gr-qc].
- (39) H. Heiselberg and C. J. Pethick, Phys. Rev. D 48, 2916 (1993).