Spin Evolution of Accreting Neutron Stars: Nonlinear Development of the R-mode
Instability
Ruxandra Bondarescu, Saul A. Teukolsky, and Ira Wasserman
Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853
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 that change slowly as the temperature and spin
of the star evolve. 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. The
evolution of the neutron star angular velocity and temperature follow easily calculated trajectories
along these sequences of quasi-stationary states. The outcome depends on whether or not the
star will reach thermal equilibrium, where the viscous heating by the three modes is equal to
the neutrino cooling (H = C curve). If, when the r-mode becomes unstable, the star spins at a
frequency below the maximum of the H = C curve, then it will reach a state of thermal equilibrium.
It can then either (1) undergo a cyclic evolution with a small cycle size 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 ≈ 106 years. If the star does not reach a state of thermal equilibrium, then a faster
thermal runaway (timescale of ≈ 100 years) occurs and the r-mode amplitude increases above the
second parametric instability threshold. Following this evolution requires more inertial modes to
be included. The sources of damping considered are shear viscosity, hyperon bulk viscosity and
viscosity within the core-crust boundary layer. We vary proprieties of the star such as the hyperon
superfluid transition temperature Tc, 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. We focus on Tc & 5 × 10
9 K where nonlinear effects are important. Wagoner [1] has
shown that a very low r-mode amplitude arises at smaller Tc. For all our bounded evolutions the
r-mode amplitude remains small ∼ 10−5. The spin frequency of accreting neutron stars is limited by
boundary layer viscosity to νmax ≈ 800Hz[Sns/(M1.4R6)]
4/11T
−2/11
. Fast rotators are allowed for
[Sns/(M1.4R6)]
4/11T
−2/11
∼ 1 and we find that in this case 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.
PACS numbers: 04.40.Dg, 04.30.Db, 97.10.Sj, 97.60.Jd
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 [2, 3],
which is driven by the gravitational radiation backreac-
tion force. Andersson [4] and Friedman and Morsink [5]
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 [6] that depends on
the angular velocity Ω and temperature T of the star.
Above a critical curve in the Ω−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 exponen-
tially, but soon it may enter a regime where other in-
ertial modes that couple to the r-mode become excited
and nonlinear effects become important. Roughly speak-
ing, nonlinear effects first become significant as the am-
plitude passes its first parametric instability threshold,
which is very low (∼ 10−5). Modeling and understand-
ing the nonlinear effects is crucial in determining (1) the
final saturation amplitude of the r-mode and (2) the lim-
iting 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 as-
sociated gravitational radiation could be detectable by
laser interferometers on Earth.
The r-mode instability has been proposed as an expla-
nation for the sub-breakup spin rates of both Low Mass
X-ray Binaries (LMXBs) [7, 8] and young, hot neutron
stars [6, 9]. The idea that gravitational radiation could
balance accretion was proposed independently by Bild-
sten [7] and Andersson et al. [8]. Cook, Shapiro and
Teukolsky [10, 11] model the recycling of pulsars to mil-
lisecond periods via accretion from a Keplerian disk onto
a bare neutron star with M = 1.4M⊙ when Ω = 0. De-
pending on the equation of state they found that spin
frequencies of between ≈ 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 millisec-
ond pulsars is 716 Hz for PSR J1748-2446ad [12, 13].
http://arxiv.org/abs/0704.0799v3
PSR B1937+21, which was discovered in 1982, was the
previous fastest known radio pulsar with a spin rate of
642 Hz [14]; 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 oscil-
lations, Chakrabarty et al. [15] claimed a cutoff limit of
νmax = 760 Hz (95% confidence); A more recent analy-
sis, which added two more pulsars to the sample, found
νmax = 730 Hz [16].
At first sight, one might conclude that mass shedding
or radial instability sets νmax, and that it is just above
the record ν = 716 Hz determined for PSR J1748-2446ad.
However, the nuclear equations of state consistent with
this picture all have rather large radii ≈ 16 − 17 km
for non-rotating 1.4 M⊙ models; see Table 1 in Cook et
al. [10]. For these equations of state, the r-mode insta-
bility should lead to ν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 de-
tection of the 716 Hz rotator is consistent with accretion
spin-up mitigated by the r-mode instability only for equa-
tions 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 estab-
lishing such constraints.
Since a physical model to follow the nonlinear phase of
the evolution was initially unavailable, Owen et al. [17]
proposed a simple one-mode evolution model in which
they assumed that nonlinear hydrodynamics effects satu-
rate 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 ampli-
tude (see Eqs. (3.16) and (3.17) in Ref. [17]). They used
this model to study the impact of the r-mode instability
on the spin evolution of young hot neutron stars assum-
ing 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 109 K
and spin down from a frequency close to the Kepler fre-
quency to about 100 Hz in a period of ∼ 1 yr [17].
Most subsequent investigations that did not perform
direct hydrodynamic simulations used the one-amplitude
model of Ref. [17] for studying the r-mode instability.
Levin [18] 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 ∼ 1; he
adopted a modified shear viscosity to match the maxi-
mum LMXB spin frequency of 330 Hz known in 1999.
Levin found that the neutron star followed a cyclic evo-
lution in the Ω − 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 (∼ 1). At this
point the spin and r-mode amplitude evolution equations
are changed, following the prescription of Ref. [17] to en-
sure 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 gravita-
tional radiation. However, following work by Arras et al.
[19] showing that nonlinear effects become significant at
small r-mode amplitude, Heyl [20] 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 α ≈ 10−5, and that
this would make some of the fastest spinning LMXBs in
our galaxy detectable by interferometers on Earth.
Jones [21] and Lindblom and Owen [22] 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 neu-
tron stars. While this additional high viscosity coefficient
could eliminate the instability altogether in newly born
neutron stars [21, 22, 23, 24], Nayyar and Owen [24] pro-
posed that it would enhance the probability of detection
of gravitational radiation from LMXBs by blocking the
thermal runaway.
The cyclic evolution found by Levin [18] and gener-
alized by Heyl [20] arises when shear or boundary layer
viscosity dominates the r-mode dissipation. In the evo-
lutionary picture of Nayyar and Owen [24], the r-mode
first becomes unstable at a temperature where shear
and boundary layer viscosity dominate, but the result-
ing 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 stabil-
ity curve at a point where the viscosity is an increas-
ing function of temperature. Such scenarios were stud-
ied by Wagoner [1] for hyperon bulk viscosty with low
hyperon superfluid transition temperature; similar evo-
lution was found for strange stars by Andersson, Jones
and Kokkotas [25]. 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. [26, 27, 28] for modes with n < 30, and hence in this
regime nonlinear effects may not play a role.
Schenk et al. [29] 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 nonlin-
ear perturbation theory via three-mode couplings. This
assumption was tested by Arras et al. [19] and Brink
et al. [26, 27, 28]. Arras et al. proposed that a turbu-
lent 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 in-
teraction of about 5000 modes via approximatively 1.3
million couplings of the 109 possible couplings among
the modes with n ≤ 30. The couplings were restricted to
mode triplets with a fractional detuning δω/(2Ω) < 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 thresh-
old 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 cal-
culation 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 dis-
sipates via boundary layer effects [30] and bulk viscos-
ity. 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 vary-
ing the free parameters of our model. In particular, we
vary: (1) the slippage factor Sns, which regulates the
boundary layer viscosity, between 0 and 1 (see for exam-
ple [31, 32, 33] 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 fdU, 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 Tc, which is believed
to be between 109− 1010 K (We use a single, effective Tc
rather than modelling its spatial variation.) We focus on
Tc & 5×109 K for which nonlinear effects are important.
For low Tc . 3 × 109 K, Wagoner [1] showed that the
evolution reaches a steady state at amplitudes below the
lowest parametric instability threshold found by Brink
et al. [28]. It is important to note that all our evolu-
tions 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. [28]. We evolve the coupled equa-
tions for the three-mode system numerically in conjunc-
tion 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 ob-
serve that the mode amplitudes settle into a series of
quasi-stationary states within a period of a few years af-
ter 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 equa-
tions (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 evo-
lution equations (six for the real and imaginary parts of
the mode amplitudes, which are complex [29]; one for the
spin, and one for the temperature) to two equations gov-
erning 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 equilib-
rium. This can be predicted by finding the curve where
the viscous heating by the three modes balances the neu-
trino cooling, referred to below as the Heating = Cooling
(H = C) curve. TheH = 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) un-
dergo 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
106 yr, and subsequently enter a regime of thermal run-
away. On the other hand, if the star cannot find a state
of thermal equilibrium, then it enters a regime of ther-
mogravitational runaway within a few hundred years of
crossing the r-mode stability curve. When this happens,
the r-mode amplitude increases beyond the second para-
metric 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 fre-
quencies and couplings are computed for an incompress-
ible 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 in-
volves 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 it-
self 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 vis-
cosity, 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 fea-
tures of our calculations ought to be robust. More sophis-
ticated 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 nu-
merous 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 viscos-
ity is relatively unimportant, soon after the star becomes
unstable its r-mode amplitude passes its first paramet-
ric 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 out-
line here (parametric instability, quasisteady evolution,
competition between heating and cooling) ought to ap-
ply 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 rota-
tional 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. IVA
discusses cyclic evolution in more detail. An evolution
that leads to an equilibrium steady state is presented
next in Sec. IVB. The two types of thermal runaway are
then discussed in Sec. IVC. 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
A. 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 τ for the time variable with
dτ = Ω dt Eq. (2.49) of Schenk et al. [29] can be rewritten
= iω̃αCα +
2iω̃ακ̃√
CβCγ , (1)
= iω̃βCβ −
2iω̃βκ̃√
= iω̃γCγ −
2iω̃γκ̃√
Here the scaled frequency ω̃j is defined to be ω̃j = ωj/Ω,
the dissipation rates of the daughter modes are γβ and γγ ,
γα is the sum of the driving and damping rates of the r-
mode γα = γGR−γαv, and the dimensionless coupling is
κ̃ = κ/(MR2Ω2). These amplitude variables are complex
and can be written in terms of the variables of Ref. [29]
as Cj(t) =
Ω(t)cj(t) (see Appendix A for a derivation
of Eqs. 1). The index j loops over the three modes j =
α, β, γ, where α labels the r-mode or parent mode and β
and γ 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 assump-
tion 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. [28]
for a full derivation). This amplitude is the parametric
instability threshold
|Cα|2 =
4κ̃2ω̃βω̃γΩ
1 + Ω2
γβ + γγ
, (2)
where the fractional detuning is δω̃ = ω̃α − ω̃β − ω̃γ .
Thorough explorations of the phase space of damped
three-mode systems were performed by Dimant [34] and
Wersinger et al. [35].
For the three modes at the lowest parametric instabil-
ity threshold, ω̃α ≈ 0.66, ω̃β ≈ 0.44, ω̃γ ≈ 0.22, κ̃ ≈ 0.19
and |δω̃| ≈ 3.82 × 10−6. Note that ω̃ is twice the w of
Brink et al. [26, 27, 28]. Here β labels the mode with
n = 13,m = −3 and γ labels the n = 14,m = 1 mode.
The amplitude the r-mode has to reach before exciting
these two daughter modes is |Cα| ≈ 1.5× 10−5
Ω [28].
We next rescale the rotational phase τ by the fractional
detuning as τ̃ = τ |δω̃| and the mode amplitudes by
|Cα|0 =
|δω̃|
ω̃βω̃γ
, |Cβ |0 =
|δω̃|
ω̃αω̃γ
, (3)
|Cγ |0 =
|δω̃|
ω̃βω̃α
which for the r-mode is, up to a factor of
Ω/Ωc,
the no-damping limit of the parametric instability thresh-
old below which no oscillations will occur. The coupled
equations become
|δω̃|
C̄α +
|δω̃|Ω̃
C̄α −
C̄βC̄γ , (4)
|δω̃|
C̄β −
|δω̃|Ω̃
C̄β −
C̄αC̄
|δω̃|
C̄γ −
|δω̃|Ω̃
C̄γ −
C̄αC̄
with C̄j = Cj/|Cj |0 and γ̃j = γj/Ωc being the newly
rescaled amplitudes and dissipation/driving rates, re-
spectively.
1. Quasi-Stationary Solution
In terms of amplitudes and phase variables Cj =
|Cj |eiφj Eqs. (4) can be rewritten as
d|C̄α|
Ω̃|δw̃|
|C̄α| −
sinφ|C̄β ||C̄γ |
, (5)
d|C̄β |
= − γ̃β
Ω̃|δw̃|
|C̄β |+
sinφ|C̄α||C̄γ |
d|C̄γ |
Ω̃|δw̃|
|C̄γ |+
sinφ|C̄α||C̄β |
|δω̃|
− cosφ
|C̄β ||C̄γ |
|C̄α|
− |C̄α||C̄γ |
|C̄β |
− |C̄β ||C̄α|
|C̄γ |
where we have defined the relative phase difference as
φ = φα − φβ − φγ . These equations have the stationary
solution
|C̄α|2 =
4γ̃βγ̃γ
Ω̃|δω̃|2
tan2 φ
, (6)
|C̄β |2 =
4γ̃αγ̃γ
Ω̃|δω̃|2
tan2 φ
|C̄γ |2 =
4γ̃αγ̃β
Ω̃|δω̃|2
tan2 φ
tanφ =
γ̃β + γ̃γ − γ̃α
Ω̃|δω̃|
Note that in the limit in which γβ+γγ >> γα the station-
ary solution for the r-mode amplitude |Cα| is the same
as the parametric instability threshold.
B. Temperature and Spin Evolution
The spin evolution equation is obtained from conser-
vation of total angular momentum J , where
J = IΩ + Jphys. (7)
Following Eq (K39-K42) of Schenk et al. [29] the physical
angular momentum of the perturbation can be written as
ΩJphys =
C⋆BCA
d3xρ[(Ω̂× ξ⋆B) · (Ω̂× ξA) (8)
− i (ω̃A + ω̃B)
ξ⋆B · (Ω̂× ξA)].
Since the eigenvectors ξA ∝ eimAφ the cross-terms will
vanish for modes with different magnetic quantum num-
bers m as
ei(mA−mB)φdφ = 0 for mA 6= mB. Eq. (8)
can be re-written for our triplet of modes as
Jphys = MR
2(kαα|Cα|2 + kββ|Cβ |2 + kγγ |Cγ |2), (9)
where kαα is defined as
kαα =
d3xρ[(Ω̂×ξ⋆α) · (Ω̂×ξα)− iω̃αξ⋆α · (Ω̂×ξα)]
and similarly for kββ and kγγ . In terms of the scaled vari-
ables C̄j = Cj/|Cj |0 (with |Cj |0 defined in Eq. (3)) the
angular momentum of the perturbation can be written
Jphys =
MR2Ωc|δω̃|2
(4k̃)2ω̃αω̃βω̃γ
(kαα|C̄α|2ω̃α (11)
+kββ|C̄β |2ω̃β + kγγω̃γ |C̄γ |2).
We chose the same normalization for the eigenfuctions
as Refs. [19, 26, 27, 28, 29] so that at unit amplitude all
modes have the same energy ǫα = MR
2Ω2. The energy
of a mode α is Eα = MR
2Ω2|cα|2 = MR2Ω|Cα|2. The
rotating frame energy is the same as the canonical energy
and physical energy [29]. The canonical angular momen-
tum and the canonical energy of the perturbation satisfy
the general relation Ec = −(ω/m)Jc [3].
Angular momentum is gained because of accretion and
lost via gravitational waves emission
= 2γGRJc rmode + Ṁ
GMR, (12)
where Jc rmode = −(mα/ωα)ǫα|cα|2 = −3MR2Ω|cα|2 =
−3MR2|Cα|2. Eq. (12) can be rewritten in terms of the
scaled variables C̄j as
= −6γ̃GR
MR2Ωc|δω̃|
(4k̃)2ω̃βω̃γ
|C̄α|2 +
ΩcΩ̃|δω̃|
. (13)
Thermal energy conservation gives the temperature evo-
lution equation
C(T )
2Ejγj +KnṀc
2 − Lν(T ), (14)
= 2MR2Ω(γα v|Cα|2 + γβ |Cβ |2
+ γγ |Cγ |2) +KnṀc2 − Lν(T ).
The three terms on the right hand side of the equa-
tion represent viscous heating, nuclear heating and neu-
trino cooling. The specific heat is taken to be C(T ) ≈
1.5 × 1038 T8 erg K−1, where T = T8 × 108 K. Nu-
clear heating occurs because of pycnonuclear reactions
and neutron emission in the inner crust [36]. At large
accretion rates such as that of the brightest LMXBs of
Ṁ ≈ 10−8M⊙/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 neu-
trino emission is assumed to regulate the temperature of
the star [36, 37]. The nuclear heating constant is taken
to be Kn ≈ 1×10−3 [36]. Following Ref. [1], we take the
neutrino luminosity to be
Lν = LdUT
8RdU(T/Tp) + LmUT
8RmU(T/Tp) (15)
+ Le−iT
8 + Ln−nT
8 + LCpT
where the constants for the modified and direct URCA re-
actions are defined by LmU = 1.0×1032 erg sec−1, LdU =
fdU × 108LmU [38, 39], and the electron-ion, neutron-
neutron neutrino bremsstrahlung and Cooper pairing of
neutrons are given by Le−i = 9.1 × 1029 erg sec−1 [36],
Ln−n ≈ 0.01LmU, LCp = 8.9 × 1031 erg sec−1 [40]. The
fraction of the star fdU that is above the density thresh-
old for direct URCA reactions is in general dependent on
the equation of state [41] and in this work we treat fdU
a free parameter with values between 0 and 1.
The proton superfluid reduction factors for the modi-
fied and direct URCA reactions are taken from Ref. [39]
(see Eqs. (32) and (51) in Ref. [39]):
RdU(T/Tp) =
0.2312 +
(0.76880)2 + (0.1438v)2
× exp
3.427−
(3.427)2 + v2
RmU(T/Tp) =
0.2414 +
(0.7586)2 + (0.1318v)2
× exp
5.339−
(5.339)2 + (2v)2
where the dimensionless gap amplitude v for the singlet
type superfluidity is given by
1.456− 0.157
+ 1.764
. (17)
Similar to Ref. [1], we use Tp = 5.0× 109 K. In terms of
the scaled variables Eq. (14) becomes
C(T )
2MR2Ω2c |δω̃|
(4κ̃)2ω̃αω̃βω̃γ
(ω̃αγ̃α v|C̄α|2 + ω̃βγ̃β |C̄β |2(18)
+ω̃γγ̃γ |C̄γ |2) +
KnṀc
2 − Lν(T )
ΩcΩ̃|δω̃|
C. 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.,
dCi/dτ̃ ≈ 0, the spin and thermal evolution equations
can be rewritten by taking J ≈ IΩ and using Eqs. (6) in
Eq. (13).
= − 6γ̃GR
Ω̃2|δω̃|
γ̃β γ̃γ
4k̃2Ĩω̃βω̃γ
tan2 φ
MR2ĨΩ̃|δω̃|
where Ĩ = I/(MR2). The thermal evolution of the sys-
tem is given by
C(T )
2MR2Ω2c
(4κ̃)2ω̃αω̃βω̃γ
γ̃αγ̃β γ̃γ
Ω̃|δω̃|
ω̃αγ̃α,v
+ ω̃β (20)
+ω̃γ)
tan2 φ
KnṀc
2 − Lν(T )
ΩcΩ̃|δω̃|
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 evolu-
tion very well throughout the unstable regime. These
equations are a minimal physical model for the effects of
nonlinear coupling on r-mode evolution.
D. Sources of Driving and Dissipation
The damping mechanisms are shear viscosity, bound-
ary layer viscosity and hyperon bulk viscosity; for modes
j = α, β, γ we write
γj v(Ω, T ) = γj sh(T ) + γj bl(Ω, T ) + γj hb(Ω, T ). (21)
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 β and
n = 14,m = 1 labeled as γ) is affected only by the vis-
cous damping. Brink et al. [26, 27, 28] determined that
this pair of modes is excited at the lowest parametric
instability threshold. Their model uses the Bryan [42]
modes of an incompressible star, which has the advan-
tage 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 (Γ1 → ∞ in
Eq. (29))
For our benchmark calculations, we adopt the neutron
star model of Owen et al. Ref. [6] (n = 1 polytrope,
M = 1.4M⊙, Ωc = 8.4 × 103 rad sec−1 and R = 12.53
km) and use their gravitational driving rate and shear
viscous damping rate for the r-mode
γGR(Ω) ≃
sec−1, (22)
γα sh(T ) ≃
where τsh = 2.56 × 106 sec. (In Sec. V we consider ap-
proximate 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
γβ sh(T ) ≃ 3.48× 10−4 sec−1
, (23)
γγ sh(T ) ≃ 4.52× 10−4 sec−1
The geometric contribution γsh/η of the individual modes
increases significantly with the degree n of the mode scal-
ing approximatively like n3 for large n (see Eq. (29) of
Brink et al. [27] 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 magni-
tude larger than that of the r-mode.
The damping due to boundary layer viscosity is calcu-
lated using Eq. (4) of Ref. [30],
γα bl(T,Ω) ≃ 0.009 sec−1 S2ns
, (24)
γβ bl(T,Ω) ≃ 0.028 sec−1 S2ns
γγ bl(T,Ω) ≃ 0.021 sec−1S2ns
Analogous to Wagoner [1], we allow the slippage fac-
tor Sns to vary. The slippage factor is defined by Refs.
[1, 31, 45] to be S2ns = (2S
n + S
s )/3, with Sn being the
fractional difference in velocity of the normal fluid be-
tween the crust and the core [31] and Ss the fractional
degree of pinning of the vortices in the crust [45]. Note
that γβ bl and γγ bl are both greater than 2 × γα bl and
can easily be comparable to γ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 [24].
This rate is taken to have a form similar to that taken
by Wagoner [1]
γα hb = fhb
t−20α τ(T )Ω̃
1 + (ω̃αΩτ(T ))2
, (25)
and for the daughter modes
γβ hb = fhb
t−20β τ(T )ω̃
1 + (ω̃βΩτ(T ))2
, (26)
and similarly foe γγ hb. The relaxation timescale
τ(T ) =
Rhb(T/Tc)
The reduction factor is taken to be the product of two
single-particle reduction factors [23, 24]
Rhb single(T/Tc) =
a5/4 + b1/2
0.5068−
0.50682 + y2
where a = 1 + 0.3118y2, b = 1 + 2.556y2 and y =
1.0− T/Tc(1.456 − 0.157
Tc/T + 1.764Tc/T ). The
constants t1 ≈ 10−4 sec and t0α ≈ 0.00058 sec are found
by fitting the results of Ref. [24]. The factor fhb allows
for physical uncertainties; we take fhb = 1 throughout
the body of the paper since Tc , which enters γj hb ex-
ponentially, 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 Γ−21 by
− ĖB j =
ξj · ∇p
. (29)
This approximation was proposed by Cutler and Lind-
blom [43] and adopted by Kokkotas and Stergioulas [44]
for the r-mode and by Brink et al. [27] for the inertial
modes. The adiabatic index Γ1 is regarded as a parame-
ter; we use Γ1 ≈ 2. The damping rate is
γj hb = −
ĖB j
, (30)
where ǫ = MR2Ω2 is the mode’s energy in the rotat-
ing frame at unit amplitude and j = β, γ. Using this
procedure, we calculate
t0β ≈ 1.4× 10−5 sec, (31)
t0γ ≈ 1.0× 10−5 sec.
III. SUMMARY OF RESULTS
Fig. 1(a) shows possible evolutionary trajectories of a
neutron star in the angular velocity-temperature Ω̃− T8
plane, where T = T8×108 K is the core temperature, and
Ω̃ = Ω/Ωc = Ω/
πGρ̄ with ρ̄ the mean density of the
neutron star. Fig. 1(b) displays the regions in fdU − Sns
in which the trajectories occur. Here fdU represents the
fraction of the star that is above the density threshold
for direct URCA reactions and Sns is the slippage factor
that reduces the relative motion between the crust and
the core taking into account the elasticity of the crust
[31]. The stability regions are shown at fixed hyperon
superfluidity temperature, Tc = 5.0 × 109 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
T8 = T8 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 Ω because of viscous dissipation
in all three modes.
Phase 3. Several trajectories are possible depending on
FIG. 1: (a)Typical trajectories for the four observed evolu-
tion scenarios are shown in the Ω̃ - T8 phase space, where
Ω̃ = Ω/Ωc. The dashed lines (H = C curves) represent
the points in the Ω̃ − T8 phase space where the dissipative
effects of the heating from the three-modes exactly compen-
sate the neutrino cooling for the given set of parameters (Sns,
fdU, Tc, ...) of each evolution. (b)The corresponding sta-
bility regions for which these scenarios occur are plotted at
fixed hyperon superfluidity temperature Tc = 5.0 × 10
while varying fdU and Sns. The position of the initial angu-
lar velocity and temperature (Ω̃in, T8 in) with respect to the
maximum of this curve determines the stability of the evo-
lution. (I) Ω̃in > Ω̃H=C max. Trajectory R1. Fast Runaway
Region. After the r-mode becomes unstable the star heats
up, does not find a thermal equilibrium state and continues
heating up until a thermogravitational runaway occurs. (II)
Ω̃in < Ω̃H=C max. The evolutions are either stable or, if there
is a runaway, it occurs on timescales comparable to the ac-
cretion timescale. The possible trajectories are (1)Trajectory
C. Cycle Region. (2) Trajectories S1 and S2. Steady State
Region. (3) Trajectory R2. Slow Runaway Region.
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 [24]. 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 T8 in at constant Ω̃ and
then the cycle repeats itself. At Tc = 5.0 × 109 K this
scenario occurs for values of Sns < 0.50 and large enough
values of fdU. However, if Tc is larger, the cycle region
in the fdU-Sns phase space increases dramatically (see
Fig. 9(a)). Note that our cycles are different from those
obtained by Levin [18] in that the spin-down phase
does not start when the r-mode amplitude saturates
(or in our case when it reaches the parametric insta-
bility 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 ∼ 105 and so the part
of the cycle in which the r-mode is unstable also lasts
longer than in Ref. [18]. Also, our cycles are narrow.
During spin-down the temperature changes by less than
20 % and Ω̃ changes by less than 10% of the initial value.
(See Sec. 2 for a detailed example.)
(ii) Steady State. For small Sns and large enough
fdU (fdU & 5 × 10−5, Sns . 0.04; see Fig. 1(b)) the
star evolves towards an Ω̃ 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 Sns and very small fdU
(Sns . 0.03, fdU < 5×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 Tc ∼ 1010, or for models with no hyperons at all,
there would be no runaway region (See Fig. 9(a) for an
fdU − Sns scenario space with a larger Tc = 6.5× 109 K
where the fast runaway region has shrunk dramatically
FIG. 2: Two cyclic trajectories in the Ω̃ − T8 plane are dis-
played for a star with Tc = 5.0 × 10
9 K and (a) fdU = 0.15
and Sns = 0.10, and (b) fdU = 0.142 and Sns = 0.35, which is
close to the border between the stable and unstable region (see
Fig. 1(b)). The thick solid line labeled as the Heating = Cool-
ing (H = C) curve is the locus of points in this phase space
where the neutrino cooling is equal to the viscous heating due
to the unstable modes. The other solid line representing the
r-mode stability curve is defined by setting the gravitational
driving rate equal to the viscous damping rate. The part of
the curve that decreases with T8 is dominated by boundary
layer and shear viscosity, while the part of the curve that has
a positive slope is dominated by hyperon bulk viscosity. In
portion a1 → b1 of the trajectory the star heats up at con-
stant Ω̃. Part b1 → c1 represents the spin down stage, which
occurs when the viscous heating is equal to the neutrino cool-
ing. c1 → d1 shows the star cooling back to the initial T8.
Segment d1 → a1 displays the accretional spin-up of the star
back to the r-mode stability curve. The cycle a2 → d2 pro-
ceeds in the same way. This cycle is close to the peak of the
H = C curve. Configurations above this peak will run away.
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 Ṁ =
10−8M⊙/yr and Tc = 5.0× 109 K.
A. 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) Sns = 0.10 and fdU = 0.15 and (C2) Sns = 0.35
and fdU = 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 T8 in ≈ 3.29 for both cases. The initial
Ω is the angular velocity that corresponds to this tem-
perature on the r-mode stability curve, which differs for
the different Sns (Ω̃in = 0.183 for C1 and Ω̃in = 0.288 for
Figs. 2(a) and (b) display the cyclic evolution for tra-
jectories C1 and C2 of Fig. 1(b). In leg a1 → b1 of the
trajectory the r-mode and, once the r-mode amplitude
increases above the first parametric instability thresh-
old, the two daughter modes it excites, viscously heat
up the star until point b1 when the neutrino cooling bal-
ances the viscous dissipation. This part of the evolu-
tion occurs at constant angular velocity over a period
of theat−up ≈ 100 yr and a total temperature change
(∆T )a1−b1 ≈ 0.80 (≈ 24% of T8 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 b1 → c1 of the
trajectory as it spins down and cools down back to the r-
mode stability curve. This spin-down stage lasts a time
tspin−downb1−c1 ≈ 23, 000 yr that is much longer than
the heat-up period. This timescale is very sensitive to
changes in the slippage factor and can reach 106 yr for
smaller values of Sns 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%, (∆Ω̃)b1−c1 ≈ 0.0066. The temperature also changes
by only about 2%, (∆T8)b1−c1 ≈ 0.08 in this spin-down
period. Segment c1 → d1 represents the cooling of the
star to the initial temperature on a timescale of ∼ 2, 000
yr. In part d1 → a1 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 tra-
jectory is the longest-lasting one, taking ≈ 200, 000 yr
at our chosen Ṁ of 10−8M⊙yr
−1. The cycle C2 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 there-
fore close to the peak of the H = C curve. Configura-
tions above this peak (e.g., with the same fdU and higher
Sns) 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 un-
stable and initially grows exponentially. Once it has in-
creased 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 exponen-
tially while the r-mode decreases. Similarly, when daugh-
ter 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 ampli-
tude 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 thresh-
old 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.
B. 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 re-
stricted to stars with small slippage factor (Sns . 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 sta-
bility curve the gravitational driving rate is greater than
the viscous damping rate and the r-mode grows exponen-
tially 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 evo-
lution and after a few oscillations the three modes settle
into their quasi-stationary states, which change only sec-
ularly as the spin and temperature of the star evolve.
The modes heat the star viscously at constant Ω̃ in seg-
ment a → b of the trajectory for theat−up ≈ 1, 100 yr.
At point b, the star reaches a state of thermal balance.
In leg b → c the star continues its evolution in thermal
equilibrium and slowly spins up due to accretion until
the angular velocity evolution also reaches an equilib-
5.12 5.13
τ [x 10
0.422
0.424
0.426
0.428
0.432
|cα Th|
0.4 0.5 0.6 0.7 0.8 0.9 1
τ [x 10
0 0.5 1 1.5 2 2.5 3 3.5 4
τ [x 10
0.425
0.435
FIG. 3: (a)The amplitudes of the r-mode |Cα| and of the
n = 13, m = −3 and n = 14, m = 1 inertial modes |Cβ | and
|Cγ | are shown as a function of time for a star that executes
a cyclic evolution (same parameters as in Fig. 2). The low-
est parametric instability threshold is also displayed. (b)The
ratio of the r-mode amplitude to the parametric instability
threshold is plotted as a function of time. It can be seen that
once the r-mode crosses the parametric instability threshold
it remains close to it for the rest of the evolution. (c)The
parametric instability threshold is displayed as a function of
time. Its value changes as the angular velocity and tempera-
ture evolve.
FIG. 4: The trajectory of a neutron star in the Ω̃− T8 phase
space is shown for a model with Tc = 5.0× 10
9 K, fdU = 0.03
and Sns = 0.03 that reaches an equilibrium steady state. The
star spins up until it crosses the r-mode stability curve and
the r-mode becomes unstable. The r-mode then quickly grows
to the first parametric instability threshold and excites the
daughter modes. In leg a → b of the trajectory the star is
viscously heated by the mode triplet until the system reaches
thermal equilibrium. Segment b → c shows the star contin-
uing to heat and spin up in thermal equilibrium until the
accretion torque is balanced by the gravitational radiation
emission. The r-mode stability curve represents the points
in phase space where the viscous driving rate is equal to the
gravitational driving rate. The H=C curve is the locus of
points where the viscous dissipation due to the mode triplet
balances the neutrino cooling.
FIG. 5: The (Ω̃, T8) initial values (region delimited by the
solid line) that lead to equilibrium steady states and their
corresponding final steady state values (region enclosed by
the dashed line) are shown. Since both the initial and final
values of T8 are low, these evolutions are roughly independent
of Tc.
rium. The timescale to reach an equilibrium steady state
is tsteady ≈ 3.5× 106 yr for this set of parameters.
Fig. 5 displays the possible initial values for the angu-
lar velocity Ω̃ and temperature T8 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 Ω also tend to be small. Thus, although an interest-
ing physical regime, this case is most likely not relevant
to recycling by accretion to create pulsars with spin fre-
quencies as large as 716 Hz. Note that a steady state
can be achieved when Sns = 0. This is the probable end
state of the problem first calculated by Levin [18]. The
reason we do not find a cycle at low Sns is twofold: (1)
the shear viscosity we are using is lower (shear viscosity
in Ref. [18] is amplified by a factor of 244), and (2) the
nonlinear couplings keep all mode amplitudes small.
C. 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 run-
away 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 contin-
ues to spin up by accretion until it climbs to the peak of
the H = C curve, thermal equilibrium fails and a run-
away occurs.
1. Fast Runaway
A typical trajectory of a star that goes through a rapid
thermal runaway is displayed in Fig. 6. This star has
Sns = 0.25 and fdU = 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 (Ω,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. [24].
Fig. 6(a) and (b) compare both the temperature evolu-
tion and the trajectory in the Ω̃−T8 plane of the star for a
0 0.5 1 1.5 2 2.5 3 3.5
τ [x 10
Full Amplitudes-Ω-T Evolution
Steady State Evolution
2nd Parametric Instability Threshold
3 3.5 4 4.5 5 5.5 6 6.5
r-mode Stability Curve
Full Amplitude-Ω-T Evolution
Steady State Evolution
2nd Parametric Instability Threshold
FIG. 6: This plot compares the full evolution resulting from
solving Eqs. (4),(13),(18) with the reduced Ω − T evolution
that assumes the amplitudes go through a series of steady
states Eqs. (19)-(20) for a model with Tc = 5.0 × 10
fdU = 0.058 and Sns = 0.25. (a) The temperature is dis-
played as a function of time for the two different methods.
(b) The angular velocity Ω̃ = Ω/Ωc is shown as a function
of temperature. The evolution occurs at constant spin fre-
quency. It can be seen that the steady-state amplitude ap-
proximation is extremely good. The ‘X’ shows the point at
which the r-mode crosses its second lowest parametric insta-
bility threshold, where additional dissipation would become
operative.
simulation solving the full set of equations to a simulation
that assumes quasi-stationary solutions for the three am-
plitudes and evolves only the angular velocity and tem-
perature 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 re-
duced 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 [24]. However, the r-mode crosses its second low-
FIG. 7: The trajectory of a neutron star in the Ω̃ − T8
phase space is shown for a model with Tc = 5.0 × 10
fdU = 4.0 × 10
−5 and Sns = 0.02 that goes through a slow
thermogravitational runaway. Portion a → b of the trajectory
shows the mode triplet heating up the neutron star through
boundary layer and shear viscosity until the system reaches
thermal equilibrium. Segment b → c represents the accre-
tional spin-up of the star in thermal equilibrium. The dotted-
dashed line is the locus of points where the viscous dissipation
of the mode triplet is equal to the neutrino cooling, and is la-
beled as the H = C curve. The star reaches the maximum of
this curve and fails to reach an equilibrium between the ac-
cretion torque and gravitational emission. It then continues
heating at constant angular velocity and crosses its second
lowest parametric instability threshold, at which point more
modes would need to be included to make the evolution accu-
rate. Eventually the star reaches the r-mode stability curve
again.
est 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.
2. Slow Runaway
In this section we examine evolutions in which the neu-
tron star has both a very small slippage factor, Sns .
0.03, and only a small percentage of the star is above the
threshold for direct URCA reactions, fdU < 5 × 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 insta-
bility 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 Ω in part a → b of the trajec-
tory on a timescale of about 5, 000 yr. At point b, it
FIG. 8: The spin-down timescale is shown as slippage fac-
tor Sns and fraction of the star subject to direct URCA fdU
for cyclic evolutions are varied for a fixed hyperon critical
temperature of Tc = 5.0 × 10
9 K. This timescale dominates
the heat-up timescale and hence represents the time the star
spends above the r-mode instability curve. It increases as the
viscosity is lowered and the star gets closer to the steady state
region.
reaches thermal equilibrium. In leg b → c of the tra-
jectory, the star continues its evolution by ascending the
H = C curve and spinning up because of accretion for
about 2 × 106 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 expo-
nentially 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 instabil-
ity 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 Sns and fdU are var-
ied. For large enough values of Sns the boundary layer
viscosity dominates. In this region of phase space the
spin-down timescale can be approximated by
tspin−down =
dΩ̃ (32)
FIG. 9: (a)The stability regions are plotted at fixed hyperon
superfluidity temperature Tc = 6.5×10
9 K, while varying fdU
and Sns. The steady state region remains roughly the same
as in Fig. 1(b), the slow run-away region disappears, and the
cycle region increases dramatically while shrinking the fast-
runaway region. (b) The spin-down timescale is shown for the
cyclic evolutions in part (a).
Ĩτ0GR
(4κ̃)2ω̃βω̃γ
|δω̃|2
|C̄α|2
<Ω̃>6
≈ 250 yr ∆νkHz
<νkHz>7
M1.4R
|cthα |
where M1.4 = M/(1.4M⊙), R6 = R/(10
6cm),
νkHz = ν/1kHz, Ĩ = 0.261 [17], the r-mode am-
plitude at its parametric instability threshold |cthα | ≈
|δω̃|/(4κ̃
ω̃βω̃γ) ≈ 1.5×10−5, and C̄α =
Ω̃|cα|/|cα|th.
This approximation agrees with spin-up timescales ob-
tained from our simulations to ∼ 25%.
The maximum ν is approximately the same as the ini-
tial 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
νmax ≈ 800Hz
M1.4R6
)4/11
. (33)
Thus, the spin-down timescale is very sensitive to
the slippage factor tspin−down ∝ S−24/11ns (∆νkHz/νkHz).
The dependences on fdU and accretion rate Ṁ are
much weaker; a rough approximation, obtained by
matching direct URCA cooling and nuclear heat-
ing, is T8 in ∝ Ṁ1/6f−1/6dU R
1.4 , and νmax ∝
dU Ṁ
−1/33R
−34/99
1.4 . The gravitational
wave amplitude measured at distance d [46, 47] is
h ≈ 1.6R
τ0GRc
Ω̃3|cα| (34)
≈ 3× 10−25
10kpc
M1.4R
Taking ν ≈ νmax gives
h ∝ S12/11ns M
−1/33
1.4 R
dU Ṁ
−1/11. (35)
The maximum distance at which sources could be de-
tected by Advanced LIGO interferometers, assuming
hmin = 10
−27, [46] is
dmax ≈ 3Mpc
10−27
M1.4R
|cthα |
≈ 1.5Mpc
10−27
S12/11ns M
−1/11
1.4 R
21/11
× T−6/118
|cthα |
Eqs. (33) and (36) 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
ν & 700 Hz, we must require (Sns/M1.4R6
T8in)
4/11 & 1.
Assuming this to be true, dmax . 1-1.5 Mpc. However,
tspin−down ≈ 1000 yr at most, making detection unlikely
for any given source. Moreover, unless Sns can differ sub-
stantially from one neutron star to another, only those
with ν 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 Sns . 0.50, assuming Tc ≈ 5.0× 109
K; Eq. (33) then implies ν . 450 Hz. This would restrict
detectable gravitational radiation to galactic sources, al-
though the duration of the unstable phase could be
longer.
Within the context of our three mode calculation,
Sns > 0.50, which is needed for explaining the fastest pul-
sars, would imply fast runaway. There are two possible
resolutions to this problem. One is that including addi-
tional modes prevents the runaway; we shall investigate
this in subsequent papers. The second is that Tc is larger,
or that neutron stars do not contain hyperons (e.g., be-
cause they are insufficiently dense). Fig. 9(a) shows the
same phase plane as Fig. 1(b) but with Tc = 6.5×109 K,
and Fig. 9(b) shows the results for tspin−down analogous
to Fig. 8. Larger Tc permits spin cycles for higher values
of Sns (and hence ν), but the time spent in the unstable
regime is shorter.
VI. CONCLUSIONS
In this paper, we model the nonlinear saturation of un-
stable 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 en-
ergy 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 slip-
page 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 evolu-
tions we find that the mode amplitudes quickly settle
into a series of quasi-stationary states that can be calcu-
lated algebraically, and depend weakly on angular veloc-
ity 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 evo-
lution 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 pos-
sible 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 Ω 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 Sns ∼ 1) and 106 yr (small Sns ∼ 0.05). Our
cycles are different from those previously found by Ref.
[18] in that our amplitudes remain small, ∼ 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 ∼ 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 thermo-
gravitational runaway on a timescale of ∼ 106 yr. These
scenarios occur for very low viscosity (Sns . 0.04). Al-
though theoretically interesting, they do not allow for
very fast rotators of ∼ 700 Hz.
Alternatively, if the star does not reach thermal equi-
librium, 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 [24]. 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 Tc & 5× 109 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 ro-
tation requires large dissipation, as has long been recog-
nized [18, 30] and these models can only achieve ν & 700
Hz if boundary layer viscosity is very large. Alterna-
tively, at lower Tc . 3 × 109 K, large rotation rates can
be achieved at r-mode amplitudes below the first para-
metric instability threshold [1]. Nayyar and Owen found
that increasing the mass of the star for the same equation
of state makes the hyperon bulk viscosity become impor-
tant at lower temperatures [24]. Conceivably, there are
accreting neutron stars with relatively low masses that
have lower central densities and small hyperon popula-
tions. These could evolve as detailed here and only spin
up to modest frequencies. Hyperons could be more im-
portant in more massive neutron stars leading to larger
spin rates and very small steady state r-mode amplitude
as found by Wagoner [1].
Our models imply small r-mode amplitudes of ∼ 10−5
and therefore gravitational radiation detectable by ad-
vanced 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 νmax ≈
800Hz[Sns/(M1.4R6)]
4/11T
−2/11
8 . In order to allow for
fast rotators of & 700 Hz in our models a large bound-
ary layer viscosity with (Sns/M1.4R6
T8in)
4/11 ∼ 1 is
required. Slippage factors of order ∼ 1 lead to time peri-
ods 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 fre-
quencies. We also note that in this model we have con-
sidered only very fast accretors with Ṁ ∼ 10−8M⊙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 [48], superfluid-
ity, 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 quanti-
tatively in a model different from ours in detail, we ex-
pect the qualitative behaviors we have found to be ro-
bust, as they are well described by quasi-stationary mode
evolutions whose slow variations are determined by com-
petitions between dissipation and neutrino cooling, and
accretion spin-up and gravitational radiation spin-down.
In our model, it seems that three mode evolution involv-
ing 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. [26, 27, 28]
that includes accretion spin-up, viscous heating and neu-
trino 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 Éanna
Flanagan for useful discussions. RB would especially
like to thank Jeandrew for useful discussions, encourage-
ment and advice at the beginning of this project, with-
out which the project would not have been started. RB
is very grateful to Gregory Daues for steady encourage-
ment 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 Ω.
The Lagrangian density as given by Eq. (A1) in Schenk
et al. [29] is
L = 1
ξ̇ · ξ̇ + 1
ξ̇ ·B · ξ − 1
ξ ·C · ξ + aext(t) · ξ, (A-1)
where the operators B · ξ = 2Ω× ξ and
ρ(C · ξ)i = −∇i(Γ1p∇jξj) +∇ip∇jξj + ρ∇iδφ (A-2)
− ∇jp∇iξj + ρξj∇j∇iφ+ ρξj∇j∇iφrot
with φrot = −(1/2)(Ω × 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
, dτ = Ωdt. (A-3)
In terms of these new variables the Lagrangian density
can be written as
L̃ = 1
ξ̃′ · ξ̃′ + 1
ξ̃′ · (B̃ · ξ̃) + (
|ξ̃|2 (A-4)
ξ̃ · C̃ · ξ̃ + aext(t)
· ξ̃,
where the primes denote derivatives with respect to τ ,
B̃ = Ω−1B and C̃ = Ω−2C. The momentum canonically
conjugate to ξ̃ is
= ξ̃′ + Ω̂× ξ̃. (A-5)
The associated Hamiltonian density is
B̃ · ξ̃
|ξ̃|2 +
ξ̃ · C̃ · ξ̃ −
· ξ̃.
(A-6)
Hamilton’s equations of motions can be written as
ζ̃′ = T · ζ̃ + F(τ), (A-7)
where
the operator T is T = T0 + T1 with
B̃2 − C̃ − 1
Ω)′′√
F(τ) =
We assume solutions of the form ζ̃(τ,x) = eiω̃tζ̃(x). Spe-
cializing to the case of no forcing term aext = 0 leads to
the eigenvalue equation
(T0 − iω̃)ζ̃(x) = 0. (A-8)
Since the operator T0 is not Hermitian it will have dis-
tinct right and left eigenvectors. Similar to Schenk et
al. [29] we label the right eigenvectors of T as ζ̃A, and
the associated eigenfrequencies as ω̃A = ωA/Ω, and the
eigenvalue equation above becomes
(T0 − iω̃A)ζ̃A(x) = 0. (A-9)
The left eigenvectors χA satisfy
0 − iω̃⋆A)χ̃A = 0, (A-10)
where
B̃2 − C̃
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 be-
tween right and left eigenvectors is
χ̃A, ζ̃B
d3xρ(x)χ̃
A · ζ̃B = δAB. (A-11)
We can expand ζ(τ,x) in this basis as
ζ(τ,x) =
CA(τ)ζA(x), (A-12)
where the coefficients CA are given by the inverse of this
mode expansion
CA(τ) =
χ̃A, ζ̃(τ,x)
. (A-13)
Using Eqs. (B-2,A-9,A-11) in Eq. (A-7) leads to the equa-
tions of motion for the mode amplitudes
C′A − iω̃ACA = g(τ)
(A-14)
+ 〈χ̃A, F (τ)〉 ,
where g(τ) = (
Ω)′′/
Ω. Following Sec. IV of Schenk
et al. [29] we replace the externally applied acceleration
by the nonlinear acceleration given by Eq. (4.2) of Ref.
[29]. The inner product can be written in terms of the
displacement variable ξ̃. The left eigenvectors are
χ̃A =
where τ̃A can be chosen to be proportional to ξ̃A because
they satisfy the same matrix equation.
τ̃A = −iξ̃A/b̃A, (A-15)
which corresponds to Eq. (A-45) in Schenk et al. [29]
with the proportionality constant b̃A = Ω
−1bA =
MR2/ω̃A (also given by Eq. (2.36) of Ref. [29]).
The equations of motion for the mode amplitudes be-
C′A − iω̃ACA =
ig(τ)
d3xξ̃⋆A · ξ̃B(A-16)
κ̃⋆ABCC
where the nonlinear coupling κ̃ABC = κABC/(MR
and κABC is explicitly give by Eq. (4.20) of Ref. [29]. The
g(τ) integral mixes only modes with mA = mB because
of the eimφ dependence of the displacement eigenvectors
ξ̃. (
dφei(mA−mB)φ = 0 if mA 6= mB.) 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.
(A-16) 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. [49], each
of the five variables is expanded about its equilibrium
(Xj)e as follows
Xj(τ̃ ) = {|C̄α|, |C̄β |, |C̄γ |, φ, T8} = (Xj)e[1 + ζj(τ̃ )]
(B-1)
where the perturbation |ζj | << 1 and j = α, β, γ, T . The
expansion leads to a first order differential equation for
each ζj
(γ̃α)e
Ω̃|δω̃|
ζα − ζβ − ζγ −
ζφ (B-2)
(γ̃β)e
Ω̃|δω̃|
ζα − ζβ + ζγ +
(γ̃γ)e
Ω̃|δω̃|
ζα + ζβ − ζγ +
φe tanφe
γ̃α + γ̃β + γ̃γ
Ω̃|δω̃|
−γ̃α − γ̃β + γ̃γ
Ω̃|δω̃|
−γ̃α + γ̃β − γ̃γ
Ω̃|δω̃|
(γ̃α − γ̃β − γ̃γ)e
Ω̃|δω̃|
MR2Ω2c γ̃αγ̃β γ̃γ
2κ̃2ω̃αω̃βω̃γΩ̃|δω̃|C(Te)T8e
tanφ2e
γ̃α v
ζα + ω̃βζβ + ω̃γζγ
+ T8e
+ ω̃β
+ ω̃γ
ΩcΩ̃|δω̃|C(Te)
where the equilibrium amplitudes |Cj |e have been written
in terms of the corresponding driving and damping rates
using Eqs. (6). Eq. (B-2) can be written in matrix form
= Aijζi. (B-3)
Let ζj ∝ exp(λτ̃ ). The determinant ||Aij − λδij || = 0
leads to the eigenvalue equation
λ5 + a4λ
4 + a3λ
3 + a2λ
2 + a1λ+ a0 = 0. (B-4)
The coefficients aj with j = 0, 4 are
a4 = 2 tanφe =
γ̃β + γ̃γ − γ̃α
Ω̃|δω̃|
, (B-5)
tanφ2e
γ̃2β + γ̃
γ + γ̃
(Ω̃|δω̃|)2
+ tanφ2e − 1,
γ̃αγ̃β γ̃γ
(Ω̃|δω̃|)3
tanφ2e
4γ̃αγ̃β γ̃γ
(Ω̃|δω̃|)3
tanφe
+ tanφ
2MR2Ω2c
κ̃2ω̃αω̃βω̃γC(Te)
(γ̃αγ̃β γ̃γ)
(Ω̃|δω̃|)4
tanφe
tanφ2e
4γ̃αγ̃β γ̃γ
(Ω̃|δω̃|)3
tanφe
Ω̃|δω̃|C(Te)
The eigenvalues can be approximated as
λ1,2 ≈ −
− ǫ± i
ǫ2 + w2
, (B-6)
λ3,4 ≈ ǫ± iw,
λ5 ≈ −
where ǫ = (a2 − a3a4)/a4 and w =
a1/a3. The system
is unstable when a2 − a3a4 > 0 or a0 < 0. The first
two eigenvalues will have a negative real part as long as
γ̃β + γ̃γ > γ̃α. If the heating compensates the cooling of
the star a0 ≈ 0 and becomes negative if the star can not
reach thermal equilibrium. The other critical stability
condition a2 − a3a4 = 0 can be written as
Ω̃|δω̃|
[1+Γβ+Γγ−(Γ2β+Γ2γ)−(Γβ−Γγ)2(Γβ+Γγ)] = 0,
(B-7)
where Γβ = γβ/γα and Γγ = γγ/γα. Note that we have
ignored the smaller terms of orderO([γ̃α/(Ω̃|δω̃|)]5). This
condition can be rewritten by defining variables D1 =
Γβ + Γγ and D2 = Γβ − Γγ
2 + 2D1 −D21 −D22 − 2D22D1 = 0. (B-8)
If D2 = 0 then the equation has one solutionD1 = 1+
for D1 > 2, which corresponds to Γ = Γβ = Γγ = 1.37
and matches the result of Wersinger et al. [35]. For the
viscosity we consider (see Sec. II D) a2 − a3a4 < 0.
[1] R. Wagoner, Astrophys. J 578, L63 (2002).
[2] S. Chandrasekhar, Phys. Rev. Lett. 24, 611 (1970).
[3] J. L. Friedman and B. F. Schutz, Astrophys. J. 222, 281
(1978). J. L. Friedman and B. F. Schutz, Astrophys. J.
221, 937 (1978).
[4] N. Andersson, Astrophys. J. 502, 708 (1998).
[5] J. Friedman and S. Morsink, Astrophys. J. 502,
714(1998).
[6] L. Lindblom, B. J. Owen, and S. M. Morsink, Phys. Rev.
Lett 80, 4843 (1998).
[7] L. Bildsten, Astrophys. J. 501, L89 (1998).
[8] N. Andersson, K. D. Kokkotas, N. Stergioulas, Astro-
phys. J. 516, 307 (1999).
[9] N. Andersson, K. Kokkotas, and B. F. Schutz, Astrophys.
J. 510, 846 (1999).
[10] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro-
phys. J 423 L117 (1994).
[11] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro-
phys. J 424 823 (1994).
[12] J. W. T. Hessels et al. Science 311 1901 (2006).
[13] J. E. Grindlay, Science 311, 1876 (2006).
[14] D. C. Backer et al., Nature 300, 615 (1982).
[15] D. Chakrabarty et al., Nature 424, 42 (2003).
[16] D. Chakrabarty, Astron. Soc. Pac. Conf. Series 328, 279
(2005).
[17] B. J. Owen et al., Phys. Rev. D 58, 084020 (1998).
[18] Y. Levin, Astrophys. J 517, 328 (1999).
[19] P. Arras et al., Astrophys. J 591, 1129 (2003).
[20] J. Heyl, Astrophys. J 574, L57 (2002).
[21] P.B. Jones, Astrophys. Lett. 5, 33 (1970). P.B. Jones,
Proc, Roy. Soc. (London) A323, 111 (1971). P.B. Jones,
Phys. Rev. Lett. 86, 1384 (2001). P.B. Jones, Phys. Rev.
D64, 084003 (2001).
[22] L. Lindblom and B. J. Owen, Phys. Rev. D65, 063006
(2002), astro-ph/0110558.
[23] P. Haensel, K. P. Levenfish, and D. G. Yakovlev, Astron.
and Astrophys. 381, 1080 (2002), astro-ph/0110575.
[24] M. Nayyar and B. J. Owen, Phys. Rev. D 73 (2006)
084001, astro-ph/0512041.
[25] N. Andersson, D. I. Jones, and K. D. Kokkotas, MNRAS
337, 1224 (2002).
[26] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys. Rev.
D70 (2004) 121501, gr-qc/0406085.
[27] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys.Rev.
D70 (2004) 124017, gr-qc/0409048.
[28] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys.Rev.
D71 (2005) 064029, gr-qc/0410072.
[29] A. K. Schenk, P. Arras, E. E. Flanagan, S. A. Teukol-
sky, I. Wasserman, Phys.Rev. D65 (2001) 024001,
gr-qc/0101092.
[30] L. Bildsten and G. Ushomirsky, Astrophys. J 529, L33
(2000).
[31] Y. Levin and G. Ushomirsky, MNRAS 322, 515 (2001).
[32] S. Yohida and U. Lee, Astrophys. J 546, 1121 (2001).
[33] K. Glampedakis and N. Andersson, astro-ph/0607105,
astro-ph/0411750.
[34] Y. S. Dimant, Phys. Rev. Lett. 84, 622 (2000).
[35] J. Wersinger, J. Finn, and E. Ott, Phys. Fluids 23, 1142
(1980).
[36] E. F. Brown, Ap. J 531, 988 (2000).
[37] H. Schatz, Phys. Rep. 294, 167 (1998).
[38] D. G. Yakovlev and K. P. Levenfish, Astron. Astrophys.
297, 717 (1995).
[39] D. G. Yakovlev and K. P. Levenfish, and Yu. A. Shibanov,
Soviet Phys.-Uspekhi, 42, 737 (1999).
[40] D. G. Yakovlev, A. D. Kaminker, and O. Y. Gnedin,
A&A, 379, L5 (2001).
[41] D. G. Yakovlev, and C. J. Pethick, Ann. Rev. Astron.
Astrophysics, 42, 169 (2004).
[42] G. Bryan, Philos. Trans. R. Soc. London A180, 187
(1889).
[43] C. Cutler and L. Lindblom, Astrophys. J 314, 234 (1987).
[44] K. D. Kokkotas and N. Stergioulas, Astron. and Astro-
phys. 341, 110 (1999).
[45] J. B. Kinney and G. Mendell, Phys.Rev. D67 024032
(2003).
[46] P. R. Brady, T. Creighton, C. Cutler, B. F. Schutz,
Phys. Rev. D 57, 2101 (1998), gr-qc/9702050. P. R.
Brady, T. Creighton, Phys. Rev. D 61, 082001 (2000),
gr-qc/9812014.
[47] B. J. Owen and L. Lindblom, Class.Quant.Grav. 19,
1247-1254 (2002), gr-qc/0111024.
[48] R. Prix, G. L. Comer, and N. Andersson, MNRAS 348,
625 (2004). N. Andersson and G. L. Comer, MNRAS
328,1129 (2001). N. Andersson, G. L. Comer and R. Prix,
MNRAS 354, 101 (2004).
[49] R. V. Wagoner, J. F. Hennawi, J. Liu, Proceedings of the
20th Texas Symposium on Relativistic Astrophysics, 781
(2001), astro-ph/0107229.
http://arxiv.org/abs/astro-ph/0110558
http://arxiv.org/abs/astro-ph/0110575
http://arxiv.org/abs/astro-ph/0512041
http://arxiv.org/abs/gr-qc/0406085
http://arxiv.org/abs/gr-qc/0409048
http://arxiv.org/abs/gr-qc/0410072
http://arxiv.org/abs/gr-qc/0101092
http://arxiv.org/abs/astro-ph/0607105
http://arxiv.org/abs/astro-ph/0411750
http://arxiv.org/abs/gr-qc/9702050
http://arxiv.org/abs/gr-qc/9812014
http://arxiv.org/abs/gr-qc/0111024
http://arxiv.org/abs/astro-ph/0107229

  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.

Spin Evolution of Accreting Neutron Stars: Nonlinear Development of the R-mode
Instability
Ruxandra Bondarescu, Saul A. Teukolsky, and Ira Wasserman
Center for Radiophysics and Space Research, Cornell University, Ithaca, NY 14853
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 that change slowly as the temperature and spin
of the star evolve. 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. The
evolution of the neutron star angular velocity and temperature follow easily calculated trajectories
along these sequences of quasi-stationary states. The outcome depends on whether or not the
star will reach thermal equilibrium, where the viscous heating by the three modes is equal to
the neutrino cooling (H = C curve). If, when the r-mode becomes unstable, the star spins at a
frequency below the maximum of the H = C curve, then it will reach a state of thermal equilibrium.
It can then either (1) undergo a cyclic evolution with a small cycle size 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 ≈ 106 years. If the star does not reach a state of thermal equilibrium, then a faster
thermal runaway (timescale of ≈ 100 years) occurs and the r-mode amplitude increases above the
second parametric instability threshold. Following this evolution requires more inertial modes to
be included. The sources of damping considered are shear viscosity, hyperon bulk viscosity and
viscosity within the core-crust boundary layer. We vary proprieties of the star such as the hyperon
superfluid transition temperature Tc, 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. We focus on Tc & 5 × 10
9 K where nonlinear effects are important. Wagoner [1] has
shown that a very low r-mode amplitude arises at smaller Tc. For all our bounded evolutions the
r-mode amplitude remains small ∼ 10−5. The spin frequency of accreting neutron stars is limited by
boundary layer viscosity to νmax ≈ 800Hz[Sns/(M1.4R6)]
4/11T
−2/11
. Fast rotators are allowed for
[Sns/(M1.4R6)]
4/11T
−2/11
∼ 1 and we find that in this case 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.
PACS numbers: 04.40.Dg, 04.30.Db, 97.10.Sj, 97.60.Jd
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 [2, 3],
which is driven by the gravitational radiation backreac-
tion force. Andersson [4] and Friedman and Morsink [5]
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 [6] that depends on
the angular velocity Ω and temperature T of the star.
Above a critical curve in the Ω−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 exponen-
tially, but soon it may enter a regime where other in-
ertial modes that couple to the r-mode become excited
and nonlinear effects become important. Roughly speak-
ing, nonlinear effects first become significant as the am-
plitude passes its first parametric instability threshold,
which is very low (∼ 10−5). Modeling and understand-
ing the nonlinear effects is crucial in determining (1) the
final saturation amplitude of the r-mode and (2) the lim-
iting 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 as-
sociated gravitational radiation could be detectable by
laser interferometers on Earth.
The r-mode instability has been proposed as an expla-
nation for the sub-breakup spin rates of both Low Mass
X-ray Binaries (LMXBs) [7, 8] and young, hot neutron
stars [6, 9]. The idea that gravitational radiation could
balance accretion was proposed independently by Bild-
sten [7] and Andersson et al. [8]. Cook, Shapiro and
Teukolsky [10, 11] model the recycling of pulsars to mil-
lisecond periods via accretion from a Keplerian disk onto
a bare neutron star with M = 1.4M⊙ when Ω = 0. De-
pending on the equation of state they found that spin
frequencies of between ≈ 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 millisec-
ond pulsars is 716 Hz for PSR J1748-2446ad [12, 13].
http://arxiv.org/abs/0704.0799v3
PSR B1937+21, which was discovered in 1982, was the
previous fastest known radio pulsar with a spin rate of
642 Hz [14]; 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 oscil-
lations, Chakrabarty et al. [15] claimed a cutoff limit of
νmax = 760 Hz (95% confidence); A more recent analy-
sis, which added two more pulsars to the sample, found
νmax = 730 Hz [16].
At first sight, one might conclude that mass shedding
or radial instability sets νmax, and that it is just above
the record ν = 716 Hz determined for PSR J1748-2446ad.
However, the nuclear equations of state consistent with
this picture all have rather large radii ≈ 16 − 17 km
for non-rotating 1.4 M⊙ models; see Table 1 in Cook et
al. [10]. For these equations of state, the r-mode insta-
bility should lead to ν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 de-
tection of the 716 Hz rotator is consistent with accretion
spin-up mitigated by the r-mode instability only for equa-
tions 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 estab-
lishing such constraints.
Since a physical model to follow the nonlinear phase of
the evolution was initially unavailable, Owen et al. [17]
proposed a simple one-mode evolution model in which
they assumed that nonlinear hydrodynamics effects satu-
rate 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 ampli-
tude (see Eqs. (3.16) and (3.17) in Ref. [17]). They used
this model to study the impact of the r-mode instability
on the spin evolution of young hot neutron stars assum-
ing 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 109 K
and spin down from a frequency close to the Kepler fre-
quency to about 100 Hz in a period of ∼ 1 yr [17].
Most subsequent investigations that did not perform
direct hydrodynamic simulations used the one-amplitude
model of Ref. [17] for studying the r-mode instability.
Levin [18] 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 ∼ 1; he
adopted a modified shear viscosity to match the maxi-
mum LMXB spin frequency of 330 Hz known in 1999.
Levin found that the neutron star followed a cyclic evo-
lution in the Ω − 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 (∼ 1). At this
point the spin and r-mode amplitude evolution equations
are changed, following the prescription of Ref. [17] to en-
sure 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 gravita-
tional radiation. However, following work by Arras et al.
[19] showing that nonlinear effects become significant at
small r-mode amplitude, Heyl [20] 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 α ≈ 10−5, and that
this would make some of the fastest spinning LMXBs in
our galaxy detectable by interferometers on Earth.
Jones [21] and Lindblom and Owen [22] 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 neu-
tron stars. While this additional high viscosity coefficient
could eliminate the instability altogether in newly born
neutron stars [21, 22, 23, 24], Nayyar and Owen [24] pro-
posed that it would enhance the probability of detection
of gravitational radiation from LMXBs by blocking the
thermal runaway.
The cyclic evolution found by Levin [18] and gener-
alized by Heyl [20] arises when shear or boundary layer
viscosity dominates the r-mode dissipation. In the evo-
lutionary picture of Nayyar and Owen [24], the r-mode
first becomes unstable at a temperature where shear
and boundary layer viscosity dominate, but the result-
ing 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 stabil-
ity curve at a point where the viscosity is an increas-
ing function of temperature. Such scenarios were stud-
ied by Wagoner [1] for hyperon bulk viscosty with low
hyperon superfluid transition temperature; similar evo-
lution was found for strange stars by Andersson, Jones
and Kokkotas [25]. 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. [26, 27, 28] for modes with n < 30, and hence in this
regime nonlinear effects may not play a role.
Schenk et al. [29] 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 nonlin-
ear perturbation theory via three-mode couplings. This
assumption was tested by Arras et al. [19] and Brink
et al. [26, 27, 28]. Arras et al. proposed that a turbu-
lent 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 in-
teraction of about 5000 modes via approximatively 1.3
million couplings of the 109 possible couplings among
the modes with n ≤ 30. The couplings were restricted to
mode triplets with a fractional detuning δω/(2Ω) < 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 thresh-
old 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 cal-
culation 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 dis-
sipates via boundary layer effects [30] and bulk viscos-
ity. 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 vary-
ing the free parameters of our model. In particular, we
vary: (1) the slippage factor Sns, which regulates the
boundary layer viscosity, between 0 and 1 (see for exam-
ple [31, 32, 33] 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 fdU, 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 Tc, which is believed
to be between 109− 1010 K (We use a single, effective Tc
rather than modelling its spatial variation.) We focus on
Tc & 5×109 K for which nonlinear effects are important.
For low Tc . 3 × 109 K, Wagoner [1] showed that the
evolution reaches a steady state at amplitudes below the
lowest parametric instability threshold found by Brink
et al. [28]. It is important to note that all our evolu-
tions 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. [28]. We evolve the coupled equa-
tions for the three-mode system numerically in conjunc-
tion 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 ob-
serve that the mode amplitudes settle into a series of
quasi-stationary states within a period of a few years af-
ter 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 equa-
tions (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 evo-
lution equations (six for the real and imaginary parts of
the mode amplitudes, which are complex [29]; one for the
spin, and one for the temperature) to two equations gov-
erning 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 equilib-
rium. This can be predicted by finding the curve where
the viscous heating by the three modes balances the neu-
trino cooling, referred to below as the Heating = Cooling
(H = C) curve. TheH = 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) un-
dergo 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
106 yr, and subsequently enter a regime of thermal run-
away. On the other hand, if the star cannot find a state
of thermal equilibrium, then it enters a regime of ther-
mogravitational runaway within a few hundred years of
crossing the r-mode stability curve. When this happens,
the r-mode amplitude increases beyond the second para-
metric 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 fre-
quencies and couplings are computed for an incompress-
ible 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 in-
volves 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 it-
self 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 vis-
cosity, 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 fea-
tures of our calculations ought to be robust. More sophis-
ticated 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 nu-
merous 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 viscos-
ity is relatively unimportant, soon after the star becomes
unstable its r-mode amplitude passes its first paramet-
ric 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 out-
line here (parametric instability, quasisteady evolution,
competition between heating and cooling) ought to ap-
ply 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 rota-
tional 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. IVA
discusses cyclic evolution in more detail. An evolution
that leads to an equilibrium steady state is presented
next in Sec. IVB. The two types of thermal runaway are
then discussed in Sec. IVC. 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
A. 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 τ for the time variable with
dτ = Ω dt Eq. (2.49) of Schenk et al. [29] can be rewritten
= iω̃αCα +
2iω̃ακ̃√
CβCγ , (1)
= iω̃βCβ −
2iω̃βκ̃√
= iω̃γCγ −
2iω̃γκ̃√
Here the scaled frequency ω̃j is defined to be ω̃j = ωj/Ω,
the dissipation rates of the daughter modes are γβ and γγ ,
γα is the sum of the driving and damping rates of the r-
mode γα = γGR−γαv, and the dimensionless coupling is
κ̃ = κ/(MR2Ω2). These amplitude variables are complex
and can be written in terms of the variables of Ref. [29]
as Cj(t) =
Ω(t)cj(t) (see Appendix A for a derivation
of Eqs. 1). The index j loops over the three modes j =
α, β, γ, where α labels the r-mode or parent mode and β
and γ 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 assump-
tion 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. [28]
for a full derivation). This amplitude is the parametric
instability threshold
|Cα|2 =
4κ̃2ω̃βω̃γΩ
1 + Ω2
γβ + γγ
, (2)
where the fractional detuning is δω̃ = ω̃α − ω̃β − ω̃γ .
Thorough explorations of the phase space of damped
three-mode systems were performed by Dimant [34] and
Wersinger et al. [35].
For the three modes at the lowest parametric instabil-
ity threshold, ω̃α ≈ 0.66, ω̃β ≈ 0.44, ω̃γ ≈ 0.22, κ̃ ≈ 0.19
and |δω̃| ≈ 3.82 × 10−6. Note that ω̃ is twice the w of
Brink et al. [26, 27, 28]. Here β labels the mode with
n = 13,m = −3 and γ labels the n = 14,m = 1 mode.
The amplitude the r-mode has to reach before exciting
these two daughter modes is |Cα| ≈ 1.5× 10−5
Ω [28].
We next rescale the rotational phase τ by the fractional
detuning as τ̃ = τ |δω̃| and the mode amplitudes by
|Cα|0 =
|δω̃|
ω̃βω̃γ
, |Cβ |0 =
|δω̃|
ω̃αω̃γ
, (3)
|Cγ |0 =
|δω̃|
ω̃βω̃α
which for the r-mode is, up to a factor of
Ω/Ωc,
the no-damping limit of the parametric instability thresh-
old below which no oscillations will occur. The coupled
equations become
|δω̃|
C̄α +
|δω̃|Ω̃
C̄α −
C̄βC̄γ , (4)
|δω̃|
C̄β −
|δω̃|Ω̃
C̄β −
C̄αC̄
|δω̃|
C̄γ −
|δω̃|Ω̃
C̄γ −
C̄αC̄
with C̄j = Cj/|Cj |0 and γ̃j = γj/Ωc being the newly
rescaled amplitudes and dissipation/driving rates, re-
spectively.
1. Quasi-Stationary Solution
In terms of amplitudes and phase variables Cj =
|Cj |eiφj Eqs. (4) can be rewritten as
d|C̄α|
Ω̃|δw̃|
|C̄α| −
sinφ|C̄β ||C̄γ |
, (5)
d|C̄β |
= − γ̃β
Ω̃|δw̃|
|C̄β |+
sinφ|C̄α||C̄γ |
d|C̄γ |
Ω̃|δw̃|
|C̄γ |+
sinφ|C̄α||C̄β |
|δω̃|
− cosφ
|C̄β ||C̄γ |
|C̄α|
− |C̄α||C̄γ |
|C̄β |
− |C̄β ||C̄α|
|C̄γ |
where we have defined the relative phase difference as
φ = φα − φβ − φγ . These equations have the stationary
solution
|C̄α|2 =
4γ̃βγ̃γ
Ω̃|δω̃|2
tan2 φ
, (6)
|C̄β |2 =
4γ̃αγ̃γ
Ω̃|δω̃|2
tan2 φ
|C̄γ |2 =
4γ̃αγ̃β
Ω̃|δω̃|2
tan2 φ
tanφ =
γ̃β + γ̃γ − γ̃α
Ω̃|δω̃|
Note that in the limit in which γβ+γγ >> γα the station-
ary solution for the r-mode amplitude |Cα| is the same
as the parametric instability threshold.
B. Temperature and Spin Evolution
The spin evolution equation is obtained from conser-
vation of total angular momentum J , where
J = IΩ + Jphys. (7)
Following Eq (K39-K42) of Schenk et al. [29] the physical
angular momentum of the perturbation can be written as
ΩJphys =
C⋆BCA
d3xρ[(Ω̂× ξ⋆B) · (Ω̂× ξA) (8)
− i (ω̃A + ω̃B)
ξ⋆B · (Ω̂× ξA)].
Since the eigenvectors ξA ∝ eimAφ the cross-terms will
vanish for modes with different magnetic quantum num-
bers m as
ei(mA−mB)φdφ = 0 for mA 6= mB. Eq. (8)
can be re-written for our triplet of modes as
Jphys = MR
2(kαα|Cα|2 + kββ|Cβ |2 + kγγ |Cγ |2), (9)
where kαα is defined as
kαα =
d3xρ[(Ω̂×ξ⋆α) · (Ω̂×ξα)− iω̃αξ⋆α · (Ω̂×ξα)]
and similarly for kββ and kγγ . In terms of the scaled vari-
ables C̄j = Cj/|Cj |0 (with |Cj |0 defined in Eq. (3)) the
angular momentum of the perturbation can be written
Jphys =
MR2Ωc|δω̃|2
(4k̃)2ω̃αω̃βω̃γ
(kαα|C̄α|2ω̃α (11)
+kββ|C̄β |2ω̃β + kγγω̃γ |C̄γ |2).
We chose the same normalization for the eigenfuctions
as Refs. [19, 26, 27, 28, 29] so that at unit amplitude all
modes have the same energy ǫα = MR
2Ω2. The energy
of a mode α is Eα = MR
2Ω2|cα|2 = MR2Ω|Cα|2. The
rotating frame energy is the same as the canonical energy
and physical energy [29]. The canonical angular momen-
tum and the canonical energy of the perturbation satisfy
the general relation Ec = −(ω/m)Jc [3].
Angular momentum is gained because of accretion and
lost via gravitational waves emission
= 2γGRJc rmode + Ṁ
GMR, (12)
where Jc rmode = −(mα/ωα)ǫα|cα|2 = −3MR2Ω|cα|2 =
−3MR2|Cα|2. Eq. (12) can be rewritten in terms of the
scaled variables C̄j as
= −6γ̃GR
MR2Ωc|δω̃|
(4k̃)2ω̃βω̃γ
|C̄α|2 +
ΩcΩ̃|δω̃|
. (13)
Thermal energy conservation gives the temperature evo-
lution equation
C(T )
2Ejγj +KnṀc
2 − Lν(T ), (14)
= 2MR2Ω(γα v|Cα|2 + γβ |Cβ |2
+ γγ |Cγ |2) +KnṀc2 − Lν(T ).
The three terms on the right hand side of the equa-
tion represent viscous heating, nuclear heating and neu-
trino cooling. The specific heat is taken to be C(T ) ≈
1.5 × 1038 T8 erg K−1, where T = T8 × 108 K. Nu-
clear heating occurs because of pycnonuclear reactions
and neutron emission in the inner crust [36]. At large
accretion rates such as that of the brightest LMXBs of
Ṁ ≈ 10−8M⊙/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 neu-
trino emission is assumed to regulate the temperature of
the star [36, 37]. The nuclear heating constant is taken
to be Kn ≈ 1×10−3 [36]. Following Ref. [1], we take the
neutrino luminosity to be
Lν = LdUT
8RdU(T/Tp) + LmUT
8RmU(T/Tp) (15)
+ Le−iT
8 + Ln−nT
8 + LCpT
where the constants for the modified and direct URCA re-
actions are defined by LmU = 1.0×1032 erg sec−1, LdU =
fdU × 108LmU [38, 39], and the electron-ion, neutron-
neutron neutrino bremsstrahlung and Cooper pairing of
neutrons are given by Le−i = 9.1 × 1029 erg sec−1 [36],
Ln−n ≈ 0.01LmU, LCp = 8.9 × 1031 erg sec−1 [40]. The
fraction of the star fdU that is above the density thresh-
old for direct URCA reactions is in general dependent on
the equation of state [41] and in this work we treat fdU
a free parameter with values between 0 and 1.
The proton superfluid reduction factors for the modi-
fied and direct URCA reactions are taken from Ref. [39]
(see Eqs. (32) and (51) in Ref. [39]):
RdU(T/Tp) =
0.2312 +
(0.76880)2 + (0.1438v)2
× exp
3.427−
(3.427)2 + v2
RmU(T/Tp) =
0.2414 +
(0.7586)2 + (0.1318v)2
× exp
5.339−
(5.339)2 + (2v)2
where the dimensionless gap amplitude v for the singlet
type superfluidity is given by
1.456− 0.157
+ 1.764
. (17)
Similar to Ref. [1], we use Tp = 5.0× 109 K. In terms of
the scaled variables Eq. (14) becomes
C(T )
2MR2Ω2c |δω̃|
(4κ̃)2ω̃αω̃βω̃γ
(ω̃αγ̃α v|C̄α|2 + ω̃βγ̃β |C̄β |2(18)
+ω̃γγ̃γ |C̄γ |2) +
KnṀc
2 − Lν(T )
ΩcΩ̃|δω̃|
C. 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.,
dCi/dτ̃ ≈ 0, the spin and thermal evolution equations
can be rewritten by taking J ≈ IΩ and using Eqs. (6) in
Eq. (13).
= − 6γ̃GR
Ω̃2|δω̃|
γ̃β γ̃γ
4k̃2Ĩω̃βω̃γ
tan2 φ
MR2ĨΩ̃|δω̃|
where Ĩ = I/(MR2). The thermal evolution of the sys-
tem is given by
C(T )
2MR2Ω2c
(4κ̃)2ω̃αω̃βω̃γ
γ̃αγ̃β γ̃γ
Ω̃|δω̃|
ω̃αγ̃α,v
+ ω̃β (20)
+ω̃γ)
tan2 φ
KnṀc
2 − Lν(T )
ΩcΩ̃|δω̃|
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 evolu-
tion very well throughout the unstable regime. These
equations are a minimal physical model for the effects of
nonlinear coupling on r-mode evolution.
D. Sources of Driving and Dissipation
The damping mechanisms are shear viscosity, bound-
ary layer viscosity and hyperon bulk viscosity; for modes
j = α, β, γ we write
γj v(Ω, T ) = γj sh(T ) + γj bl(Ω, T ) + γj hb(Ω, T ). (21)
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 β and
n = 14,m = 1 labeled as γ) is affected only by the vis-
cous damping. Brink et al. [26, 27, 28] determined that
this pair of modes is excited at the lowest parametric
instability threshold. Their model uses the Bryan [42]
modes of an incompressible star, which has the advan-
tage 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 (Γ1 → ∞ in
Eq. (29))
For our benchmark calculations, we adopt the neutron
star model of Owen et al. Ref. [6] (n = 1 polytrope,
M = 1.4M⊙, Ωc = 8.4 × 103 rad sec−1 and R = 12.53
km) and use their gravitational driving rate and shear
viscous damping rate for the r-mode
γGR(Ω) ≃
sec−1, (22)
γα sh(T ) ≃
where τsh = 2.56 × 106 sec. (In Sec. V we consider ap-
proximate 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
γβ sh(T ) ≃ 3.48× 10−4 sec−1
, (23)
γγ sh(T ) ≃ 4.52× 10−4 sec−1
The geometric contribution γsh/η of the individual modes
increases significantly with the degree n of the mode scal-
ing approximatively like n3 for large n (see Eq. (29) of
Brink et al. [27] 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 magni-
tude larger than that of the r-mode.
The damping due to boundary layer viscosity is calcu-
lated using Eq. (4) of Ref. [30],
γα bl(T,Ω) ≃ 0.009 sec−1 S2ns
, (24)
γβ bl(T,Ω) ≃ 0.028 sec−1 S2ns
γγ bl(T,Ω) ≃ 0.021 sec−1S2ns
Analogous to Wagoner [1], we allow the slippage fac-
tor Sns to vary. The slippage factor is defined by Refs.
[1, 31, 45] to be S2ns = (2S
n + S
s )/3, with Sn being the
fractional difference in velocity of the normal fluid be-
tween the crust and the core [31] and Ss the fractional
degree of pinning of the vortices in the crust [45]. Note
that γβ bl and γγ bl are both greater than 2 × γα bl and
can easily be comparable to γ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 [24].
This rate is taken to have a form similar to that taken
by Wagoner [1]
γα hb = fhb
t−20α τ(T )Ω̃
1 + (ω̃αΩτ(T ))2
, (25)
and for the daughter modes
γβ hb = fhb
t−20β τ(T )ω̃
1 + (ω̃βΩτ(T ))2
, (26)
and similarly foe γγ hb. The relaxation timescale
τ(T ) =
Rhb(T/Tc)
The reduction factor is taken to be the product of two
single-particle reduction factors [23, 24]
Rhb single(T/Tc) =
a5/4 + b1/2
0.5068−
0.50682 + y2
where a = 1 + 0.3118y2, b = 1 + 2.556y2 and y =
1.0− T/Tc(1.456 − 0.157
Tc/T + 1.764Tc/T ). The
constants t1 ≈ 10−4 sec and t0α ≈ 0.00058 sec are found
by fitting the results of Ref. [24]. The factor fhb allows
for physical uncertainties; we take fhb = 1 throughout
the body of the paper since Tc , which enters γj hb ex-
ponentially, 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 Γ−21 by
− ĖB j =
ξj · ∇p
. (29)
This approximation was proposed by Cutler and Lind-
blom [43] and adopted by Kokkotas and Stergioulas [44]
for the r-mode and by Brink et al. [27] for the inertial
modes. The adiabatic index Γ1 is regarded as a parame-
ter; we use Γ1 ≈ 2. The damping rate is
γj hb = −
ĖB j
, (30)
where ǫ = MR2Ω2 is the mode’s energy in the rotat-
ing frame at unit amplitude and j = β, γ. Using this
procedure, we calculate
t0β ≈ 1.4× 10−5 sec, (31)
t0γ ≈ 1.0× 10−5 sec.
III. SUMMARY OF RESULTS
Fig. 1(a) shows possible evolutionary trajectories of a
neutron star in the angular velocity-temperature Ω̃− T8
plane, where T = T8×108 K is the core temperature, and
Ω̃ = Ω/Ωc = Ω/
πGρ̄ with ρ̄ the mean density of the
neutron star. Fig. 1(b) displays the regions in fdU − Sns
in which the trajectories occur. Here fdU represents the
fraction of the star that is above the density threshold
for direct URCA reactions and Sns is the slippage factor
that reduces the relative motion between the crust and
the core taking into account the elasticity of the crust
[31]. The stability regions are shown at fixed hyperon
superfluidity temperature, Tc = 5.0 × 109 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
T8 = T8 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 Ω because of viscous dissipation
in all three modes.
Phase 3. Several trajectories are possible depending on
FIG. 1: (a)Typical trajectories for the four observed evolu-
tion scenarios are shown in the Ω̃ - T8 phase space, where
Ω̃ = Ω/Ωc. The dashed lines (H = C curves) represent
the points in the Ω̃ − T8 phase space where the dissipative
effects of the heating from the three-modes exactly compen-
sate the neutrino cooling for the given set of parameters (Sns,
fdU, Tc, ...) of each evolution. (b)The corresponding sta-
bility regions for which these scenarios occur are plotted at
fixed hyperon superfluidity temperature Tc = 5.0 × 10
while varying fdU and Sns. The position of the initial angu-
lar velocity and temperature (Ω̃in, T8 in) with respect to the
maximum of this curve determines the stability of the evo-
lution. (I) Ω̃in > Ω̃H=C max. Trajectory R1. Fast Runaway
Region. After the r-mode becomes unstable the star heats
up, does not find a thermal equilibrium state and continues
heating up until a thermogravitational runaway occurs. (II)
Ω̃in < Ω̃H=C max. The evolutions are either stable or, if there
is a runaway, it occurs on timescales comparable to the ac-
cretion timescale. The possible trajectories are (1)Trajectory
C. Cycle Region. (2) Trajectories S1 and S2. Steady State
Region. (3) Trajectory R2. Slow Runaway Region.
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 [24]. 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 T8 in at constant Ω̃ and
then the cycle repeats itself. At Tc = 5.0 × 109 K this
scenario occurs for values of Sns < 0.50 and large enough
values of fdU. However, if Tc is larger, the cycle region
in the fdU-Sns phase space increases dramatically (see
Fig. 9(a)). Note that our cycles are different from those
obtained by Levin [18] in that the spin-down phase
does not start when the r-mode amplitude saturates
(or in our case when it reaches the parametric insta-
bility 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 ∼ 105 and so the part
of the cycle in which the r-mode is unstable also lasts
longer than in Ref. [18]. Also, our cycles are narrow.
During spin-down the temperature changes by less than
20 % and Ω̃ changes by less than 10% of the initial value.
(See Sec. 2 for a detailed example.)
(ii) Steady State. For small Sns and large enough
fdU (fdU & 5 × 10−5, Sns . 0.04; see Fig. 1(b)) the
star evolves towards an Ω̃ 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 Sns and very small fdU
(Sns . 0.03, fdU < 5×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 Tc ∼ 1010, or for models with no hyperons at all,
there would be no runaway region (See Fig. 9(a) for an
fdU − Sns scenario space with a larger Tc = 6.5× 109 K
where the fast runaway region has shrunk dramatically
FIG. 2: Two cyclic trajectories in the Ω̃ − T8 plane are dis-
played for a star with Tc = 5.0 × 10
9 K and (a) fdU = 0.15
and Sns = 0.10, and (b) fdU = 0.142 and Sns = 0.35, which is
close to the border between the stable and unstable region (see
Fig. 1(b)). The thick solid line labeled as the Heating = Cool-
ing (H = C) curve is the locus of points in this phase space
where the neutrino cooling is equal to the viscous heating due
to the unstable modes. The other solid line representing the
r-mode stability curve is defined by setting the gravitational
driving rate equal to the viscous damping rate. The part of
the curve that decreases with T8 is dominated by boundary
layer and shear viscosity, while the part of the curve that has
a positive slope is dominated by hyperon bulk viscosity. In
portion a1 → b1 of the trajectory the star heats up at con-
stant Ω̃. Part b1 → c1 represents the spin down stage, which
occurs when the viscous heating is equal to the neutrino cool-
ing. c1 → d1 shows the star cooling back to the initial T8.
Segment d1 → a1 displays the accretional spin-up of the star
back to the r-mode stability curve. The cycle a2 → d2 pro-
ceeds in the same way. This cycle is close to the peak of the
H = C curve. Configurations above this peak will run away.
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 Ṁ =
10−8M⊙/yr and Tc = 5.0× 109 K.
A. 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) Sns = 0.10 and fdU = 0.15 and (C2) Sns = 0.35
and fdU = 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 T8 in ≈ 3.29 for both cases. The initial
Ω is the angular velocity that corresponds to this tem-
perature on the r-mode stability curve, which differs for
the different Sns (Ω̃in = 0.183 for C1 and Ω̃in = 0.288 for
Figs. 2(a) and (b) display the cyclic evolution for tra-
jectories C1 and C2 of Fig. 1(b). In leg a1 → b1 of the
trajectory the r-mode and, once the r-mode amplitude
increases above the first parametric instability thresh-
old, the two daughter modes it excites, viscously heat
up the star until point b1 when the neutrino cooling bal-
ances the viscous dissipation. This part of the evolu-
tion occurs at constant angular velocity over a period
of theat−up ≈ 100 yr and a total temperature change
(∆T )a1−b1 ≈ 0.80 (≈ 24% of T8 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 b1 → c1 of the
trajectory as it spins down and cools down back to the r-
mode stability curve. This spin-down stage lasts a time
tspin−downb1−c1 ≈ 23, 000 yr that is much longer than
the heat-up period. This timescale is very sensitive to
changes in the slippage factor and can reach 106 yr for
smaller values of Sns 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%, (∆Ω̃)b1−c1 ≈ 0.0066. The temperature also changes
by only about 2%, (∆T8)b1−c1 ≈ 0.08 in this spin-down
period. Segment c1 → d1 represents the cooling of the
star to the initial temperature on a timescale of ∼ 2, 000
yr. In part d1 → a1 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 tra-
jectory is the longest-lasting one, taking ≈ 200, 000 yr
at our chosen Ṁ of 10−8M⊙yr
−1. The cycle C2 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 there-
fore close to the peak of the H = C curve. Configura-
tions above this peak (e.g., with the same fdU and higher
Sns) 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 un-
stable and initially grows exponentially. Once it has in-
creased 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 exponen-
tially while the r-mode decreases. Similarly, when daugh-
ter 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 ampli-
tude 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 thresh-
old 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.
B. 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 re-
stricted to stars with small slippage factor (Sns . 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 sta-
bility curve the gravitational driving rate is greater than
the viscous damping rate and the r-mode grows exponen-
tially 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 evo-
lution and after a few oscillations the three modes settle
into their quasi-stationary states, which change only sec-
ularly as the spin and temperature of the star evolve.
The modes heat the star viscously at constant Ω̃ in seg-
ment a → b of the trajectory for theat−up ≈ 1, 100 yr.
At point b, the star reaches a state of thermal balance.
In leg b → c the star continues its evolution in thermal
equilibrium and slowly spins up due to accretion until
the angular velocity evolution also reaches an equilib-
5.12 5.13
τ [x 10
0.422
0.424
0.426
0.428
0.432
|cα Th|
0.4 0.5 0.6 0.7 0.8 0.9 1
τ [x 10
0 0.5 1 1.5 2 2.5 3 3.5 4
τ [x 10
0.425
0.435
FIG. 3: (a)The amplitudes of the r-mode |Cα| and of the
n = 13, m = −3 and n = 14, m = 1 inertial modes |Cβ | and
|Cγ | are shown as a function of time for a star that executes
a cyclic evolution (same parameters as in Fig. 2). The low-
est parametric instability threshold is also displayed. (b)The
ratio of the r-mode amplitude to the parametric instability
threshold is plotted as a function of time. It can be seen that
once the r-mode crosses the parametric instability threshold
it remains close to it for the rest of the evolution. (c)The
parametric instability threshold is displayed as a function of
time. Its value changes as the angular velocity and tempera-
ture evolve.
FIG. 4: The trajectory of a neutron star in the Ω̃− T8 phase
space is shown for a model with Tc = 5.0× 10
9 K, fdU = 0.03
and Sns = 0.03 that reaches an equilibrium steady state. The
star spins up until it crosses the r-mode stability curve and
the r-mode becomes unstable. The r-mode then quickly grows
to the first parametric instability threshold and excites the
daughter modes. In leg a → b of the trajectory the star is
viscously heated by the mode triplet until the system reaches
thermal equilibrium. Segment b → c shows the star contin-
uing to heat and spin up in thermal equilibrium until the
accretion torque is balanced by the gravitational radiation
emission. The r-mode stability curve represents the points
in phase space where the viscous driving rate is equal to the
gravitational driving rate. The H=C curve is the locus of
points where the viscous dissipation due to the mode triplet
balances the neutrino cooling.
FIG. 5: The (Ω̃, T8) initial values (region delimited by the
solid line) that lead to equilibrium steady states and their
corresponding final steady state values (region enclosed by
the dashed line) are shown. Since both the initial and final
values of T8 are low, these evolutions are roughly independent
of Tc.
rium. The timescale to reach an equilibrium steady state
is tsteady ≈ 3.5× 106 yr for this set of parameters.
Fig. 5 displays the possible initial values for the angu-
lar velocity Ω̃ and temperature T8 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 Ω also tend to be small. Thus, although an interest-
ing physical regime, this case is most likely not relevant
to recycling by accretion to create pulsars with spin fre-
quencies as large as 716 Hz. Note that a steady state
can be achieved when Sns = 0. This is the probable end
state of the problem first calculated by Levin [18]. The
reason we do not find a cycle at low Sns is twofold: (1)
the shear viscosity we are using is lower (shear viscosity
in Ref. [18] is amplified by a factor of 244), and (2) the
nonlinear couplings keep all mode amplitudes small.
C. 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 run-
away 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 contin-
ues to spin up by accretion until it climbs to the peak of
the H = C curve, thermal equilibrium fails and a run-
away occurs.
1. Fast Runaway
A typical trajectory of a star that goes through a rapid
thermal runaway is displayed in Fig. 6. This star has
Sns = 0.25 and fdU = 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 (Ω,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. [24].
Fig. 6(a) and (b) compare both the temperature evolu-
tion and the trajectory in the Ω̃−T8 plane of the star for a
0 0.5 1 1.5 2 2.5 3 3.5
τ [x 10
Full Amplitudes-Ω-T Evolution
Steady State Evolution
2nd Parametric Instability Threshold
3 3.5 4 4.5 5 5.5 6 6.5
r-mode Stability Curve
Full Amplitude-Ω-T Evolution
Steady State Evolution
2nd Parametric Instability Threshold
FIG. 6: This plot compares the full evolution resulting from
solving Eqs. (4),(13),(18) with the reduced Ω − T evolution
that assumes the amplitudes go through a series of steady
states Eqs. (19)-(20) for a model with Tc = 5.0 × 10
fdU = 0.058 and Sns = 0.25. (a) The temperature is dis-
played as a function of time for the two different methods.
(b) The angular velocity Ω̃ = Ω/Ωc is shown as a function
of temperature. The evolution occurs at constant spin fre-
quency. It can be seen that the steady-state amplitude ap-
proximation is extremely good. The ‘X’ shows the point at
which the r-mode crosses its second lowest parametric insta-
bility threshold, where additional dissipation would become
operative.
simulation solving the full set of equations to a simulation
that assumes quasi-stationary solutions for the three am-
plitudes and evolves only the angular velocity and tem-
perature 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 re-
duced 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 [24]. However, the r-mode crosses its second low-
FIG. 7: The trajectory of a neutron star in the Ω̃ − T8
phase space is shown for a model with Tc = 5.0 × 10
fdU = 4.0 × 10
−5 and Sns = 0.02 that goes through a slow
thermogravitational runaway. Portion a → b of the trajectory
shows the mode triplet heating up the neutron star through
boundary layer and shear viscosity until the system reaches
thermal equilibrium. Segment b → c represents the accre-
tional spin-up of the star in thermal equilibrium. The dotted-
dashed line is the locus of points where the viscous dissipation
of the mode triplet is equal to the neutrino cooling, and is la-
beled as the H = C curve. The star reaches the maximum of
this curve and fails to reach an equilibrium between the ac-
cretion torque and gravitational emission. It then continues
heating at constant angular velocity and crosses its second
lowest parametric instability threshold, at which point more
modes would need to be included to make the evolution accu-
rate. Eventually the star reaches the r-mode stability curve
again.
est 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.
2. Slow Runaway
In this section we examine evolutions in which the neu-
tron star has both a very small slippage factor, Sns .
0.03, and only a small percentage of the star is above the
threshold for direct URCA reactions, fdU < 5 × 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 insta-
bility 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 Ω in part a → b of the trajec-
tory on a timescale of about 5, 000 yr. At point b, it
FIG. 8: The spin-down timescale is shown as slippage fac-
tor Sns and fraction of the star subject to direct URCA fdU
for cyclic evolutions are varied for a fixed hyperon critical
temperature of Tc = 5.0 × 10
9 K. This timescale dominates
the heat-up timescale and hence represents the time the star
spends above the r-mode instability curve. It increases as the
viscosity is lowered and the star gets closer to the steady state
region.
reaches thermal equilibrium. In leg b → c of the tra-
jectory, the star continues its evolution by ascending the
H = C curve and spinning up because of accretion for
about 2 × 106 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 expo-
nentially 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 instabil-
ity 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 Sns and fdU are var-
ied. For large enough values of Sns the boundary layer
viscosity dominates. In this region of phase space the
spin-down timescale can be approximated by
tspin−down =
dΩ̃ (32)
FIG. 9: (a)The stability regions are plotted at fixed hyperon
superfluidity temperature Tc = 6.5×10
9 K, while varying fdU
and Sns. The steady state region remains roughly the same
as in Fig. 1(b), the slow run-away region disappears, and the
cycle region increases dramatically while shrinking the fast-
runaway region. (b) The spin-down timescale is shown for the
cyclic evolutions in part (a).
Ĩτ0GR
(4κ̃)2ω̃βω̃γ
|δω̃|2
|C̄α|2
<Ω̃>6
≈ 250 yr ∆νkHz
<νkHz>7
M1.4R
|cthα |
where M1.4 = M/(1.4M⊙), R6 = R/(10
6cm),
νkHz = ν/1kHz, Ĩ = 0.261 [17], the r-mode am-
plitude at its parametric instability threshold |cthα | ≈
|δω̃|/(4κ̃
ω̃βω̃γ) ≈ 1.5×10−5, and C̄α =
Ω̃|cα|/|cα|th.
This approximation agrees with spin-up timescales ob-
tained from our simulations to ∼ 25%.
The maximum ν is approximately the same as the ini-
tial 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
νmax ≈ 800Hz
M1.4R6
)4/11
. (33)
Thus, the spin-down timescale is very sensitive to
the slippage factor tspin−down ∝ S−24/11ns (∆νkHz/νkHz).
The dependences on fdU and accretion rate Ṁ are
much weaker; a rough approximation, obtained by
matching direct URCA cooling and nuclear heat-
ing, is T8 in ∝ Ṁ1/6f−1/6dU R
1.4 , and νmax ∝
dU Ṁ
−1/33R
−34/99
1.4 . The gravitational
wave amplitude measured at distance d [46, 47] is
h ≈ 1.6R
τ0GRc
Ω̃3|cα| (34)
≈ 3× 10−25
10kpc
M1.4R
Taking ν ≈ νmax gives
h ∝ S12/11ns M
−1/33
1.4 R
dU Ṁ
−1/11. (35)
The maximum distance at which sources could be de-
tected by Advanced LIGO interferometers, assuming
hmin = 10
−27, [46] is
dmax ≈ 3Mpc
10−27
M1.4R
|cthα |
≈ 1.5Mpc
10−27
S12/11ns M
−1/11
1.4 R
21/11
× T−6/118
|cthα |
Eqs. (33) and (36) 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
ν & 700 Hz, we must require (Sns/M1.4R6
T8in)
4/11 & 1.
Assuming this to be true, dmax . 1-1.5 Mpc. However,
tspin−down ≈ 1000 yr at most, making detection unlikely
for any given source. Moreover, unless Sns can differ sub-
stantially from one neutron star to another, only those
with ν 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 Sns . 0.50, assuming Tc ≈ 5.0× 109
K; Eq. (33) then implies ν . 450 Hz. This would restrict
detectable gravitational radiation to galactic sources, al-
though the duration of the unstable phase could be
longer.
Within the context of our three mode calculation,
Sns > 0.50, which is needed for explaining the fastest pul-
sars, would imply fast runaway. There are two possible
resolutions to this problem. One is that including addi-
tional modes prevents the runaway; we shall investigate
this in subsequent papers. The second is that Tc is larger,
or that neutron stars do not contain hyperons (e.g., be-
cause they are insufficiently dense). Fig. 9(a) shows the
same phase plane as Fig. 1(b) but with Tc = 6.5×109 K,
and Fig. 9(b) shows the results for tspin−down analogous
to Fig. 8. Larger Tc permits spin cycles for higher values
of Sns (and hence ν), but the time spent in the unstable
regime is shorter.
VI. CONCLUSIONS
In this paper, we model the nonlinear saturation of un-
stable 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 en-
ergy 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 slip-
page 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 evolu-
tions we find that the mode amplitudes quickly settle
into a series of quasi-stationary states that can be calcu-
lated algebraically, and depend weakly on angular veloc-
ity 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 evo-
lution 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 pos-
sible 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 Ω 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 Sns ∼ 1) and 106 yr (small Sns ∼ 0.05). Our
cycles are different from those previously found by Ref.
[18] in that our amplitudes remain small, ∼ 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 ∼ 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 thermo-
gravitational runaway on a timescale of ∼ 106 yr. These
scenarios occur for very low viscosity (Sns . 0.04). Al-
though theoretically interesting, they do not allow for
very fast rotators of ∼ 700 Hz.
Alternatively, if the star does not reach thermal equi-
librium, 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 [24]. 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 Tc & 5× 109 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 ro-
tation requires large dissipation, as has long been recog-
nized [18, 30] and these models can only achieve ν & 700
Hz if boundary layer viscosity is very large. Alterna-
tively, at lower Tc . 3 × 109 K, large rotation rates can
be achieved at r-mode amplitudes below the first para-
metric instability threshold [1]. Nayyar and Owen found
that increasing the mass of the star for the same equation
of state makes the hyperon bulk viscosity become impor-
tant at lower temperatures [24]. Conceivably, there are
accreting neutron stars with relatively low masses that
have lower central densities and small hyperon popula-
tions. These could evolve as detailed here and only spin
up to modest frequencies. Hyperons could be more im-
portant in more massive neutron stars leading to larger
spin rates and very small steady state r-mode amplitude
as found by Wagoner [1].
Our models imply small r-mode amplitudes of ∼ 10−5
and therefore gravitational radiation detectable by ad-
vanced 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 νmax ≈
800Hz[Sns/(M1.4R6)]
4/11T
−2/11
8 . In order to allow for
fast rotators of & 700 Hz in our models a large bound-
ary layer viscosity with (Sns/M1.4R6
T8in)
4/11 ∼ 1 is
required. Slippage factors of order ∼ 1 lead to time peri-
ods 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 fre-
quencies. We also note that in this model we have con-
sidered only very fast accretors with Ṁ ∼ 10−8M⊙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 [48], superfluid-
ity, 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 quanti-
tatively in a model different from ours in detail, we ex-
pect the qualitative behaviors we have found to be ro-
bust, as they are well described by quasi-stationary mode
evolutions whose slow variations are determined by com-
petitions between dissipation and neutrino cooling, and
accretion spin-up and gravitational radiation spin-down.
In our model, it seems that three mode evolution involv-
ing 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. [26, 27, 28]
that includes accretion spin-up, viscous heating and neu-
trino 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 Éanna
Flanagan for useful discussions. RB would especially
like to thank Jeandrew for useful discussions, encourage-
ment and advice at the beginning of this project, with-
out which the project would not have been started. RB
is very grateful to Gregory Daues for steady encourage-
ment 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 Ω.
The Lagrangian density as given by Eq. (A1) in Schenk
et al. [29] is
L = 1
ξ̇ · ξ̇ + 1
ξ̇ ·B · ξ − 1
ξ ·C · ξ + aext(t) · ξ, (A-1)
where the operators B · ξ = 2Ω× ξ and
ρ(C · ξ)i = −∇i(Γ1p∇jξj) +∇ip∇jξj + ρ∇iδφ (A-2)
− ∇jp∇iξj + ρξj∇j∇iφ+ ρξj∇j∇iφrot
with φrot = −(1/2)(Ω × 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
, dτ = Ωdt. (A-3)
In terms of these new variables the Lagrangian density
can be written as
L̃ = 1
ξ̃′ · ξ̃′ + 1
ξ̃′ · (B̃ · ξ̃) + (
|ξ̃|2 (A-4)
ξ̃ · C̃ · ξ̃ + aext(t)
· ξ̃,
where the primes denote derivatives with respect to τ ,
B̃ = Ω−1B and C̃ = Ω−2C. The momentum canonically
conjugate to ξ̃ is
= ξ̃′ + Ω̂× ξ̃. (A-5)
The associated Hamiltonian density is
B̃ · ξ̃
|ξ̃|2 +
ξ̃ · C̃ · ξ̃ −
· ξ̃.
(A-6)
Hamilton’s equations of motions can be written as
ζ̃′ = T · ζ̃ + F(τ), (A-7)
where
the operator T is T = T0 + T1 with
B̃2 − C̃ − 1
Ω)′′√
F(τ) =
We assume solutions of the form ζ̃(τ,x) = eiω̃tζ̃(x). Spe-
cializing to the case of no forcing term aext = 0 leads to
the eigenvalue equation
(T0 − iω̃)ζ̃(x) = 0. (A-8)
Since the operator T0 is not Hermitian it will have dis-
tinct right and left eigenvectors. Similar to Schenk et
al. [29] we label the right eigenvectors of T as ζ̃A, and
the associated eigenfrequencies as ω̃A = ωA/Ω, and the
eigenvalue equation above becomes
(T0 − iω̃A)ζ̃A(x) = 0. (A-9)
The left eigenvectors χA satisfy
0 − iω̃⋆A)χ̃A = 0, (A-10)
where
B̃2 − C̃
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 be-
tween right and left eigenvectors is
χ̃A, ζ̃B
d3xρ(x)χ̃
A · ζ̃B = δAB. (A-11)
We can expand ζ(τ,x) in this basis as
ζ(τ,x) =
CA(τ)ζA(x), (A-12)
where the coefficients CA are given by the inverse of this
mode expansion
CA(τ) =
χ̃A, ζ̃(τ,x)
. (A-13)
Using Eqs. (B-2,A-9,A-11) in Eq. (A-7) leads to the equa-
tions of motion for the mode amplitudes
C′A − iω̃ACA = g(τ)
(A-14)
+ 〈χ̃A, F (τ)〉 ,
where g(τ) = (
Ω)′′/
Ω. Following Sec. IV of Schenk
et al. [29] we replace the externally applied acceleration
by the nonlinear acceleration given by Eq. (4.2) of Ref.
[29]. The inner product can be written in terms of the
displacement variable ξ̃. The left eigenvectors are
χ̃A =
where τ̃A can be chosen to be proportional to ξ̃A because
they satisfy the same matrix equation.
τ̃A = −iξ̃A/b̃A, (A-15)
which corresponds to Eq. (A-45) in Schenk et al. [29]
with the proportionality constant b̃A = Ω
−1bA =
MR2/ω̃A (also given by Eq. (2.36) of Ref. [29]).
The equations of motion for the mode amplitudes be-
C′A − iω̃ACA =
ig(τ)
d3xξ̃⋆A · ξ̃B(A-16)
κ̃⋆ABCC
where the nonlinear coupling κ̃ABC = κABC/(MR
and κABC is explicitly give by Eq. (4.20) of Ref. [29]. The
g(τ) integral mixes only modes with mA = mB because
of the eimφ dependence of the displacement eigenvectors
ξ̃. (
dφei(mA−mB)φ = 0 if mA 6= mB.) 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.
(A-16) 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. [49], each
of the five variables is expanded about its equilibrium
(Xj)e as follows
Xj(τ̃ ) = {|C̄α|, |C̄β |, |C̄γ |, φ, T8} = (Xj)e[1 + ζj(τ̃ )]
(B-1)
where the perturbation |ζj | << 1 and j = α, β, γ, T . The
expansion leads to a first order differential equation for
each ζj
(γ̃α)e
Ω̃|δω̃|
ζα − ζβ − ζγ −
ζφ (B-2)
(γ̃β)e
Ω̃|δω̃|
ζα − ζβ + ζγ +
(γ̃γ)e
Ω̃|δω̃|
ζα + ζβ − ζγ +
φe tanφe
γ̃α + γ̃β + γ̃γ
Ω̃|δω̃|
−γ̃α − γ̃β + γ̃γ
Ω̃|δω̃|
−γ̃α + γ̃β − γ̃γ
Ω̃|δω̃|
(γ̃α − γ̃β − γ̃γ)e
Ω̃|δω̃|
MR2Ω2c γ̃αγ̃β γ̃γ
2κ̃2ω̃αω̃βω̃γΩ̃|δω̃|C(Te)T8e
tanφ2e
γ̃α v
ζα + ω̃βζβ + ω̃γζγ
+ T8e
+ ω̃β
+ ω̃γ
ΩcΩ̃|δω̃|C(Te)
where the equilibrium amplitudes |Cj |e have been written
in terms of the corresponding driving and damping rates
using Eqs. (6). Eq. (B-2) can be written in matrix form
= Aijζi. (B-3)
Let ζj ∝ exp(λτ̃ ). The determinant ||Aij − λδij || = 0
leads to the eigenvalue equation
λ5 + a4λ
4 + a3λ
3 + a2λ
2 + a1λ+ a0 = 0. (B-4)
The coefficients aj with j = 0, 4 are
a4 = 2 tanφe =
γ̃β + γ̃γ − γ̃α
Ω̃|δω̃|
, (B-5)
tanφ2e
γ̃2β + γ̃
γ + γ̃
(Ω̃|δω̃|)2
+ tanφ2e − 1,
γ̃αγ̃β γ̃γ
(Ω̃|δω̃|)3
tanφ2e
4γ̃αγ̃β γ̃γ
(Ω̃|δω̃|)3
tanφe
+ tanφ
2MR2Ω2c
κ̃2ω̃αω̃βω̃γC(Te)
(γ̃αγ̃β γ̃γ)
(Ω̃|δω̃|)4
tanφe
tanφ2e
4γ̃αγ̃β γ̃γ
(Ω̃|δω̃|)3
tanφe
Ω̃|δω̃|C(Te)
The eigenvalues can be approximated as
λ1,2 ≈ −
− ǫ± i
ǫ2 + w2
, (B-6)
λ3,4 ≈ ǫ± iw,
λ5 ≈ −
where ǫ = (a2 − a3a4)/a4 and w =
a1/a3. The system
is unstable when a2 − a3a4 > 0 or a0 < 0. The first
two eigenvalues will have a negative real part as long as
γ̃β + γ̃γ > γ̃α. If the heating compensates the cooling of
the star a0 ≈ 0 and becomes negative if the star can not
reach thermal equilibrium. The other critical stability
condition a2 − a3a4 = 0 can be written as
Ω̃|δω̃|
[1+Γβ+Γγ−(Γ2β+Γ2γ)−(Γβ−Γγ)2(Γβ+Γγ)] = 0,
(B-7)
where Γβ = γβ/γα and Γγ = γγ/γα. Note that we have
ignored the smaller terms of orderO([γ̃α/(Ω̃|δω̃|)]5). This
condition can be rewritten by defining variables D1 =
Γβ + Γγ and D2 = Γβ − Γγ
2 + 2D1 −D21 −D22 − 2D22D1 = 0. (B-8)
If D2 = 0 then the equation has one solutionD1 = 1+
for D1 > 2, which corresponds to Γ = Γβ = Γγ = 1.37
and matches the result of Wersinger et al. [35]. For the
viscosity we consider (see Sec. II D) a2 − a3a4 < 0.
[1] R. Wagoner, Astrophys. J 578, L63 (2002).
[2] S. Chandrasekhar, Phys. Rev. Lett. 24, 611 (1970).
[3] J. L. Friedman and B. F. Schutz, Astrophys. J. 222, 281
(1978). J. L. Friedman and B. F. Schutz, Astrophys. J.
221, 937 (1978).
[4] N. Andersson, Astrophys. J. 502, 708 (1998).
[5] J. Friedman and S. Morsink, Astrophys. J. 502,
714(1998).
[6] L. Lindblom, B. J. Owen, and S. M. Morsink, Phys. Rev.
Lett 80, 4843 (1998).
[7] L. Bildsten, Astrophys. J. 501, L89 (1998).
[8] N. Andersson, K. D. Kokkotas, N. Stergioulas, Astro-
phys. J. 516, 307 (1999).
[9] N. Andersson, K. Kokkotas, and B. F. Schutz, Astrophys.
J. 510, 846 (1999).
[10] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro-
phys. J 423 L117 (1994).
[11] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro-
phys. J 424 823 (1994).
[12] J. W. T. Hessels et al. Science 311 1901 (2006).
[13] J. E. Grindlay, Science 311, 1876 (2006).
[14] D. C. Backer et al., Nature 300, 615 (1982).
[15] D. Chakrabarty et al., Nature 424, 42 (2003).
[16] D. Chakrabarty, Astron. Soc. Pac. Conf. Series 328, 279
(2005).
[17] B. J. Owen et al., Phys. Rev. D 58, 084020 (1998).
[18] Y. Levin, Astrophys. J 517, 328 (1999).
[19] P. Arras et al., Astrophys. J 591, 1129 (2003).
[20] J. Heyl, Astrophys. J 574, L57 (2002).
[21] P.B. Jones, Astrophys. Lett. 5, 33 (1970). P.B. Jones,
Proc, Roy. Soc. (London) A323, 111 (1971). P.B. Jones,
Phys. Rev. Lett. 86, 1384 (2001). P.B. Jones, Phys. Rev.
D64, 084003 (2001).
[22] L. Lindblom and B. J. Owen, Phys. Rev. D65, 063006
(2002), astro-ph/0110558.
[23] P. Haensel, K. P. Levenfish, and D. G. Yakovlev, Astron.
and Astrophys. 381, 1080 (2002), astro-ph/0110575.
[24] M. Nayyar and B. J. Owen, Phys. Rev. D 73 (2006)
084001, astro-ph/0512041.
[25] N. Andersson, D. I. Jones, and K. D. Kokkotas, MNRAS
337, 1224 (2002).
[26] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys. Rev.
D70 (2004) 121501, gr-qc/0406085.
[27] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys.Rev.
D70 (2004) 124017, gr-qc/0409048.
[28] J. Brink, S. A. Teukolsky, and I. Wasserman, Phys.Rev.
D71 (2005) 064029, gr-qc/0410072.
[29] A. K. Schenk, P. Arras, E. E. Flanagan, S. A. Teukol-
sky, I. Wasserman, Phys.Rev. D65 (2001) 024001,
gr-qc/0101092.
[30] L. Bildsten and G. Ushomirsky, Astrophys. J 529, L33
(2000).
[31] Y. Levin and G. Ushomirsky, MNRAS 322, 515 (2001).
[32] S. Yohida and U. Lee, Astrophys. J 546, 1121 (2001).
[33] K. Glampedakis and N. Andersson, astro-ph/0607105,
astro-ph/0411750.
[34] Y. S. Dimant, Phys. Rev. Lett. 84, 622 (2000).
[35] J. Wersinger, J. Finn, and E. Ott, Phys. Fluids 23, 1142
(1980).
[36] E. F. Brown, Ap. J 531, 988 (2000).
[37] H. Schatz, Phys. Rep. 294, 167 (1998).
[38] D. G. Yakovlev and K. P. Levenfish, Astron. Astrophys.
297, 717 (1995).
[39] D. G. Yakovlev and K. P. Levenfish, and Yu. A. Shibanov,
Soviet Phys.-Uspekhi, 42, 737 (1999).
[40] D. G. Yakovlev, A. D. Kaminker, and O. Y. Gnedin,
A&A, 379, L5 (2001).
[41] D. G. Yakovlev, and C. J. Pethick, Ann. Rev. Astron.
Astrophysics, 42, 169 (2004).
[42] G. Bryan, Philos. Trans. R. Soc. London A180, 187
(1889).
[43] C. Cutler and L. Lindblom, Astrophys. J 314, 234 (1987).
[44] K. D. Kokkotas and N. Stergioulas, Astron. and Astro-
phys. 341, 110 (1999).
[45] J. B. Kinney and G. Mendell, Phys.Rev. D67 024032
(2003).
[46] P. R. Brady, T. Creighton, C. Cutler, B. F. Schutz,
Phys. Rev. D 57, 2101 (1998), gr-qc/9702050. P. R.
Brady, T. Creighton, Phys. Rev. D 61, 082001 (2000),
gr-qc/9812014.
[47] B. J. Owen and L. Lindblom, Class.Quant.Grav. 19,
1247-1254 (2002), gr-qc/0111024.
[48] R. Prix, G. L. Comer, and N. Andersson, MNRAS 348,
625 (2004). N. Andersson and G. L. Comer, MNRAS
328,1129 (2001). N. Andersson, G. L. Comer and R. Prix,
MNRAS 354, 101 (2004).
[49] R. V. Wagoner, J. F. Hennawi, J. Liu, Proceedings of the
20th Texas Symposium on Relativistic Astrophysics, 781
(2001), astro-ph/0107229.
http://arxiv.org/abs/astro-ph/0110558
http://arxiv.org/abs/astro-ph/0110575
http://arxiv.org/abs/astro-ph/0512041
http://arxiv.org/abs/gr-qc/0406085
http://arxiv.org/abs/gr-qc/0409048
http://arxiv.org/abs/gr-qc/0410072
http://arxiv.org/abs/gr-qc/0101092
http://arxiv.org/abs/astro-ph/0607105
http://arxiv.org/abs/astro-ph/0411750
http://arxiv.org/abs/gr-qc/9702050
http://arxiv.org/abs/gr-qc/9812014
http://arxiv.org/abs/gr-qc/0111024
http://arxiv.org/abs/astro-ph/0107229
