Spin Evolution of Accreting Neutron Stars: Nonlinear Development of the R-mode Instability

TL;DR

This paper models the nonlinear development of r-mode instability in accreting neutron stars, revealing how it influences their spin evolution, potential thermal runaway scenarios, and implications for gravitational wave detection.

Contribution

It introduces a coupled mode and evolution model that simplifies to key equations, providing new insights into the saturation and thermal behavior of r-mode instability in neutron stars.

Findings

Mode amplitudes settle into quasi-stationary states

Spin frequency limited to about 800 Hz by boundary layer viscosity

Gravitational waves from LMXBs could be detectable by advanced LIGO

Abstract

The nonlinear saturation of the r-mode instability and its effects on the spin evolution of Low Mass X-ray Binaries (LMXBs) are modeled using the triplet of modes at the lowest parametric instability threshold. We solve numerically the coupled equations for the three mode amplitudes in conjunction with the spin and temperature evolution equations. We observe that very quickly the mode amplitudes settle into quasi-stationary states. Once these states are reached, the mode amplitudes can be found algebraically and the system of equations is reduced from eight to two equations: spin and temperature evolution. Eventually, the system may reach thermal equilibrium and either (1) undergo a cyclic evolution with a frequency change of at most 10%, (2) evolve toward a full equilibrium state in which the accretion torque balances the gravitational radiation emission, or (3) enter a thermogravitational runaway on a very long timescale of about $10^6$ years. Alternatively, a faster thermal runaway (timescale of about 100 years) may occur. The sources of damping considered are shear viscosity, hyperon bulk viscosity and boundary layer viscosity. We vary proprieties of the star such as the hyperon superfluid transition temperature T_c, the fraction of the star that is above the threshold for direct URCA reactions, and slippage factor, and map the different scenarios we obtain to ranges of these parameters. For all our bound evolutions the r-mode amplitude remains small $\sim 10^{-5}$. The spin frequency is limited by boundary layer viscosity to $\nu_{max} \sim 800 Hz [S_{ns}/(M_{1.4} R_6)]^{4/11} T_8^{-2/11}$. We find that for $\nu > 700$ Hz the r-mode instability would be active for about 1 in 1000 LMXBs and that only the gravitational waves from LMXBs in the local group of galaxies could be detected by advanced LIGO interferometers.

Figures

Equations

$$\frac{d{C}_{\alpha }}{d\tau }$$

$$\frac{d{C}_{\beta }}{d\tau }$$

$$\frac{d{C}_{\gamma }}{d\tau }$$

$$|C_{\alpha }{|}^2=\frac{{\gamma }_{\beta }{\gamma }_{\gamma }}{4{\tilde{\kappa }}^2{\tilde{\omega }}_{\beta }{\tilde{\omega }}_{\gamma }\Omega }\left[1+{\Omega }^2{\left(\frac{\delta \tilde{\omega }}{{\gamma }_{\beta }+{\gamma }_{\gamma }}\right)}^2\right],$$

$$|C_{\alpha }{|}_0$$

$$|C_{\gamma }{|}_0$$

$$\frac{d{\bar{C}}_{\alpha }}{d\tilde{\tau }}$$

$$\frac{d{\bar{C}}_{\beta }}{d\tilde{\tau }}$$

$$\frac{d{\bar{C}}_{\gamma }}{d\tilde{\tau }}$$

$$\frac{d|{\bar{C}}_{\alpha }|}{d\tilde{\tau }}$$

$$\frac{d|{\bar{C}}_{\beta }|}{d\tilde{\tau }}$$

$$\frac{d|{\bar{C}}_{\gamma }|}{d\tilde{\tau }}$$

$$\frac{d\varphi }{d\tilde{\tau }}$$

$$|{\bar{C}}_{\alpha }{|}^2=\frac{4{\tilde{\gamma }}_{\beta }{\tilde{\gamma }}_{\gamma }}{\tilde{\Omega }|\delta \tilde{\omega }{|}^2}\left(1+\frac{1}{{\tan }^2\varphi }\right),$$

$$|{\bar{C}}_{\beta }{|}^2=\frac{4{\tilde{\gamma }}_{\alpha }{\tilde{\gamma }}_{\gamma }}{\tilde{\Omega }|\delta \tilde{\omega }{|}^2}\left(1+\frac{1}{{\tan }^2\varphi }\right),$$

$$|{\bar{C}}_{\gamma }{|}^2=\frac{4{\tilde{\gamma }}_{\alpha }{\tilde{\gamma }}_{\beta }}{\tilde{\Omega }|\delta \tilde{\omega }{|}^2}\left(1+\frac{1}{{\tan }^2\varphi }\right),$$

$$\tan \varphi =\frac{{\tilde{\gamma }}_{\beta }+{\tilde{\gamma }}_{\gamma }-{\tilde{\gamma }}_{\alpha }}{\tilde{\Omega }|\delta \tilde{\omega }|}.$$

$$J=I\Omega +J_{\mathrm{phys}}.$$

$$\Omega J_{\mathrm{phys}}$$

$$J_{\mathrm{phys}}=M{R}^2(k_{\alpha \alpha }|C_{\alpha }{|}^2+k_{\beta \beta }|C_{\beta }{|}^2+k_{\gamma \gamma }|C_{\gamma }{|}^2),$$

$$k_{\alpha \alpha }=\frac{1}{M{R}^2}\int d^{3}x\rho [(\hat{\Omega }\times {\xi }_{\alpha }^{\star })\cdot (\hat{\Omega }\times {\xi }_{\alpha })-i{\tilde{\omega }}_{\alpha }{\xi }_{\alpha }^{\star }\cdot (\hat{\Omega }\times {\xi }_{\alpha })]$$

$$J_{\mathrm{phys}}=\frac{M{R}^2{\Omega }_c|\delta \tilde{\omega }{|}^2}{(4\tilde{k}{)}^2{\tilde{\omega }}_{\alpha }{\tilde{\omega }}_{\beta }{\tilde{\omega }}_{\gamma }}(k_{\alpha \alpha }|{\bar{C}}_{\alpha }{|}^2{\tilde{\omega }}_{\alpha }$$

$$+k_{\beta \beta }|{\bar{C}}_{\beta }{|}^2{\tilde{\omega }}_{\beta }+k_{\gamma \gamma }{\tilde{\omega }}_{\gamma }|{\bar{C}}_{\gamma }{|}^2).$$

$$\frac{dJ}{dt}=2{\gamma }_{GR}{J}_{c\mathrm{rmode}}+\dot{M}\sqrt{GMR},$$

$$\frac{dJ}{d\tilde{\tau }}=-\frac{6{\tilde{\gamma }}_{\mathrm{GR}}}{\tilde{\Omega }}\frac{M{R}^2{\Omega }_c|\delta \tilde{\omega }|}{(4\tilde{k}{)}^2{\tilde{\omega }}_{\beta }{\tilde{\omega }}_{\gamma }}|{\bar{C}}_{\alpha }{|}^2+\frac{\dot{M}\sqrt{GMR}}{{\Omega }_c\tilde{\Omega }|\delta \tilde{\omega }|}.$$

$$C(T)\frac{dT}{dt}$$

$$L_{\nu }$$

$$R_{\mathrm{dU}}(T/T_p)$$

$$R_{\mathrm{mU}}(T/T_p)$$

$$v=\sqrt{1-\frac{T}{T_p}}\left(1.456-0.157\sqrt{\frac{T_p}{T}}+1.764\frac{T_p}{T}\right).$$

$$C(T)\frac{dT}{d\tilde{\tau }}=\frac{2M{R}^2{\Omega }_c^2|\delta \tilde{\omega }|}{(4\tilde{\kappa }{)}^2{\tilde{\omega }}_{\alpha }{\tilde{\omega }}_{\beta }{\tilde{\omega }}_{\gamma }}({\tilde{\omega }}_{\alpha }{\tilde{\gamma }}_{\alpha v}|{\bar{C}}_{\alpha }{|}^2+{\tilde{\omega }}_{\beta }{\tilde{\gamma }}_{\beta }|{\bar{C}}_{\beta }{|}^2$$

$$+{\tilde{\omega }}_{\gamma }{\tilde{\gamma }}_{\gamma }|{\bar{C}}_{\gamma }{|}^2)+\frac{K_n\dot{M}{c}^2-L_{\nu }(T)}{{\Omega }_c\tilde{\Omega }|\delta \tilde{\omega }|}.$$

$$\frac{d\tilde{\Omega }}{d\tilde{\tau }}=-\frac{6{\tilde{\gamma }}_{\mathrm{GR}}}{{\tilde{\Omega }}^2|\delta \tilde{\omega }|}\frac{{\tilde{\gamma }}_{\beta }{\tilde{\gamma }}_{\gamma }}{4{\tilde{k}}^2\tilde{I}{\tilde{\omega }}_{\beta }{\tilde{\omega }}_{\gamma }}{k}_{\alpha \alpha }\left(1+\frac{1}{{\tan }^2\varphi }\right)$$

$$+\frac{\dot{M}}{{\Omega }_c^2}\frac{\sqrt{GMR}}{M{R}^2\tilde{I}\tilde{\Omega }|\delta \tilde{\omega }|},$$

$$C(T)\frac{dT}{d\tilde{\tau }}=\frac{2M{R}^2{\Omega }_c^2}{(4\tilde{\kappa }{)}^2{\tilde{\omega }}_{\alpha }{\tilde{\omega }}_{\beta }{\tilde{\omega }}_{\gamma }}\frac{{\tilde{\gamma }}_{\alpha }{\tilde{\gamma }}_{\beta }{\tilde{\gamma }}_{\gamma }}{\tilde{\Omega }|\delta \tilde{\omega }|}(\frac{{\tilde{\omega }}_{\alpha }{\tilde{\gamma }}_{\alpha ,v}}{{\tilde{\gamma }}_{\alpha }}+{\tilde{\omega }}_{\beta }$$

$$+{\tilde{\omega }}_{\gamma })\left(1+\frac{1}{{\tan }^2\varphi }\right)+\frac{K_n\dot{M}{c}^2-L_{\nu }(T)}{{\Omega }_c\tilde{\Omega }|\delta \tilde{\omega }|}.$$

$${\gamma }_{jv}(\Omega ,T)={\gamma }_{j\mathrm{sh}}(T)+{\gamma }_{j\mathrm{bl}}(\Omega ,T)+{\gamma }_{j\mathrm{hb}}(\Omega ,T).$$

$${\gamma }_{GR}(\Omega )$$

Full text

I Introduction

R-modes are oscillations in rotating fluids that are due to the Coriolis effect. They are subject to the classical Chandrashekar-Friedman-Shutz (CFS) instability C ; FS , which is driven by the gravitational radiation backreaction force. Andersson nils and Friedman and Morsink Sharon showed that, in the absence of fluid dissipation, r-modes are linearly unstable at all rotation rates. However, in real stars there is a competition between internal viscous dissipation and gravitational driving LOM that depends on the angular velocity $\Omega$ and temperature $T$ of the star. Above a critical curve in the $\Omega-T$ plane the $n=3,m=2$ mode, referred to as ’the r-mode’ in this work, becomes unstable. At first, an unstable r-mode grows exponentially, but soon it may enter a regime where other inertial modes that couple to the r-mode become excited and nonlinear effects become important. Roughly speaking, nonlinear effects first become significant as the amplitude passes its first parametric instability threshold, which is very low ($\sim 10^{-5}$). Modeling and understanding the nonlinear effects is crucial in determining (1) the final saturation amplitude of the r-mode and (2) the limiting spin frequency that neutron stars can achieve. The r-mode amplitude and the duration of the instability are among the main factors that determine whether the associated gravitational radiation could be detectable by laser interferometers on Earth.

The r-mode instability has been proposed as an explanation for the sub-breakup spin rates of both Low Mass X-ray Binaries (LMXBs) Lars ; NKS and young, hot neutron stars LOM ; AKS . The idea that gravitational radiation could balance accretion was proposed independently by Bildsten Lars and Andersson et al. NKS . Cook, Shapiro and Teukolsky saul1 ; saul2 model the recycling of pulsars to millisecond periods via accretion from a Keplerian disk onto a bare neutron star with $M=1.4M_{\odot}$ when $\Omega=0$. Depending on the equation of state they found that spin frequencies of between $\approx 670$ Hz and 1600 Hz could be achieved before mass shedding or radial instability set in (these calculations predated the realization that the r-mode instability could limit the spin frequency). For comparison, the highest observed spin rate of millisecond pulsars is 716 Hz for PSR J1748-2446ad fsharp1 ; fsharp2 . PSR B1937+21, which was discovered in 1982, was the previous fastest known radio pulsar with a spin rate of $642\;\rm{Hz}$ dsharp ; that this “speed” record stood for 24 years suggests that neutron stars rotating this fast are rare. Moreover, based on a Bayesian statistical analysis of the spin frequencies of the 11 nuclear-powered millisecond pulsars whose spin periods are known from burst oscillations, Chakrabarty et al. deepto1 claimed a cutoff limit of $\nu_{\rm{max}}=760\;\rm{Hz}$ (95% confidence); A more recent analysis, which added two more pulsars to the sample, found $\nu_{\rm{max}}=730\;\rm{Hz}$ deepto2 .

At first sight, one might conclude that mass shedding or radial instability sets $\nu_{\rm{max}}$, and that it is just above the record $\nu=716$ Hz determined for PSR J1748-2446ad. However, the nuclear equations of state consistent with this picture all have rather large radii $\approx 16-17$ km for non-rotating 1.4 $M_{\odot}$ models; see Table 1 in Cook et al. saul1 . For these equations of state, the r-mode instability should lead to $\nu_{\rm{max}}$ somewhat below 716 Hz; see Eq. (33) in Sec. V below. Thus, the r-mode instability may prevent recycling by accretion from reaching mass shedding or radial instability. In other words, the detection of the 716 Hz rotator is consistent with accretion spin-up mitigated by the r-mode instability only for equations of state for which mass shedding or radial instability would permit even faster rotation. Ultimately, this may be turned into useful constraints on nuclear equations of state. However, at present the uncertainty in the physics of internal dissipation is a significant hindrance in establishing such constraints.

Since a physical model to follow the nonlinear phase of the evolution was initially unavailable, Owen et al. Owen1998 proposed a simple one-mode evolution model in which they assumed that nonlinear hydrodynamics effects saturate the r-mode amplitude at some arbitrarily fixed value. According to their model, once this maximum allowed amplitude is achieved, the r-mode amplitude remains constant and the star spins down at this fixed amplitude (see Eqs. (3.16) and (3.17) in Ref. Owen1998 ). They used this model to study the impact of the r-mode instability on the spin evolution of young hot neutron stars assuming normal matter. In their calculation they include the effects of shear viscosity and n-p-e bulk viscosity. They found that the star would cool to approximately $10^{9}$ K and spin down from a frequency close to the Kepler frequency to about $100\;\rm{Hz}$ in a period of $\sim 1$ yr Owen1998 .

Most subsequent investigations that did not perform direct hydrodynamic simulations used the one-amplitude model of Ref. Owen1998 for studying the r-mode instability. Levin levin used this model to study the limiting effects of the r-mode instability on the spin evolution of LMXBs, assuming an r-mode saturation amplitude of $\sim 1$; he adopted a modified shear viscosity to match the maximum LMXB spin frequency of $330\;\rm{Hz}$ known in 1999. Levin found that the neutron star followed a cyclic evolution in the $\Omega-T$ phase plane. The star spins up for several million years until it crosses the r-mode stability curve, whereupon the r-mode becomes unstable and the star is viscously heated for a fraction of a year until the r-mode reaches its saturation amplitude ($\sim 1$). At this point the spin and r-mode amplitude evolution equations are changed, following the prescription of Ref. Owen1998 to ensure constant amplitude. The star then spins down by emitting gravitational radiation for another fraction of a year until it crosses the r-mode stability curve again and the instability shuts off. The time period during which the r-mode is unstable was found to be about $10^{-6}$ times shorter than the spin-up time, and Levin concluded that it is unlikely that any neutron stars in LMXBs in our galaxy are currently spinning down and emitting gravitational radiation. However, following work by Arras et al. arras showing that nonlinear effects become significant at small r-mode amplitude, Heyl Heyl varied the saturation amplitude, and found that the duration of the spin-down depends sensitively on it. He predicted that the unstable phase could be as much as 30% of the cyclic evolution for an r-mode saturation amplitude of $\alpha\approx 10^{-5}$, and that this would make some of the fastest spinning LMXBs in our galaxy detectable by interferometers on Earth.

Jones Jones and Lindblom and Owen LO pointed out that if the star contains exotic particles such as hyperons (massive nucleons where an up or down quark is replaced with a strange quark), internal processes could lead to a very high coefficient of bulk viscosity in the cores of neutron stars. While this additional high viscosity coefficient could eliminate the instability altogether in newly born neutron stars Jones ; LO ; HLY ; mohit , Nayyar and Owen mohit proposed that it would enhance the probability of detection of gravitational radiation from LMXBs by blocking the thermal runaway.

The cyclic evolution found by Levin levin and generalized by Heyl Heyl arises when shear or boundary layer viscosity dominates the r-mode dissipation. In the evolutionary picture of Nayyar and Owen mohit , the r-mode first becomes unstable at a temperature where shear and boundary layer viscosity dominate, but the resulting thermal runaway halts once hyperon bulk viscosity becomes dominant. The key feature behind the runaway is that shear and boundary layer viscosities both decrease with increasing temperature, so the instability speeds up as the star grows hotter. However, if the bulk viscosity is sufficiently large the star can cross the r-mode stability curve at a point where the viscosity is an increasing function of temperature. Such scenarios were studied by Wagoner Wagoner for hyperon bulk viscosty with low hyperon superfluid transition temperature; similar evolution was found for strange stars by Andersson, Jones and Kokkotas AJK . In this picture, the star evolves near the r-mode stability curve until an equilibrium between accretion spin-up and gravitational radiation spin-down is achieved. The value of the r-mode amplitude remains below the lowest instability threshold found by Brink et al. Jeandrew1 ; Jeandrew2 ; Jeandrew3 for modes with $n<30$, and hence in this regime nonlinear effects may not play a role.

Schenk et al. Schenk developed a formalism to study the nonlinear interaction of the r-mode with other inertial modes. They assumed a small r-mode amplitude and treated the oscillations of the modes with weakly nonlinear perturbation theory via three-mode couplings. This assumption was tested by Arras et al. arras and Brink et al. Jeandrew1 ; Jeandrew2 ; Jeandrew3 . Arras et al. proposed that a turbulent cascade would develop in the strong driving regime. They estimated that r-mode amplitude was small and could have values between $10^{-1}-10^{-4}$. Brink et al. modeled the star as incompressible and calculated the coupling coefficients analytically. They computed the interaction of about 5000 modes via approximatively 1.3 million couplings of the $10^{9}$ possible couplings among the modes with $n\leq 30$. The couplings were restricted to mode triplets with a fractional detuning $\delta\omega/(2\Omega)<0.002$ since near-resonances promote modal excitation at very small amplitudes. Brink et al. showed that the nonlinear evolution saturates at a very small amplitude, generally comparable to the lowest parametric instability threshold that controls the initiation of energy sharing among the sea of inertial modes. However, Brink et al. did not model accretion spin-up or neutrino cooling in their calculation and only included minimal dissipation via shear viscosity.

In this paper we begin a more complete study of the saturation of the r-mode instability including accretion spin up and neutrino cooling. We use a simple model in which we parameterize uncertain properties of the star such as the rate at which it cools via neutrino emission and the rate at which the energy in inertial modes dissipates via boundary layer effects BU and bulk viscosity. In order to exhibit the variety of possible nonlinear behaviors, we explore a range of models with different neutrino cooling and viscous heating coefficients by varying the free parameters of our model. In particular, we vary: (1) the slippage factor $S_{\rm{ns}}$, which regulates the boundary layer viscosity, between 0 and 1 (see for example LU ; YL ; GA for some models of the interaction between the oscillating fluid core and an elastic crust) ; (2) the fraction of the star that is above the density threshold for direct URCA reactions $f_{\rm{dU}}$, which is taken to be between 0 (0% of the star cools via direct URCA) and 1 (100% of the star is subjected to direct URCA reactions), and in general depends on the equation of state used; and (3) the hyperon superfluidity temperature $T_{c}$, which is believed to be between $10^{9}-10^{10}$ K (We use a single, effective $T_{c}$ rather than modelling its spatial variation.) We focus on $T_{c}\gtrsim 5\times 10^{9}$ K for which nonlinear effects are important. For low $T_{c}\lesssim 3\times 10^{9}$ K, Wagoner Wagoner showed that the evolution reaches a steady state at amplitudes below the lowest parametric instability threshold found by Brink et al. Jeandrew3 . It is important to note that all our evolutions start on the part of the r-mode stability curve that decreases with temperature and that the bulk viscosity does not play a role in any of our bound evolutions.

We include three modes: the r-mode at $n=3$ and the two inertial modes at $n=13$ and $n=14$ that become unstable at the lowest parametric instability threshold found by Brink et al. Jeandrew3 . We evolve the coupled equations for the three-mode system numerically in conjunction with the spin and temperature evolution equations. The lowest parametric instability threshold provides a physical cutoff for the r-mode amplitude. In all cases we investigate, the growth of the r-mode is initially halted by energy transfer to the two daughter modes. We observe that the mode amplitudes settle into a series of quasi-stationary states within a period of a few years after the spin frequency of the star has increased above the r-mode stability curve. These quasi-stationary states are algebraic solutions of the three-mode amplitude equations (see Eqs. (6)) and change slowly as the spin and the temperature of the star evolve. Using these solutions for the mode amplitudes, one can reduce the eight evolution equations (six for the real and imaginary parts of the mode amplitudes, which are complex Schenk ; one for the spin, and one for the temperature) to two equations governing the rotational frequency and the temperature of the star. Our work can be regarded as a minimal physical model for modeling amplitude saturation realistically.

The outcome of the evolution is crucially dependent on whether the star can reach a state of thermal equilibrium. This can be predicted by finding the curve where the viscous heating by the three modes balances the neutrino cooling, referred to below as the Heating = Cooling ($H=C$) curve. The $H=C$ curve can be calculated prior to carrying out an evolution using the quasi-stationary solutions for the mode amplitudes. If the spin frequency of the star upon becoming unstable is below the peak of the $H=C$ curve, then the star will reach a state of thermal equilibrium. When such a state is reached we find several possible scenarios. The star can: (1) undergo a cyclic evolution; (2) reach a true equilibrium in which the accretion torque is balanced by the rate of loss of angular momentum via gravitational radiation; or (3) evolve in thermal equilibrium until it reaches the peak of the $H=C$ curve, which occurs on a timescale of about $10^{6}$ yr, and subsequently enter a regime of thermal runaway. On the other hand, if the star cannot find a state of thermal equilibrium, then it enters a regime of thermogravitational runaway within a few hundred years of crossing the r-mode stability curve. When this happens, the r-mode amplitude increases beyond the second parametric instability, and more inertial modes would need to be included to correctly model the nonlinear effects. This will be done in a later paper.

This paper focuses on showing how nonlinear mode couplings affect the evolution of the temperature and spin frequency of a neutron star once it becomes prone to the r-mode CFS instability. We do this in the context of three mode coupling, which may be sufficient for large enough dissipation. To illustrate the types of behavior that arise, we adopt a very specific model in which the mode frequencies and couplings are computed for an incompressible star, modes damp via shear viscosity, boundary layer viscosity and hyperon bulk viscosity, and the star cools via a mixture of fast and slow processes. This model involves several parameters that are uncertain, and we vary these to find ‘phase diagrams’ in which different generic types of behavior are expected. Moreover, the model itself is simplified: (1) A more realistic treatment of the modes could include buoyant forces, and also mixtures of superfluids or of superfluid and normal fluid in different regions. (2) Dissipation rates, particularly from bulk viscosity, depend on the composition of high density nuclear matter, which could differ from what we assume.

Nevertheless, although the quantitative details may differ from what we compute, we believe that many features of our calculations ought to be robust. More sophisticated treatment of the modes of the star will still find a dense set of modes confined to a relatively small range of frequencies. Most importantly, this set will exhibit numerous three mode resonances, which is the prerequisite for strong nonlinear effects at small mode amplitudes. Thus, whenever the unstable r-mode can pass its lowest parametric instability threshold, it must start exciting its daughters. Whether or not that occurs depends on the temperature dependence of the dissipation rate of the r-mode; for the models considered here, where bulk viscosity is relatively unimportant, soon after the star becomes unstable its r-mode amplitude passes its first parametric instability threshold. Once that happens, the generic types of behavior we find - cycles, steady states, slow and fast runaway - ought to follow suit. The details of when different behaviors arise will depend on the precise features of the stellar model, but the principles we outline here (parametric instability, quasisteady evolution, competition between heating and cooling) ought to apply quite generally.

In Sec. II we describe the evolution equations of the three modes, the angular frequency and the temperature of the neutron star. We first show how the equations of motion for the modes of Schenk *et al. *couple to the rotational frequency of the star in the limit of slow rotation. We then give a short review of the parametric instability threshold and the quasi-stationary solutions of the three-mode system. The thermal and spin evolution of the star is discussed next. This is followed by a description of the driving and damping rates used. Sec. III provides an overview of the results, which includes a discussion of each evolution scenario and of the initial conditions and input physics that lead to each scenario. Sec. IV.1 discusses cyclic evolution in more detail. An evolution that leads to an equilibrium steady state is presented next in Sec. IV.2. The two types of thermal runaway are then discussed in Sec. IV.3. The prospects for detecting gravitational radiation for the evolutions in which the three-mode system correctly models the nonlinear effects are considered in Sec. V. We summarize the results in the conclusion. Appendix A sketches a derivation of the equations of motion for the three modes and Appendix B contains a stability analysis of the evolution equations around the thermal equilibrium state.

II Evolution Equations

II.1 Three mode system: coupling to uniform rotation

In this section we review the equations of motion for the three-mode system in the limit of slow rotation. In terms of rotational phase $\tau$ for the time variable with $d\tau=\Omega\;dt$ Eq. (2.49) of Schenk et al. Schenk can be rewritten as

[TABLE]

Here the scaled frequency $\tilde{\omega}_{j}$ is defined to be $\tilde{\omega}_{j}=\omega_{j}/\Omega$, the dissipation rates of the daughter modes are $\gamma_{\beta}$ and $\gamma_{\gamma}$, $\gamma_{\alpha}$ is the sum of the driving and damping rates of the r-mode $\gamma_{\alpha}=\gamma_{GR}-\gamma_{\alpha\,v}$, and the dimensionless coupling is $\tilde{\kappa}=\kappa/(MR^{2}\Omega^{2})$. These amplitude variables are complex and can be written in terms of the variables of Ref.  Schenk as $C_{j}(t)=\sqrt{\Omega(t)}c_{j}(t)$ (see Appendix A for a derivation of Eqs. 1). The index $j$ loops over the three modes $j=\alpha,\beta,\gamma$, where $\alpha$ labels the r-mode or parent mode and $\beta$ and $\gamma$ label the two daughter modes in the mode triplet.

When the daughter mode amplitudes are much smaller than that of the parent mode, one can approximate the parent mode amplitude as constant. Under this assumption one performs a linear stability analysis on Eqs. (1) and finds the r-mode amplitude when the two daughter modes become unstable (see Eqs. (B5-B7) of Ref. Jeandrew3 for a full derivation). This amplitude is the parametric instability threshold

[TABLE]

where the fractional detuning is $\delta\tilde{\omega}=\tilde{\omega}_{\alpha}-\tilde{\omega}_{\beta}-\tilde{\omega}_{\gamma}$. Thorough explorations of the phase space of damped three-mode systems were performed by Dimant Dimant and Wersinger et al. Wer .

For the three modes at the lowest parametric instability threshold, $\tilde{\omega}_{\alpha}\approx 0.66$, $\tilde{\omega}_{\beta}\approx 0.44$, $\tilde{\omega}_{\gamma}\approx 0.22$, $\tilde{\kappa}\approx 0.19$ and $|\delta\tilde{\omega}|\approx 3.82\times 10^{-6}$. Note that $\tilde{\omega}$ is twice the $w$ of Brink et al. Jeandrew1 ; Jeandrew2 ; Jeandrew3 . Here $\beta$ labels the mode with $n=13,m=-3$ and $\gamma$ labels the $n=14,m=1$ mode. The amplitude the r-mode has to reach before exciting these two daughter modes is $|C_{\alpha}|\approx 1.5\times 10^{-5}\sqrt{\Omega}$ Jeandrew3 .

We next rescale the rotational phase $\tau$ by the fractional detuning as $\tilde{\tau}=\tau|\delta\tilde{\omega}|$ and the mode amplitudes by

[TABLE]

which for the r-mode is, up to a factor of $\sqrt{\tilde{\Omega}}=\sqrt{\Omega/\Omega_{c}}$, the no-damping limit of the parametric instability threshold below which no oscillations will occur. The coupled equations become

[TABLE]

with $\bar{C}_{j}=C_{j}/|C_{j}|_{0}$ and $\tilde{\gamma}_{j}=\gamma_{j}/\Omega_{c}$ being the newly rescaled amplitudes and dissipation/driving rates, respectively.

II.1.1 Quasi-Stationary Solution

In terms of amplitudes and phase variables $C_{j}=|C_{j}|e^{i\phi_{j}}$ Eqs. (4) can be rewritten as

[TABLE]

where we have defined the relative phase difference as $\phi=\phi_{\alpha}-\phi_{\beta}-\phi_{\gamma}$. These equations have the stationary solution

[TABLE]

Note that in the limit in which $\gamma_{\beta}+\gamma_{\gamma}>>\gamma_{\alpha}$ the stationary solution for the r-mode amplitude $|C_{\alpha}|$ is the same as the parametric instability threshold.

II.2 Temperature and Spin Evolution

The spin evolution equation is obtained from conservation of total angular momentum $J$, where

[TABLE]

Following Eq (K39-K42) of Schenk et al. Schenk the physical angular momentum of the perturbation can be written as

[TABLE]

Since the eigenvectors $\xi_{A}\propto e^{im_{A}\phi}$ the cross-terms will vanish for modes with different magnetic quantum numbers $m$ as $\int e^{i(m_{A}-m_{B})\phi}d\phi=0$ for $m_{A}\neq m_{B}$. Eq. (II.2) can be re-written for our triplet of modes as

[TABLE]

where $k_{\alpha\alpha}$ is defined as

[TABLE]

and similarly for $k_{\beta\beta}$ and $k_{\gamma\gamma}$. In terms of the scaled variables $\bar{C}_{j}=C_{j}/|C_{j}|_{0}$ (with $|C_{j}|_{0}$ defined in Eq. (3)) the angular momentum of the perturbation can be written as

[TABLE]

We chose the same normalization for the eigenfuctions as Refs. Schenk ; arras ; Jeandrew1 ; Jeandrew2 ; Jeandrew3 so that at unit amplitude all modes have the same energy $\epsilon_{\alpha}=MR^{2}\Omega^{2}$. The energy of a mode $\alpha$ is $E_{\alpha}=MR^{2}\Omega^{2}|c_{\alpha}|^{2}=MR^{2}\Omega|C_{\alpha}|^{2}$. The rotating frame energy is the same as the canonical energy and physical energy Schenk . The canonical angular momentum and the canonical energy of the perturbation satisfy the general relation $E_{c}=-(\omega/m)J_{c}$ FS .

Angular momentum is gained because of accretion and lost via gravitational waves emission

[TABLE]

where $J_{c\;\rm{rmode}}=-(m_{\alpha}/\omega_{\alpha})\epsilon_{\alpha}|c_{\alpha}|^{2}=-3MR^{2}\Omega|c_{\alpha}|^{2}=-3MR^{2}|C_{\alpha}|^{2}$. Eq. (12) can be rewritten in terms of the scaled variables $\bar{C}_{j}$ as

[TABLE]

Thermal energy conservation gives the temperature evolution equation

[TABLE]

The three terms on the right hand side of the equation represent viscous heating, nuclear heating and neutrino cooling. The specific heat is taken to be $C(T)\approx 1.5\times 10^{38}\;T_{8}\;\rm{erg\;K}^{-1}$, where $T=T_{8}\times 10^{8}$ K. Nuclear heating occurs because of pycnonuclear reactions and neutron emission in the inner crust Brown . At large accretion rates such as that of the brightest LMXBs of $\dot{M}\approx 10^{-8}M_{\odot}/\rm{yr}$, the accreted helium and hydrogen burns stably and most of the heat released in the crust is conducted into the core of the neutron star, where neutrino emission is assumed to regulate the temperature of the star Brown ; Schatz . The nuclear heating constant is taken to be $K_{n}\approx 1\times 10^{-3}$ Brown . Following Ref. Wagoner , we take the neutrino luminosity to be

[TABLE]

where the constants for the modified and direct URCA reactions are defined by $L_{\rm{mU}}=1.0\times 10^{32}\;\rm{erg\;sec}^{-1}$, $L_{\rm{dU}}=f_{\rm{dU}}\times 10^{8}L_{\rm{mU}}$ Yakovlev ; Yakovlev2 , and the electron-ion, neutron-neutron neutrino bremsstrahlung and Cooper pairing of neutrons are given by $L_{\rm{e-i}}=9.1\times 10^{29}\;\rm{erg\;sec}^{-1}$ Brown , $L_{\rm{n-n}}\approx 0.01L_{\rm{mU}}$, $L_{\rm{Cp}}=8.9\times 10^{31}\;\rm{erg\;sec}^{-1}$ Yakovlev3 . The fraction of the star $f_{\rm{dU}}$ that is above the density threshold for direct URCA reactions is in general dependent on the equation of state YP and in this work we treat $f_{\rm{dU}}$ a free parameter with values between 0 and 1.

The proton superfluid reduction factors for the modified and direct URCA reactions are taken from Ref. Yakovlev2 (see Eqs. (32) and (51) in Ref. Yakovlev2 ):

[TABLE]

where the dimensionless gap amplitude $v$ for the singlet type superfluidity is given by

[TABLE]

Similar to Ref. Wagoner , we use $T_{p}=5.0\times 10^{9}$ K. In terms of the scaled variables Eq. (II.2) becomes

[TABLE]

II.3 Temperature and Spin Evolution with the Mode Amplitudes in Quasi-Stationary States

Assuming that the amplitudes evolve through a series of spin- and temperature-dependent steady states, i.e., $dC_{i}/d\tilde{\tau}\approx 0$, the spin and thermal evolution equations can be rewritten by taking $J\approx I\Omega$ and using Eqs. (6) in Eq. (13).

[TABLE]

where $\tilde{I}=I/(MR^{2})$. The thermal evolution of the system is given by

[TABLE]

By setting the right hand side of the above equation to zero, one can find the Heating = Cooling ($H=C$) curve. Below, we find that Eqs. (19)-(20) describe the evolution very well throughout the unstable regime. These equations are a minimal physical model for the effects of nonlinear coupling on r-mode evolution.

II.4 Sources of Driving and Dissipation

The damping mechanisms are shear viscosity, boundary layer viscosity and hyperon bulk viscosity; for modes $j=\alpha,\beta,\gamma$ we write

[TABLE]

The r-mode is driven by gravitational radiation and damped by these dissipation mechanisms, while the pair of daughter modes ($n=13,m=-3$ labeled as $\beta$ and $n=14,m=1$ labeled as $\gamma$) is affected only by the viscous damping. Brink et al. Jeandrew1 ; Jeandrew2 ; Jeandrew3 determined that this pair of modes is excited at the lowest parametric instability threshold. Their model uses the Bryan Bryan modes of an incompressible star, which has the advantage that the mode eigenfrequencies (and eigenfunctions) are known analytically. This enables them to find near resonances efficiently. We are using their results, but we include more realistic effects such as bulk viscosity, whose effect vanishes in the incompressible limit ($\Gamma_{1}\to\infty$ in Eq. (29))

For our benchmark calculations, we adopt the neutron star model of Owen et al. Ref. LOM ($n=1$ polytrope, $M=1.4M_{\odot}$, $\Omega_{c}=8.4\times 10^{3}\;\rm{rad}\;\rm{sec}^{-1}$ and $R=12.53$ km) and use their gravitational driving rate and shear viscous damping rate for the r-mode

[TABLE]

where $\tau_{\rm{sh}}=2.56\times 10^{6}\;\rm{sec}$. (In Sec. V we consider approximate scalings with M and R.)

The damping rate due to shear viscosity for the two daughter modes is calculated using the Bryan modes for a star with the same mass and radius

[TABLE]

The geometric contribution $\gamma_{\rm{sh}}/\eta$ of the individual modes increases significantly with the degree $n$ of the mode scaling approximatively like $n^{3}$ for large $n$ (see Eq. (29) of Brink et al. Jeandrew2 for an analytic fit to the shear damping rates computed for the 5,000 modes in their network), and hence the inertial modes with $n=13$ and $n=14$ have shear damping rates about three orders of magnitude larger than that of the r-mode.

The damping due to boundary layer viscosity is calculated using Eq. (4) of Ref. BU ,

[TABLE]

Analogous to Wagoner Wagoner , we allow the slippage factor $S_{\rm{ns}}$ to vary. The slippage factor is defined by Refs. Wagoner ; LU ; KM to be $S_{\rm{ns}}^{2}=(2S_{\rm{n}}^{2}+S_{\rm{s}}^{2})/3$, with $S_{n}$ being the fractional difference in velocity of the normal fluid between the crust and the core LU and $S_{\rm{s}}$ the fractional degree of pinning of the vortices in the crust KM . Note that $\gamma_{\beta\;\rm{bl}}$ and $\gamma_{\gamma\;\rm{bl}}$ are both greater than $2\times\gamma_{\alpha\;\rm{bl}}$ and can easily be comparable to $\gamma_{\rm{GR}}$ in the unstable regime.

The damping rate due to bulk viscosity produced by out-of-equilibrium hyperon reactions for the r-mode is found by fitting the results of Nayyar and Owen mohit . This rate is taken to have a form similar to that taken by Wagoner Wagoner

[TABLE]

and for the daughter modes

[TABLE]

and similarly foe $\gamma_{\gamma\;\rm{hb}}$. The relaxation timescale

[TABLE]

The reduction factor is taken to be the product of two single-particle reduction factors HLY ; mohit

[TABLE]

where $a=1+0.3118y^{2}$, $b=1+2.556y^{2}$ and $y=\sqrt{1.0-T/T_{c}}(1.456-0.157\sqrt{T_{c}/T}+1.764T_{c}/T)$. The constants $t_{1}\approx 10^{-4}\;\rm{sec}$ and $t_{0\alpha}\approx 0.00058\;\rm{sec}$ are found by fitting the results of Ref. mohit . The factor $f_{\rm{hb}}$ allows for physical uncertainties; we take $f_{\rm{hb}}=1$ throughout the body of the paper since $T_{c}$ , which enters $\gamma_{j\;\rm{hb}}$ exponentially, is also uncertain. For the daughter modes, the dissipation energy due to bulk viscosity is calculated using the modes for the incompressible star. In the slow rotation limit, it is given to leading order in $\Gamma_{1}^{-2}$ by

[TABLE]

This approximation was proposed by Cutler and Lindblom CL and adopted by Kokkotas and Stergioulas KS for the r-mode and by Brink et al. Jeandrew2 for the inertial modes. The adiabatic index $\Gamma_{1}$ is regarded as a parameter; we use $\Gamma_{1}\approx 2$. The damping rate is

[TABLE]

where $\epsilon=MR^{2}\Omega^{2}$ is the mode’s energy in the rotating frame at unit amplitude and $j=\beta,\gamma$. Using this procedure, we calculate

[TABLE]

III Summary of Results

Fig. 1(a) shows possible evolutionary trajectories of a neutron star in the angular velocity-temperature $\tilde{\Omega}-T_{8}$ plane, where $T=T_{8}\times 10^{8}$ K is the core temperature, and $\tilde{\Omega}=\Omega/\Omega_{c}=\Omega/\sqrt{\pi G\bar{\rho}}$ with $\bar{\rho}$ the mean density of the neutron star. Fig. 1(b) displays the regions in $f_{\rm{dU}}-S_{\rm{ns}}$ in which the trajectories occur. Here $f_{\rm{dU}}$ represents the fraction of the star that is above the density threshold for direct URCA reactions and $S_{\rm{ns}}$ is the slippage factor that reduces the relative motion between the crust and the core taking into account the elasticity of the crust LU . The stability regions are shown at fixed hyperon superfluidity temperature, $T_{c}=5.0\times 10^{9}$ K. The initial part of the evolution is similar in all scenarios and can be divided into phases.

Phase 0. Spin up below the r-mode stability curve at $T_{8}=T_{8\,\rm{in}}$ such that nuclear heating balances neutrino cooling.

Phase 1. Linear regime. The r-mode amplitude grows exponentially. The phase ends when the r-mode reaches the parametric instability.

Phase 2. The triplet coupling leads to quasi-steady mode amplitudes. The star is secularly heated at approximately constant $\Omega$ because of viscous dissipation in all three modes.

Phase 3. Several trajectories are possible depending on how the previous phase ends.

a. Fast Runaway. The star fails to reach thermal equilibrium when the trajectory passes over the peak of the Heating = Cooling ($H=C$) curve. This leads to rapid runaway. The daughter modes damp eventually as bulk viscosity becomes important, and the r-mode grows exponentially until the trajectory hits the r-mode stability curve again. This scenario ends as predicted by Nayyar and Owen mohit . However, the r-mode passes its second parametric instability threshold soon after it starts growing again. This requires the inclusion of more modes to follow the evolution, which is the subject of future work.

b. The star reaches thermal equilibrium. There are then three possibilities:

(i) Cycle. The star cools and spins down slowly, descending the $H=C$ curve until it crosses the r-mode stability curve again. At this point the instability shuts off. The star cools back to $T_{8\,\rm{in}}$ at constant $\tilde{\Omega}$ and then the cycle repeats itself. At $T_{c}=5.0\times 10^{9}$ K this scenario occurs for values of $S_{\rm{ns}}<0.50$ and large enough values of $f_{\rm{dU}}$. However, if $T_{c}$ is larger, the cycle region in the $f_{\rm{dU}}$-$S_{\rm{ns}}$ phase space increases dramatically (see Fig. 9(a)). Note that our cycles are different from those obtained by Levin levin in that the spin-down phase does not start when the r-mode amplitude saturates (or in our case when it reaches the parametric instability threshold), but rather when the system reaches thermal equilibrium. The r-mode amplitude does not grow significantly above its first parametric instability threshold, remaining close to $\sim 10^{5}$ and so the part of the cycle in which the r-mode is unstable also lasts longer than in Ref. levin . Also, our cycles are narrow. During spin-down the temperature changes by less than 20 % and $\tilde{\Omega}$ changes by less than 10% of the initial value. (See Sec. 2 for a detailed example.)

(ii) Steady State. For small $S_{\rm{ns}}$ and large enough $f_{\rm{dU}}$ ($f_{\rm{dU}}\gtrsim 5\times 10^{-5}$, $S_{\rm{ns}}\lesssim 0.04$; see Fig. 1(b)) the star evolves towards an $\tilde{\Omega}$ equilibrium. The trajectory either ascends or descends the $H=C$ curve (spins up and heats or spins down and cools). The evolution stops when the accretion torque equals the gravitational radiation emission.

(iii) Slow Runaway. For small $S_{\rm{ns}}$ and very small $f_{\rm{dU}}$ ($S_{\rm{ns}}\lesssim 0.03$, $f_{\rm{dU}}<5\times 10^{-5}$) the star ascends the $H=C$ curve until the peak is overcome and subsequently a runaway occurs. The daughter modes eventually damp and the r-mode grows exponentially until it crosses its second parametric instability threshold and more modes need to be included.

Bulk viscosity only affects the runaway evolutions; the cyclic and steady state evolutions found here would be the same if there were no hyperon bulk viscosity. For large $T_{c}\sim 10^{10}$, or for models with no hyperons at all, there would be no runaway region (See Fig. 9(a) for an $f_{\rm{dU}}-S_{\rm{ns}}$ scenario space with a larger $T_{c}=6.5\times 10^{9}$ K where the fast runaway region has shrunk dramatically and the slow runaway region has disappeared.)

IV Possible Evolution Scenarios

In this section we examine examples of the different types of evolution in more detail. We assume $\dot{M}=10^{-8}M_{\odot}$/yr and $T_{c}=5.0\times 10^{9}$ K.

IV.1 Cyclic Evolution

In this sub-section we present the features of typical cyclic trajectories of neutron stars in the angular velocity temperature plane in more detail. We focus on two cases: (C1) $S_{\rm{ns}}=0.10$ and $f_{\rm{dU}}=0.15$ and (C2) $S_{\rm{ns}}=0.35$ and $f_{\rm{dU}}=0.142$. In this scenario the 3-mode system is sufficient to model the nonlinear effects and successfully stops the thermal runaway. The numerical evolution is started once the star reaches the r-mode stability curve. The initial temperature of the star is at the point where nuclear heating equals neutrino cooling in Eq. (18) that is approximately $T_{8\,\rm{in}}\approx 3.29$ for both cases. The initial $\Omega$ is the angular velocity that corresponds to this temperature on the r-mode stability curve, which differs for the different $S_{\rm{ns}}$ ($\tilde{\Omega}_{\rm{in}}=0.183$ for $C_{1}$ and $\tilde{\Omega}_{\rm{in}}=0.288$ for $C_{2}$).

Figs.  2(a) and (b) display the cyclic evolution for trajectories $C_{1}$ and $C_{2}$ of Fig. 1(b). In leg $a_{1}\to b_{1}$ of the trajectory the r-mode and, once the r-mode amplitude increases above the first parametric instability threshold, the two daughter modes it excites, viscously heat up the star until point $b_{1}$ when the neutrino cooling balances the viscous dissipation. This part of the evolution occurs at constant angular velocity over a period of $t_{\rm{heat-up}}\approx 100\;\rm{yr}$ and a total temperature change $(\Delta T)_{a_{1}-b_{1}}\approx 0.80$ ($\approx 24$% of $T_{8\,\rm{in}}$). The points where the viscous heating compensates the neutrino cooling are represented by the Heating = Cooling ($H=C$) curve. This is determined by setting Eq. 18 to zero and using the quasi-stationary solutions given by Eq. (6) for the three modes on the right hand side. The star continues to evolve on the $H=C$ curve for part $b_{1}\to c_{1}$ of the trajectory as it spins down and cools down back to the r-mode stability curve. This spin-down stage lasts a time $t_{\rm{spin-down}\,b_{1}-c_{1}}\approx 23,000\;\rm{yr}$ that is much longer than the heat-up period. This timescale is very sensitive to changes in the slippage factor and can reach $10^{6}$ yr for smaller values of $S_{\rm{ns}}$ that are close to boundary of the steady state region. The cycle is very narrow in angular velocity with a total angular velocity change of less than 4%, $(\Delta\tilde{\Omega})_{b_{1}-c_{1}}\approx 0.0066$. The temperature also changes by only about 2%, $(\Delta T_{8})_{b_{1}-c_{1}}\approx 0.08$ in this spin-down period. Segment $c_{1}\to d_{1}$ represents the cooling of the star to the initial temperature on a timescale of $\sim 2,000$ yr. In part $d_{1}\to a_{1}$ the star spins up by accretion at constant temperature back to the original crossing point on the r-mode stability curve. This last part of the trajectory is the longest-lasting one, taking $\approx 200,000$ yr at our chosen $\dot{M}$ of $10^{-8}M_{\odot}\rm{yr}^{-1}$. The cycle $C_{2}$ in Fig.  2(b) proceeds in a similar fashion. It is important to note that this configuration is close to the border between the “FAST RUNAWAY” and “CYCLE” regions and therefore close to the peak of the $H=C$ curve. Configurations above this peak (e.g., with the same $f_{\rm{dU}}$ and higher $S_{\rm{ns}}$) will go through a fast runaway.

Fig. 3(a) shows the evolution of the three modes in the first few years after the star first reaches the r-mode stability curve. In this region the r-mode is unstable and initially grows exponentially. Once it has increased above the first parametric instability threshold the daughter modes are excited. The oscillations of the three modes display some of the typical dynamics of a driven three-mode system. When the r-mode transfers energy to the daughter modes they increase exponentially while the r-mode decreases. Similarly, when daughter modes decrease the r-mode increases. The viscosity damps the oscillations and the r-mode amplitude settles at a value close to the parametric instability threshold. Fig. 3(b) displays the evolution of the r-mode amplitude divided by the parametric instability threshold on a longer timescale. It can be seen that the r-mode never grows significantly beyond this first threshold. Fig. 3(c) shows the evolution of the parametric instability threshold as a function of time. The threshold increases as the temperature increases and the star is viscously heated by the three modes. When the star spins down in thermal equilibrium, the threshold decreases to a value close to its initial value.

IV.2 Steady State Evolution

This sub-section focuses on evolutions that lead to a steady equilibrium state in which the rate of accretion of angular momentum is balanced by the rate of loss via gravitational radiation emission. This scenario is restricted to stars with small slippage factor ($S_{\rm{ns}}\lesssim 0.04$, see Fig. 1(b)) and boundary layer viscosity. A typical trajectory of a star that reaches such an equilibrium is shown in Fig. 4. As always, we start the evolution at the point on the r-mode stability curve at which the nuclear heating balances neutrino cooling. Above the r-mode stability curve the gravitational driving rate is greater than the viscous damping rate and the r-mode grows exponentially until nonlinear effects become important. In this case, as in the cyclic evolution, the triplet of modes at the lowest parametric instability threshold is sufficient to stop the thermal runaway. The r-mode remains close to the first instability threshold for the length of the evolution and after a few oscillations the three modes settle into their quasi-stationary states, which change only secularly as the spin and temperature of the star evolve. The modes heat the star viscously at constant $\tilde{\Omega}$ in segment $a\to b$ of the trajectory for $t_{\rm{heat-up}}\approx 1,100\;\rm{yr}$. At point $b$, the star reaches a state of thermal balance. In leg $b\to c$ the star continues its evolution in thermal equilibrium and slowly spins up due to accretion until the angular velocity evolution also reaches an equilibrium. The timescale to reach an equilibrium steady state is $t_{\rm{steady}}\approx 3.5\times 10^{6}\;\rm{yr}$ for this set of parameters.

Fig. 5 displays the possible initial values for the angular velocity $\tilde{\Omega}$ and temperature $T_{8}$ of the star that lead to a balancing between the accreted angular momentum and the angular momentum emitted in gravitational waves. The fraction of the star that is above the threshold for direct URCA reactions and the slippage factor are varied within the corresponding “STEADY STATE” region of Fig. 1(b). The final equilibrium values are also displayed and cluster in a narrower region than the initial values. Because viscosity is so small in this regime, the values of $\Omega$ also tend to be small. Thus, although an interesting physical regime, this case is most likely not relevant to recycling by accretion to create pulsars with spin frequencies as large as 716 Hz. Note that a steady state can be achieved when $S_{\rm{ns}}=0$. This is the probable end state of the problem first calculated by Levin levin . The reason we do not find a cycle at low $S_{\rm{ns}}$ is twofold: (1) the shear viscosity we are using is lower (shear viscosity in Ref. levin is amplified by a factor of 244), and (2) the nonlinear couplings keep all mode amplitudes small.

IV.3 Thermal Runaway Evolutions

We now consider evolutions in which the three-mode system is not sufficient to halt the thermal runaway. We observe two such scenarios. In the first scenario, the star is unable to reach thermal equilibrium. The runaway occurs on a period much shorter than the accretion timescale and so the whole evolution is at approximately constant angular frequency. In the second scenario, the star reaches a state of thermal equilibrium but the spin evolution does not reach a steady state. The star continues to spin up by accretion until it climbs to the peak of the $H=C$ curve, thermal equilibrium fails and a runaway occurs.

IV.3.1 Fast Runaway

A typical trajectory of a star that goes through a rapid thermal runaway is displayed in Fig. 6. This star has $S_{\rm{ns}}=0.25$ and $f_{\rm{dU}}=0.058$. Initially, the growth of the r-mode is halted by the two daughter modes once the lowest parametric instability threshold is crossed, and the three modes settle in the ($\Omega$,$T$)-dependent quasi-stationary states of Eqs. (6). They viscously heat up the star until hyperon bulk viscosity becomes important for the daughter modes. As the amplitudes of the daughter modes decrease the coupling is no longer strong enough to drain enough energy to stop the growth of the r-mode. The daughter modes are completely damped and the r-mode increases exponentially. The system goes back to the one-mode evolution described by Ref. mohit .

Fig. 6(a) and (b) compare both the temperature evolution and the trajectory in the $\tilde{\Omega}-T_{8}$ plane of the star for a simulation solving the full set of equations to a simulation that assumes quasi-stationary solutions for the three amplitudes and evolves only the angular velocity and temperature of the star. It can be seen that the steady state approximation is very good until the thermal runaway occurs. Afterward, the temperature evolution of the reduced equations is offset slightly from the quasi-steady result and intersects the r-mode instability curve sooner. This evolution is similar to that described by Nayyar and Owen mohit . However, the r-mode crosses its second lowest parametric instability much earlier in the evolution (see the ‘X’ in the figure), and at that point more modes need to be included to model the instability accurately. Thus, we cannot be sure that a runaway must occur in this case. We shall return to this issue in a subsequent paper.

IV.3.2 Slow Runaway

In this section we examine evolutions in which the neutron star has both a very small slippage factor, $S_{\rm{ns}}\lesssim 0.03$, and only a small percentage of the star is above the threshold for direct URCA reactions, $f_{\rm{dU}}<5\times 10^{-5}$. A trajectory for this kind of evolution is displayed in Fig. 7. After the star crosses the r-mode stability curve, the r-mode increases beyond the first parametric instability threshold, and its growth is temporarily stopped by energy transfer to the daughter modes. As in the previous scenarios, the star is viscously heated by the mode triplet at constant $\Omega$ in part $a\to b$ of the trajectory on a timescale of about $5,000$ yr. At point $b$, it reaches thermal equilibrium. In leg $b\to c$ of the trajectory, the star continues its evolution by ascending the $H=C$ curve and spinning up because of accretion for about $2\times 10^{6}$ yr without finding an equilibrium state for the angular momentum evolution. Once it reaches the peak of the $H=C$ curve, the cooling is no longer sufficient to stop the temperature from increasing exponentially and a thermal runaway occurs. The cross mark ‘X’ on the trajectory shows the point at which the r-mode amplitude crosses its second lowest parametric instability threshold. At this stage more inertial modes need to be included to model the rest of this evolution correctly. As for the cases that evolve to steady states, these long-timescale runaways tend to occur at low spin rates.

V Probability of Detection

Fig. 8 shows how the time the star spends above the r-mode stability curve changes when $S_{\rm{ns}}$ and $f_{\rm{dU}}$ are varied. For large enough values of $S_{\rm{ns}}$ the boundary layer viscosity dominates. In this region of phase space the spin-down timescale can be approximated by

[TABLE]

where $M_{1.4}=M/(1.4M_{\odot})$, $R_{6}=R/(10^{6}$cm), $\nu_{\rm{kHz}}=\nu/1kHz$, $\tilde{I}=0.261$ Owen1998 , the r-mode amplitude at its parametric instability threshold $|c_{\alpha}^{\rm{th}}|\approx|\delta\tilde{\omega}|/(4\tilde{\kappa}\sqrt{\tilde{\omega}_{\beta}\tilde{\omega}_{\gamma}})\approx 1.5\times 10^{-5}$, and $\bar{C}_{\alpha}=\sqrt{\tilde{\Omega}}|c_{\alpha}|/|c_{\alpha}|^{\rm{th}}$. This approximation agrees with spin-up timescales obtained from our simulations to $\sim 25\%$.

The maximum $\nu$ is approximately the same as the initial frequency, and can be determined by equating the driving and damping rate of the r-mode, since it is on the r-mode stability curve

[TABLE]

Thus, the spin-down timescale is very sensitive to the slippage factor $t_{\rm{spin-down}}\propto S_{\rm{ns}}^{-24/11}(\Delta\nu_{\rm{kHz}}/\nu_{\rm{kHz}})$. The dependences on $f_{\rm{dU}}$ and accretion rate $\dot{M}$ are much weaker; a rough approximation, obtained by matching direct URCA cooling and nuclear heating, is $T_{8\,\rm{in}}$ $\propto$ $\dot{M}^{1/6}f_{\rm{dU}}^{-1/6}R_{6}^{-1/6}M_{1.4}^{-1/9}$, and $\nu_{\rm{max}}\propto S_{\rm{ns}}^{4/11}f_{\rm{dU}}^{1/33}\dot{M}^{-1/33}R_{6}^{-1/3}M_{1.4}^{-34/99}$. The gravitational wave amplitude measured at distance d brady ; Owen2001 is

[TABLE]

Taking $\nu\approx\nu_{\rm{max}}$ gives

[TABLE]

The maximum distance at which sources could be detected by Advanced LIGO interferometers, assuming $h_{\rm{min}}=10^{-27}$, brady is

[TABLE]

Eqs. (33) and (V) imply that gravitational radiation from the r-mode instability may only be detectable for sources in the Local Group of galaxies. Eq. (33) implies that for accretion to be able to spin up neutron stars to $\nu\gtrsim 700$ Hz, we must require $(S_{\rm{ns}}/M_{1.4}R_{6}\sqrt{T_{8\rm{in}}})^{4/11}\gtrsim 1$. Assuming this to be true, $d_{\rm{max}}\lesssim 1$-1.5 Mpc. However, $t_{\rm{spin-down}}\approx 1000$ yr at most, making detection unlikely for any given source. Moreover, unless $S_{\rm{ns}}$ can differ substantially from one neutron star to another, only those with $\nu$ given by Eq. (33) can be r-mode unstable. Slower rotators, including almost all LMXBs, are still in their stable spin-up phases.

Still more seriously, Fig. 1(b) shows that spin cycles are only possible for $S_{\rm{ns}}\lesssim 0.50$, assuming $T_{c}\approx 5.0\times 10^{9}$ K; Eq. (33) then implies $\nu\lesssim 450$ Hz. This would restrict detectable gravitational radiation to galactic sources, although the duration of the unstable phase could be longer.

Within the context of our three mode calculation, $S_{\rm{ns}}>0.50$, which is needed for explaining the fastest pulsars, would imply fast runaway. There are two possible resolutions to this problem. One is that including additional modes prevents the runaway; we shall investigate this in subsequent papers. The second is that $T_{c}$ is larger, or that neutron stars do not contain hyperons (e.g., because they are insufficiently dense). Fig. 9(a) shows the same phase plane as Fig. 1(b) but with $T_{c}=6.5\times 10^{9}$ K, and Fig. 9(b) shows the results for $t_{\rm{spin-down}}$ analogous to Fig. 8. Larger $T_{c}$ permits spin cycles for higher values of $S_{\rm{ns}}$ (and hence $\nu$), but the time spent in the unstable regime is shorter.

VI Conclusions

In this paper, we model the nonlinear saturation of unstable r-modes of accreting neutron stars using the triplet of modes formed from the $n=3,m=2$ r-mode and the the first two near resonant modes that become unstable ($n=13,m=-3$ and $n=14,m=1$) by coupling to the r-mode. This is the first treatment of the spin and thermal evolution including the nonlinear saturation of the r-mode instability to provide a physical cutoff by energy transfer to other modes in the system. The model includes neutrino cooling and shear, boundary layer and hyperon bulk viscosity. We allow for some uncertainties in neutron star physics that is not yet understood by varying the superfluid transition temperature, the slippage factor that regulates the boundary layer viscosity, and the fraction of the star that is above the density threshold for direct URCA reactions. In all our evolutions we find that the mode amplitudes quickly settle into a series of quasi-stationary states that can be calculated algebraically, and depend weakly on angular velocity and temperature. The evolution continues along these sequences of quasi-steady states as long as the r-mode is in the unstable regime. The spin and temperature of the neutron star can follow several possible trajectories depending on interior physics. The first part of the evolution is the same for all types of trajectories: the star viscously heats up at constant angular velocity.

If thermal equilibrium is reached, we find several possible scenarios. The star may follow a cyclic evolution, and spin down and cool in thermal equilibrium until the r-mode enters the stable regime. It subsequently cools at constant $\Omega$ until it reaches the initial temperature. At this point the star starts spinning up by accretion until the r-mode becomes unstable again and the cycle is repeated. The time the star spends in the unstable regime is found to vary between a few hundred years (large $S_{\rm{ns}}\sim 1$) and $10^{6}$ yr (small $S_{\rm{ns}}\sim 0.05$). Our cycles are different from those previously found by Ref. levin in that our amplitudes remain small, $\sim 10^{-5}$, which slows the viscous heating and causes the star to spend more time in the regime where the r-mode instability is active. Furthermore, we find that the star stops heating when it reaches thermal equilibrium and not when the r-mode reaches a maximum value. The cycles we find are narrow with the spin frequency of the star changing less than 10% even in the case of high spin rates $\sim 750$ Hz. Other possible trajectories are an evolution toward a full steady state in which the accretion torque balances the gravitational radiation emission, and a very slow thermogravitational runaway on a timescale of $\sim 10^{6}$ yr. These scenarios occur for very low viscosity ($S_{\rm{ns}}\lesssim 0.04$). Although theoretically interesting, they do not allow for very fast rotators of $\sim 700$ Hz.

Alternatively, if the star does not reach thermal equilibrium, we find that it continues heating up at constant spin frequency until it enters a regime in which the r-mode is no longer unstable. This evolution is similar to that predicted by Nayyar and Owen mohit . However, the r-mode grows above its second parametric instability threshold fairly early in its evolution and at this point more inertial modes should be excited and the three-mode model becomes insufficient. Modeling this scenario accurately is subject of future work.

We have focused on cases with $T_{c}\gtrsim 5\times 10^{9}$ K. These are cases for which the nonlinear effects are substantial. In this regime, hyperon bulk viscosity is not important except for thermal runaways where we expect other mode couplings, ignored here, to play important roles. Fast rotation requires large dissipation, as has long been recognized BU ; levin and these models can only achieve $\nu\gtrsim 700$ Hz if boundary layer viscosity is very large. Alternatively, at lower $T_{c}\lesssim 3\times 10^{9}$ K, large rotation rates can be achieved at r-mode amplitudes below the first parametric instability threshold Wagoner . Nayyar and Owen found that increasing the mass of the star for the same equation of state makes the hyperon bulk viscosity become important at lower temperatures mohit . Conceivably, there are accreting neutron stars with relatively low masses that have lower central densities and small hyperon populations. These could evolve as detailed here and only spin up to modest frequencies. Hyperons could be more important in more massive neutron stars leading to larger spin rates and very small steady state r-mode amplitude as found by Wagoner Wagoner .

Our models imply small r-mode amplitudes of $\sim 10^{-5}$ and therefore gravitational radiation detectable by advanced LIGO interferometers only in the local group of galaxies up to a distance of a few Mpc. The r-mode instability puts a fairly stringent limit on the spin frequencies of accreting neutron stars of $\nu_{\rm{max}}\approx 800Hz[S_{\rm{ns}}/(M_{1.4}R_{6})]^{4/11}T_{8}^{-2/11}$. In order to allow for fast rotators of $\gtrsim 700$ Hz in our models a large boundary layer viscosity with $(S_{\rm{ns}}/M_{1.4}R_{6}\sqrt{T_{8\rm{in}}})^{4/11}\sim 1$ is required. Slippage factors of order $\sim 1$ lead to time periods on which the r-mode is unstable with a timescale of at most 1000 yr, which is about $10^{-3}$ times shorter than the accretion timescale. This would mean that only about 1 in 1000 LMXBs in the galaxy are possible LIGO sources. However, lower slippage factors lead to a longer duration of the gravitational wave emission, but also lower frequencies. We also note that in this model we have considered only very fast accretors with $\dot{M}\sim 10^{-8}M_{\odot}\rm{yr}^{-1}$ and most LMXBs in our galaxy accrete at slower rates. Investigations with more accurate nuclear heating models are a subject for future work.

Our analysis could be made more realistic in several ways, such as by including the effects of magnetic fields, compressibility, multi-fluid composition AnderssonSup , superfluidity, superconductivity, etc. These features would render the model more realistic, but its generic features ought to persist, since the upshot would still be a dense set of mode frequencies exhibiting three mode resonances and parametric instabilities with low threshold amplitudes. Although the behavior of the star would differ quantitatively in a model different from ours in detail, we expect the qualitative behaviors we have found to be robust, as they are well described by quasi-stationary mode evolutions whose slow variations are determined by competitions between dissipation and neutrino cooling, and accretion spin-up and gravitational radiation spin-down. In our model, it seems that three mode evolution involving interactions of the r-mode with two daughters at the lowest parametric instability threshold is often sufficient to quench the instability. Our treatment is inadequate to follow what happens when the system runs away; for this, coupling to additional modes is essential. For this regime, a generalization of the work of Brink et al. Jeandrew1 ; Jeandrew2 ; Jeandrew3 that includes accretion spin-up, viscous heating and neutrino cooling would be needed. Such a calculation is formidable even in a “simple” model involving coupled inertial modes of an incompressible star.

Acknowledgments

It is a pleasure to thank Jeandrew Brink and Éanna Flanagan for useful discussions. RB would especially like to thank Jeandrew for useful discussions, encouragement and advice at the beginning of this project, without which the project would not have been started. RB is very grateful to Gregory Daues for steady encouragement and support for the duration of this project, and also to Gabrielle Allen and Ed Seidel. This research was funded by grants NSF AST-0307273, NSF AST-0606710 and NSF PHY-0354631.

APPENDIX A

This appendix will sketch the derivation of Eqs. (1) from the Lagrangian density. We follow closely Appendix A in Schenk et al., which contains the derivation of the equations of motion for constant $\Omega$.

The Lagrangian density as given by Eq. (A1) in Schenk et al. Schenk is

[TABLE]

where the operators ${\bf B}\cdot{\bf\xi}=2{\bf\Omega}\times{\bf\xi}$ and

[TABLE]

with $\phi_{\rm{rot}}=-(1/2)({\bf\Omega}\times{\bf x})^{2}$. We are interested in a situation where the uniform angular velocity of the star changes slowly on the timescale of the rotation period itself. In order to remove the time dependence we define the new displacement and time variables

[TABLE]

In terms of these new variables the Lagrangian density can be written as

[TABLE]

where the primes denote derivatives with respect to $\tau$, ${\bf\tilde{B}}=\Omega^{-1}B$ and ${\bf\tilde{C}}=\Omega^{-2}C$. The momentum canonically conjugate to ${\bf\tilde{\xi}}$ is

[TABLE]

The associated Hamiltonian density is

[TABLE]

Hamilton’s equations of motions can be written as

[TABLE]

where

[TABLE]

the operator $T$ is $T=T_{0}+T_{1}$ with

[TABLE]

and

[TABLE]

and

[TABLE]

We assume solutions of the form $\tilde{\zeta}(\tau,{\bf x})=e^{i\tilde{\omega}t}\tilde{\zeta}({\bf x})$. Specializing to the case of no forcing term $a_{\rm{ext}}=0$ leads to the eigenvalue equation

[TABLE]

Since the operator $T_{0}$ is not Hermitian it will have distinct right and left eigenvectors. Similar to Schenk et al. Schenk we label the right eigenvectors of T as $\tilde{\zeta}_{A}$, and the associated eigenfrequencies as $\tilde{\omega}_{A}=\omega_{A}/\Omega$, and the eigenvalue equation above becomes

[TABLE]

The left eigenvectors $\chi_{A}$ satisfy

[TABLE]

where

[TABLE]

For simplicity, in this appendix we specialize to the case of no Jordan chains when the set of right eigenvectors forms a complete basis. The orthonormality relation between right and left eigenvectors is

[TABLE]

We can expand $\zeta(\tau,{\bf x})$ in this basis as

[TABLE]

where the coefficients $C_{A}$ are given by the inverse of this mode expansion

[TABLE]

Using Eqs. (APPENDIX B,A-9,A-11) in Eq. (A-7) leads to the equations of motion for the mode amplitudes

[TABLE]

where $g(\tau)=(\sqrt{\Omega})^{\prime\prime}/\sqrt{\Omega}$. Following Sec. IV of Schenk et al. Schenk we replace the externally applied acceleration by the nonlinear acceleration given by Eq. (4.2) of Ref. Schenk . The inner product can be written in terms of the displacement variable $\tilde{\xi}$. The left eigenvectors are

[TABLE]

where ${\bf\tilde{\tau}}_{A}$ can be chosen to be proportional to ${\bf\tilde{\xi}}_{A}$ because they satisfy the same matrix equation.

[TABLE]

which corresponds to Eq. (A-45) in Schenk et al. Schenk with the proportionality constant $\tilde{b}_{A}=\Omega^{-1}b_{A}=MR^{2}/\tilde{\omega}_{A}$ (also given by Eq. (2.36) of Ref. Schenk ).

The equations of motion for the mode amplitudes become

[TABLE]

where the nonlinear coupling $\tilde{\kappa}_{ABC}=\kappa_{ABC}/(MR^{2}\Omega^{2})$ and $\kappa_{ABC}$ is explicitly give by Eq. (4.20) of Ref. Schenk . The $g(\tau)$ integral mixes only modes with $m_{A}=m_{B}$ because of the $e^{im\phi}$ dependence of the displacement eigenvectors $\tilde{\xi}$. ($\int_{0}^{2\pi}\;d\phi e^{i(m_{A}-m_{B})\phi}=0$ if $m_{A}\neq m_{B}$.) So, this term will be zero for our mode triplet. Also, in the case of a single mode triplet there is only one coupling and Eqs. (APPENDIX A) take the form of Eqs. (1).

APPENDIX B

In this appendix we study the behavior of the mode amplitudes and temperature near equilibrium assuming constant angular velocity. We are performing a first order expansion of Eqs. (5) and (18). Similar to Ref. WHL , each of the five variables is expanded about its equilibrium $(X_{j})_{e}$ as follows

[TABLE]

where the perturbation $|\zeta_{j}|<<1$ and $j=\alpha,\beta,\gamma,T$. The expansion leads to a first order differential equation for each $\zeta_{j}$

[TABLE]

where the equilibrium amplitudes $|C_{j}|_{e}$ have been written in terms of the corresponding driving and damping rates using Eqs. (6). Eq. (APPENDIX B) can be written in matrix form as

[TABLE]

Let $\zeta_{j}\propto\exp(\lambda\tilde{\tau})$. The determinant $||A^{ij}-\lambda\delta^{ij}||=0$ leads to the eigenvalue equation

[TABLE]

The coefficients $a_{j}$ with $j=0,4$ are

[TABLE]

The eigenvalues can be approximated as

[TABLE]

where $\epsilon=(a_{2}-a_{3}a_{4})/a_{4}$ and $w=\sqrt{a_{1}/a_{3}}$. The system is unstable when $a_{2}-a_{3}a_{4}>0$ or $a_{0}<0$. The first two eigenvalues will have a negative real part as long as $\tilde{\gamma}_{\beta}+\tilde{\gamma}_{\gamma}>\tilde{\gamma}_{\alpha}$. If the heating compensates the cooling of the star $a_{0}\approx 0$ and becomes negative if the star can not reach thermal equilibrium. The other critical stability condition $a_{2}-a_{3}a_{4}=0$ can be written as

[TABLE]

where $\Gamma_{\beta}=\gamma_{\beta}/\gamma_{\alpha}$ and $\Gamma_{\gamma}=\gamma_{\gamma}/\gamma_{\alpha}$. Note that we have ignored the smaller terms of order $O([\tilde{\gamma}_{\alpha}/(\tilde{\Omega}|\delta\tilde{\omega}|)]^{5})$. This condition can be rewritten by defining variables $D_{1}=\Gamma_{\beta}+\Gamma_{\gamma}$ and $D_{2}=\Gamma_{\beta}-\Gamma_{\gamma}$

[TABLE]

If $D_{2}=0$ then the equation has one solution $D_{1}=1+\sqrt{3}$ for $D_{1}>2$, which corresponds to $\Gamma=\Gamma_{\beta}=\Gamma_{\gamma}=1.37$ and matches the result of Wersinger et al. Wer . For the viscosity we consider (see Sec. II.4) $a_{2}-a_{3}a_{4}<0$.