AMR simulations of the low T/|W | bar-mode
instability of neutron stars
Pablo Cerdá-Durán, Vicent Quilis, and José A. Font
Departamento de Astronomı́a y Astrof́ısica, Universidad de Valencia, Dr. Moliner
50, 46100 Burjassot (Valencia), Spain
Abstract
It has been recently argued through numerical work that rotating stars with a
high degree of differential rotation are dynamically unstable against bar-mode de-
formation, even for values of the ratio of rotational kinetic energy to gravitational
potential energy as low as O(0.01). This may have implications for gravitational
wave astronomy in high-frequency sources such as core collapse supernovae. In this
paper we present high-resolution simulations, performed with an adaptive mesh re-
finement hydrodynamics code, of such low T/|W | bar-mode instability. The complex
morphological features involved in the nonlinear dynamics of the instability are re-
vealed in our simulations, which show that the excitation of Kelvin-Helmholtz-like
fluid modes outside the corotation radius of the star leads to the saturation of the
bar-mode deformation. While the overall trends reported in an earlier investigation
are confirmed by our work, we also find that numerical resolution plays an impor-
tant role during the long-term, nonlinear behaviour of the instability, which has
implications on the dynamics of rotating stars and on the attainable amplitudes of
the associated gravitational wave signals.
Key words: gravitational waves, hydrodynamics, instabilities, stars: neutron stars:
rotation
PACS: 97.60.Jd, 04.30.-w, 95.30.Lz
1 Introduction
Neutron stars following a core collapse supernova are rotating at birth and
can be subject to various nonaxisymmetric instabilities (see e.g. [1] for a re-
view). Among those, if the rotation rate is high enough so that the ratio of
rotational kinetic energy T to gravitational potential energy W , β ≡ T/|W |,
exceeds the critical value βd ∼ 0.27, inferred from studies with incompressible
Maclaurin spheroids, the star is subject to a dynamical bar-mode (l = m = 2
Preprint submitted to Elsevier 30 July 2021
http://arxiv.org/abs/0704.0356v1
f -mode) instability driven by hydrodynamics and gravity. Its study is highly
motivated nowadays as such an instability bears important implications in
the prospects of detection of gravitational radiation from newly-born rapidly
rotating neutron stars.
Simulations of the dynamical bar-mode instability are available in the litera-
ture, both using simplified models based on equilibrium stellar configurations
perturbed with suitable eigenfunctions [2,3,4,5], and more involved models for
the core collapse scenario [6,7,8,9], and in either case both in Newtonian grav-
ity and general relativity. Due to its superior simplicity the former approach
has received much more attention, notwithstanding that the conclusions drawn
from perturbed stellar models may not be straightforwardly extended to the
collapse scenario.
Newtonian simulations of triaxial instabilities following core collapse were first
performed by [6]. These showed that the bar-mode instability sets in when
β ≫ 0.27 and when the progenitor rotates rapidly and highly differentially.
Such conditions are met when the (artificial) depletion of internal energy to
trigger the collapse is large enough to produce a very compact core for which a
significant spun-up can be achieved. More recently, three-dimensional simula-
tions of the core collapse of rotating polytropes in general relativity have been
performed by [7]. These authors studied the evolution of the bar-mode insta-
bility starting with axisymmetric core collapse initial models which reached
values of β ∼ 0.27 during the infall phase. These simulations showed that the
maximum value of β achieved during collapse and bounce depends strongly
on the velocity profile, the total mass of the initial core, and on the equa-
tion of state. In agreement with the findings from the Newtonian simula-
tions of [6], the bar-mode instability sets in if the progenitor rotates rapidly
(0.01 ≤ β ≤ 0.02) and has a high degree of differential rotation. In addition,
the artificial depletion of pressure and internal energy to trigger the collapse,
leading to a compact core which subsequently spins up, also plays a key role
in general relativity for a noticeable growth of the bar-mode instability.
Whether the requirements inferred from numerical simulations are at all met
by the collapse progenitors remains unclear. As shown by [10] magnetic torques
can spin down the core of the progenitor, which leads to slowly rotating neu-
tron stars at birth (∼ 10 − 15ms). The most recent, state-of-the-art compu-
tations of the evolution of massive stars, which include angular momentum
redistribution by magnetic torques and spin estimates of neutron stars at
birth [11,12], lead to core collapse progenitors which do not seem to rotate
fast enough to guarantee the unambiguous growth of the canonical bar-mode
instability. Rapidly-rotating cores might be produced by an appropriate mix-
ture of high progenitor mass (M > 25M⊙) and low metallicity (N. Stergioulas,
private communication). In such case the progenitor could by-pass the Red
Supergiant phase in which the differential rotation of the core produces a
magnetic field by dynamo action which couples the core to the outer layers
of the star, transporting angular momentum outwards and spinning down the
core. According to [13] about 1% of all stars with M > 10M⊙ will produce
rapidly-rotating cores.
On the other hand, Newtonian simulations of the bar-mode instability from
perturbed equilibrium models of rotating stars have shown that βd ∼ 0.27
independent of the stiffness of the equation of state provided the star is not
strongly differentially rotating. The relativistic simulations of [5] yielded a
value of β ∼ 0.24 − 0.25 for the onset of the instability, while the dynamics
of the process closely resembles that found in Newtonian theory, i.e. unstable
models with large enough β develop spiral arms following the formation of
bars, ejecting mass and redistributing the angular momentum. As the degree
of differential rotation becomes higher Newtonian simulations have also shown
that βd can be as low as 0.14 [14]. More recently [15,16] have reported that
rotating stars with an extreme degree of differential rotation are dynamically
unstable against bar-mode deformation even for values of β of O(0.01).
Given its recent discovery and its potential astrophysical implications for post-
bounce core collapse dynamics and gravitational wave astronomy, we present
in this paper high resolution simulations of such low T/|W | bar-mode in-
stabilities. This work is further motivated in the light of the few numerical
simulations available in the literature. Our main goal is to revisit the simula-
tions by [15] on the low T/|W | bar-mode instability, and particularly to check
how sensitive the onset and development of the instability is to numerical
issues such as grid resolution. To this aim we perform Newtonian hydrody-
namical simulations of a subset of models analyzed by [15] using an adaptive
mesh refinement (AMR) code [17] which allows us to perform such three di-
mensional simulations with the highest resolution ever used. Our simulations
reveal the complex morphological features involved in the nonlinear dynamics
of the instability, where the excitation of Kelvin-Helmholtz-like fluid modes
influences the saturation of the bar-mode deformation. We advance that while
the overall trends found by [15] are confirmed by our work, the resolution
employed in the simulations does play a key role for the long-term behaviour
of the instability and for the nonlinear dynamics of rotating stars, which has
implications on the attainable amplitudes of the associated gravitational wave
signals. We note that we plan to upgrade the existing AMR code to account
for the effects of magnetic fields in order to attempt the current study in a
more realistic setup. The present work is a step towards that goal.
The paper is organized as follows: Section 2 gives a brief overview of the
equations to solve. Their solution is outlined in Section 3 which also contains
the bare details of the AMR code. The results of the simulations are discussed
in Section 4. Finally Section 5 presents our conclusions.
2 Mathematical framework
The evolution of a self-gravitating ideal fluid in the Newtonian limit is de-
scribed by the hydrodynamics equations and Poisson’s equation:
+∇ · (ρv) = 0 (1)
+ (v · ∇)v = −1
∇p−∇φ (2)
+∇ · [(E + p)v] = −ρv∇φ (3)
∇2φ = 4πGρ (4)
where x, v = dx
= (vx, vy, vz), and φ(t,x) are, respectively, the Eulerian
coordinates, the velocity, and the Newtonian gravitational potential. The total
energy density, E = ρǫ + 1
ρv2 , is defined as the sum of the thermal energy,
ρǫ, where ρ is the mass density and ǫ is the specific internal energy, and the
kinetic energy (where v2 = v2x + v
y + v
z). Pressure gradients and gravitational
forces are the responsible for the evolution. An equation of state p = p(ρ, ǫ)
closes the system. We use an ideal gas equation of state p = (Γ − 1)ρǫ with
Γ = 2.
The hydrodynamics equations, Eqs. (1–3), can be rewritten in flux-conservative
form:
∂f(u)
∂g(u)
∂h(u)
= s(u) (5)
where u is the vector of unknowns (conserved variables):
u = [ρ, ρvx, ρvy, ρvz, E] . (6)
The three flux functions Fα ≡ {f , g,h} in the spatial directions x, y, z, respec-
tively, are defined by
f(u) =
ρvx, ρv
x + p, ρvxvy, ρvxvz, (E + p)vx
g(u)=
ρvy, ρvxvy, ρv
y + p, ρvyvz, (E + p)vy
h(u) =
ρvz, ρvxvz, ρvyvz, ρv
z + p, (E + p)vz
and the source terms s are given by
Table 1
Overview of the initial models and results of the simulations. The rows report
the name of the model, the ratio of equatorial-to-polar radii (re/rp), the degree of
differential rotation (Â), the ratio of kinetic to potential energy (T/|W |), the size of
the computational grid (L) and the location of the corotation radius (rc) for the two
resolutions used: high (AMR H) and low (AMR L). In models R1H and R2H the
corotation radius lies outside the star. The real (frequency) and imaginary (growth
rate) parts of the bar-mode σ2 are shown, for the low and high resolution simulation
in comparison with the numerical results and linear analysis by [15]. Note that for
model D3 no linear analysis results are available.
Model D1 D2 D3 R1 R2
re/rp 0.805 0.605 0.305 0.305 0.255
Â 0.3 0.3 0.3 1.0 1.0
T/|W | 0.039 0.085 0.149 0.253 0.275
L/re 4.06 3.73 3.21 4.25 4.03
rc/re AMR L 0.38 0.47 0.58
AMR H 0.36 0.48 0.56 - -
Re(σ2)/Ω0 AMR L 0.76 0.58 0.41
AMR H 0.81 0.55 0.43 - 0.82
Shibata 0.80 0.60 0.45 0.92 0.75
linear 0.80 0.58 - 0.92 0.75
Im(σ2)/Ω0 AMR L 0.0042 0.0154 0.0200
AMR H 0.0089 0.0190 0.0240 0.0005 0.1960
Shibata 0.009-0.013 0.019-0.021 0.013 <0.002 0.23
linear 0.015 0.021 - <0.002 0.20
s(u) =
, ρvx
− ρvy
− ρvz
. (10)
System (5) is a three-dimensional hyperbolic system of conservation laws with
sources s(u).
3 Numerical approach
For our study of the low T/|W | bar-mode instability we perform high-resolution
simulations of rotating neutron stars using a Newtonian AMR hydrodynamics
code called MASCLET [17]. The implementation of the AMR technique in the
code follows the procedure developed by [18]. The hydrodynamics equations
are solved using a high-resolution shock-capturing scheme based upon Roe’s
Riemann solver and second-order cell reconstruction procedures, while Pois-
son’s equation for the gravitational field is solved using multigrid techniques.
The accuracy and performance of the MASCLET code has been assessed in a
number of tests [17]. We note that the code was originally designed for cos-
mological applications, and here it is applied to simulations of self-gravitating
stellar objects for the first time.
The simulations are performed with two different grid resolutions. The low
resolution grid consists of a box of size L with 1283 zones, yielding a fixed
resolution of L/128. We note that the effective resolution of our coarse grid
is comparable to that used by [15]. Correspondingly, the high resolution grid
consists of a base coarse grid of 1283 cells, and one level of refinement composed
of patches with maximum size of 643 cells (323 coarse cells). This yields a grid
resolution on the finest grid of L/256. This resolution is enough to resolve the
structures simulated, and hence no deeper refinement levels are needed. The
patches are dynamically allocated covering those regions of the star where the
highest resolution is required (highest densities). Typically only one patch is
needed for spheroidal models, and 4-8 in models with toroidal topology. The
use of AMR techniques in our high resolution simulations, allows us to save
about a factor 4 in CPU time and memory with respect to a unigrid simulation
with 2563 cells. No symmetries are imposed in the simulations. To the best of
our knowledge, in the investigations of the bar-mode instability performed by
previous groups, grid resolutions as high as the ones we use here were never
employed.
As customary in grid-based codes [19,20] the vacuum surrounding the star
is filled with a tenuous numerical atmosphere with density ρ/ρmax ≈ 10−12
and zero velocities, ρmax being the maximun density. Every grid cell with
ρ/ρmax < 10
−6 is reset to the atmosphere values. A correct treatment of the
atmosphere is essential for an accurate description of the stellar dynamics and
correct computation of the growth rates of unstable modes. We have checked
that values for the atmosphere higher than those we chose or a free evolution of
the atmosphere altogether, lead to remarkable changes in the mode behaviour,
growth rates, and frequencies. We have also checked that lower values for the
atmosphere do not produce those changes, which ensures that our evolutions
are not affected by the atmosphere values used in the simulations.
4 Results
4.1 Initial data
Differentially rotating stellar models in equilibrium are built according to the
method of [21], and used as initial data for the AMR evolution code. The stars
obey a polytropic equation of state P = KρΓ with index Γ = 2. As [15] the
profile of the angular velocity Ω is given by
(̟/re)2 + Â2
, (11)
where re is the equatorial radius of the star, Ω0 is the central angular ve-
locity, ̟ is the distance to the rotation axis, and Â parametrizes the degree
of differential rotation, from Â ≪ 1 for highly differentially rotating stars to
Â → ∞ for rigidly rotating stars. For comparison purposes these parameters
are chosen as in some of the models of [15], and are summarized in Table 1.
Models labelled D rotate with a high degree of differential rotation, as Â = 0.3,
and may therefore be subject to the low T/|W | bar-mode instability. We also
consider models almost rigidly rotating, labelled R, prone to experience the
“classical” bar-mode instability. Labels L and H in the models refer to low
and high resolution respectively.
Following [15] we perturb the initial density profile ρ(0) according to
ρ = ρ(0)
1 + δ
x2 − y2
, (12)
the perturbation of the pressure given by the equation of state accordingly.
A perturbation amplitude δ = 0.1 is used in all our simulations. As we show
below this form of the perturbation excites the l = m = 2 bar-mode. In
addition, grid discretization can leak small amounts of energy to all other
possible modes, which could in principle grow provided they were unstable
and the simulations were carried on for sufficiently long times.
4.2 Stability analysis
To compare with [15] we calculate the distortion parameters η+ and η× (and
η = (η2+ + η
)1/2) defined as
0 0.5 1 1.5 2
Ω / Ω0
0 0.21 0.43 0.65 0.86 1.1 1.3 1.5 1.7 1.9
Fig. 1. Power spectra of Am from m = 1 to m = 8 for model D3H.
Ixx − Iyy
Ixx + Iyy
, η× ≡
Ixx + Iyy
, (13)
where Iij(i, j = x, y, z) is the mass-quadrupole moment
Iij =
dx3ρ xixj . (14)
For the study of the growth rate and interaction of the different angular modes
within the star is useful to calculate the global quantity
dx3ρ(x) e−imϕ, (15)
and Am ≡ Am/A0. We follow the time evolution of modes with m ranging
from 1 to 8. Since our initial equilibrium models are axisymmetric and have
equatorial plane symmetry, all Am are zero initially, but once perturbed all
0 20 40 60 80 100
Fig. 2. Evolution of η for models R1H (upper panel) and R2H (lower panel). Expo-
nential fits to the peaks in the growing phase are overplotted as solid lines.
initial models exhibit a dominant m = 2 component. Assuming that the modes
behave as e−i(σmt−mϕ), the real part of σm can be obtained by Fourier trans-
forming Am. In particular Re(σ2), the bar-mode frequency, can be extracted
from either A2 or η as both represent the same mode. This is the dominant
mode in all our simulations and its frequency and growth rate are given in
Table 1. The latter corresponds to the imaginary part of σ2, which is calcu-
lated fitting an exponential to the peak values of η in the growing phase of
the evolution until the modes saturate. Other modes are also identified in the
simulations for values of Am with lower amplitudes. We have checked that
these modes are harmonics of the l = m = 2 mode so that they follow to good
accuracy the relation σm = mσp, σp being the pattern frequency, calculated
as σp = σ2/2. This is shown for model D3H in Fig. 1 which displays the spec-
trum of Am from m = 1 to m = 8 (in arbitrary units). The vertical dashed
lines in this figure indicate the location of the integer multiples of the pattern
frequency σp, their values indicated on the axis at the top of the figure. Each
spectrum for each mode is normalized to its own maximum for plotting pur-
poses. Note that the lower the mode amplitude the noisier the spectrum and
the less accurate the relation σm = mσp.
For the models of our sample subject to the “clasical” bar-mode deformation
(R1H and R2H), our simulations yield a value of β between 0.253 and 0.275,
in good agreement with the critical value for the onset of the dynamical bar-
mode instability. Model R1H is stable and model R2H is unstable. The growth
rates and frequencies reported in Table 1 agree with those of [15]. Note that
0 100 200 300
Fig. 3. Evolution of η for models D1 (upper panel), D2 (central panel) and D3 (lower
panel). Dashed lines correspond to low resolution and solid lines to high resolution.
Exponential fits to the peaks in the growing phase are overplotted as solid lines.
for model R1H, which is stable, the frequency for the m = 2 mode cannot be
computed. The time evolution of η for these two models is displayed in Fig. 2.
For the unstable model R2H, our simulations show the formation of a bar
which saturates for values of η+ and η× close to 1, i. e. in the full nonlinear
regime.
Fig. 3 shows the time evolution of η for models D in our sample, prone to suffer
the low T/|W | bar-mode instability. Solid lines correspond to high resolution
simulations and dashed lines to low resolution. For all three models the pattern
frequencies σp are such that there exists a corotation radius inside the star,
i.e. a radius at which the bar-mode rotates with the same angular velocity as
the fluid. The location of the corotation radius for all models of our sample
is reported in Table 1. As recently discussed by [22] the existence of such
corotation radius is a potential requirement for the ocurrence of the instability.
As becomes clear from Fig. 3, all models are unstable but grid resolution has
an important effect on the saturation of the instability once the nonlinear
phase has been reached, as well as in the long-term dynamics of the stars.
In the linear phase of models D1H and D2H, the growth rates and frequencies
agree with the results of [15] in both, the numerical simulations and the linear
analysis (see Table 1). In the linear phase of model D3H, our frequencies are
similar to the numerical results of [15], although our growth rates are about
a factor two larger. We emphasize that no results are reported in the linear
analysis for this model in the work of [15], and therefore this discrepancy can
be an effect of the resolution used or of the characteristics of each numerical
code. Increasing resolution leads to similar results in the frequencies but to
higher growth rates.
In the nonlinear phase, models D1 and D3 behave similarly for the two resolu-
tions used (see Fig. 3), and also similarly to the results by [15] (compare with
Fig. 3 of that paper). For model D2 we observe a radical change of behavior in
the nonlinear phase of the mode evolution depending on the grid resolution.
This has implications on the long-term dynamics of the star and, in particu-
lar, on the attainable amplitudes of the gravitational radiation emitted, as we
discuss below.
It is worth mentioning the possibility that the unstable mode at the start of
model D2H might excite some other mode in the corotation band, which could
not otherwise be excited for lower grid resolution. As discussed by [23,24] in
their study of differentially rotating shells, there are many zero-step modes in
the band, so that the whole continuous spectrum could potentially be excited.
In such case these modes would have very slow power-law growth.
For all our models we have checked mass conservation along the evolution. The
worst results are obtained for model D3H, for which mass is conserved within
2.5% error when the instability saturates. At the end of the simulation (after
48 orbital periods and 25000 iterations in the coarsest grid) the error has grown
to only 6%. For all other models mass conservation is even more accurate. Note
that these errors are within the round-off error of the code, and it is not related
to the conservation properties of the numerical scheme itself. For a regular grid
with 1283 cells and a simulation employing 25000 iterations, the accumulated
round-off error (binomial distribution) using single-precision arithmetics, is
about
1283 × 25000 × 10−8 = 0.0023 = 0.23%. Correspondingly, for a 2563
grid (with twice the number of iterations for the simulation) the error is about
0.9%. Taking into account that this error affects the nonlinear evolution of the
system, it is not surprising to have an error at the level of a few percent by
the end of our high resolution simulations, for all conserved quantities.
Figure 4 shows the evolution of Am for model D3 and for m ranging from
1 to 8 for our two resolutions. According to this figure, the only two modes
1e-06
0.0001
0 100 200
1e-06
0.0001
m=1, 3, 5, 7
m=1, 3, 5, 6, 7, 8
Fig. 4. Evolution of |Am| for model D3 with low resolution (top) and high resolution
(bottom). The m = 2 mode is represented with thick solid line, m = 4 with thin
solid line, m = 6 with dashed line, m = 8 with dot-dashed line, and all other odd
m with dotted lines.
relevant for the dynamics of the star are m = 2 and m = 4. All other modes
have smaller amplitudes and play no role in the dynamics. Note that for odd
m modes, the value of the integrated quantity Am, if close to zero, is extremely
sensitive to very small numerical asymmetries, which are induced by the patch
creation scheme of our AMR code. This explains the resolution differences in
the initial values for odd m modes in Fig. 4 (at t = 0 they start off at 10−8
level for the low resolution simulation), although they saturate at the same
value irrespective of the resolution.
An important diagnosis for the accuracy of the results is the location of the
center of mass during an evolution. The round-off error of the numerical code
imposes controlled errors in mass and linear momentum, which results in tiny
displacements of the center of mass. However small (one numerical cell in
our runs) this unphysical displacement may hinder the correct analysis of the
0 100 200 300
1e-06
0.0001
Fig. 5. Effects of the artificial displacement of the center of mass (of only one
numerical cell) on the time evolution of |A1| for model D2H. The thin solid line
shows a fictitiuos evolution resulting from the numerical artifact originated by the
center of mass displacement.
mode growth rates. For this reason all integrated quantities shown in Fig. 4
are computed after correcting for the displacement of the center of mass,
xnew = xold − xCM, in a post-processing stage of the data analysis. Were this
not done, a one-armed m = 1 mode would grow much faster than it should to
bring up fictitious features in the plots. This is shown for model D2H in Fig. 5.
The thick solid line in this figure corresponds to the evolution of the m = 1
mode taking into account the correction for the center of mass displacement,
while the thin solid line is the corresponding evolution of this mode without
the correction.
4.3 Gravitational waves
The growth and saturation of the instability is also imprinted on the gravita-
tional waves emitted. The gravitational waveforms h+ and h× for models D1,
D2, and D3, computed using the standard quadrupole formula, are shown in
Fig. 6. For a source of mass M located at a distance R those waveforms can
be calculated from the dimensionless waveform amplitudes a+ and a× as
h+,× = a+,×
sin2 θ
, (16)
using G = c = 1 units. The resulting chirp-like signal in all the models, partic-
ularly apparent for model D2L, indicates the presence of a bipolar distribution
of mass within the star (see Sec. 4.4).
-0.05
0 100 200 300
Fig. 6. Gravitational waves for models D1 to D3 extracted using the standard
quadrupole formula. Thick (thin) solid lines correspond to low (high) resolution.
Only the dimensionless waveform amplitude a+ is plotted.
As mentioned before, the effects of grid resolution on the evolution of the
nonlinear phase of the bar-mode are imprinted on the gravitational waveforms.
Thick solid lines in Fig. 6 are the waveforms which correspond to the low-
resolution models, and thin solid lines to the high-resolution counterparts.
The evolution of η for model D3, displayed in Fig. 3, shows little deviations
with grid resolution, and this translates into very similar gravitational wave
patterns (bottom panel of Fig. 6), the differences becoming more noticeable in
the nonlinear phase following saturation (Ω0t ≥ 75). For model D1 (top panel),
the differences also become more apparent at later times during the evolution,
in good agreement with the dissimilar behaviour of the matter dynamics in
this model, as encoded in the evolution of η in Fig. 3. As happens for model
D3 the first few cycles of the gravitational waveform, when the mode is still
in the linear phase, are accurately captured for both resolutions.
The major dependence of the waveform on the grid resolution is found for
model D2. Again, the linear phase for the growth of the bar deformation is
Ω 0 t=
70.40.0
Fig. 7. Snapshots of the density, vorticity, and specific angular momentum, for
model D3H, at three representative instants of the evolution. All snapshots show
slices of the stars in the equatorial plane. Quantities are normalized as follows:
max, rew
s , and l
ϕ/(rev
s ), where v
s is the initial velocity at the surface
of the star.
accurately captured irrespective of the resolution (and agrees with the per-
turbative results of [15]). This is signalled in the perfect overlapping of both
gravitational waveforms during the first three cycles (see the middle panel of
Fig. 6). However, the different nonlinear dynamics of the bar-mode deforma-
tion for this model, shown in the middle panel of Fig. 3, is severely imprinted
on the gravitational waveform. Model D2H emits gravitational waves which
have roughly one order of magnitude smaller amplitude than those computed
for the corresponding low resolution model.
4.4 Morphology
We next describe the morphological features encountered during the evolution
of some representative models. Fig. 7 shows three snaphsots of the evolution
of model D3H for the density (top), the azimuthal component of the vorticity,
~wϕ = (∇ × ~v)ϕ (middle), and the specific angular momentum, ~l = ~r × ~v
(bottom). From left to right the snapshots correspond to the initial time (Ω0t =
0), a time when the bar-mode instability is growing (Ω0t = 33.6), and the time
when the instability saturates (Ω0t = 70.4). Only the equatorial plane of the
stars is shown in all these plots. Animations of all simulations performed are
available at www.uv.es/∼cerdupa/bars/. We note that our AMR code is able
to dynamically place patches (e. g. between 4 and 8 in the D3H model) and
evolves the system with continuous matching between patches, as exemplified
in Fig. 7.
The evolution of model D3H shows that as the m = 2 mode grows the star
develops an ellipsoidal shape which remains spinning beyond saturation. Since
the low β m = 2 mode saturates at lower values (η ∼ 0.1) than the classical
bar-mode instability (η ∼ 1), no clear bars are visible in the density plot. At
late times (Ω0t > 100) a “boxy” structure becomes apparent as the m = 4
mode has grown to almost similar amplitude as the m = 2 mode (see anima-
tions and Fig. 4). No other global features can be seen, consistent with the
fact that |Am| ≪ 1 for all modes other than m = 2 and 4. The vorticity plot
shows that the m = 2 mode at Ω0t = 33.6 adopts the form of a two-armed
spiral winding up around the central parts of the star. As the mode begins to
saturate (Ω0t = 70.4) the spirals break apart into the outer layers in a turbu-
lent flow reminiscent of the (shear) Kelvin-Helmholtz instability, and shock as
they reach the atmosphere. These trends are also visible in the specific angular
momentum plot.
The presence of a corotation radius, at r/re = 0.56 for model D3H, seems
to play a role in the growth and saturation of the instability, in agreement
with the recent findings of [25]. As the bar-mode grows, pressure waves carry
angular momentum outside the corotation radius, which is deposited in the
outer layers of the star. This excites Kelvin-Helmholtz-like instabilities in the
fluid that break the mode outside the corotation radius. When this happens the
m = 2 instability stops growing and no more angular momentum is extracted.
Figure 8 shows late-time snapshots of the equatorial plane distribution of the
density perturbation, i.e. (ρ−ρ(0))/ρ(0)max, for models D2H and D3H. The times
are chosen well inside the nonlinear and saturation phase of the instability.
This figure helps to interpret the mode dynamics and its saturation along the
lines mentioned before: During the evolution the density perturbations are
shed in waves from the center towards the outer layers of the star. At late
times, when the instability saturates, such shedding stops, and the density
Fig. 8. Snapshots of the density perturbation at the equatorial plane for models D2H
and D3H. The white solid curves indicate the location of the corotation radius. The
white dashed boxes indicate the location of the patches for model D3H.
perturbation reaches the largest values outside the corotation radius (depicted
with white solid lines in Fig. 8), for either model.
We note in passing that the corotation radius in all our high resolution models
lies well inside the outer boundary of the finest box set up by the AMR
refinement pattern. (see, e.g. the white dashed boxes depicted in the right
panel of Fig. 8 indicating the location of the AMR patches for model D3H)
This rules out the possibility of a numerical artifact resulting from the patch
creation scheme of our AMR code being the cause for the different long-term
evolution between low and high resolution models, particularly noticeable for
model D2 in Fig. 3.
Finally, Fig. 9 shows a comparison between models D2L and D2H at Ω0t = 101
(i.e. well within the nonlinear phase), to highlight the effects of the numer-
ical resolution on the morphology. From top to bottom this panel shows a
schlieren plot (|∇ log ρ|) , ~wϕ, and ~l. The resolution differences in the evolu-
tion of model D2 become apparent from this figure. In particular, the “boxy”
structure becomes much more clearly visible in the low resolution simulation
(D2L), indicating an excessive growth rate of the m = 4 mode. The presence
of pressure waves is emphasized in the schlieren plot, very accurately captured
in model D2H. Those waves, once the flow is driven to turbulence past the
corotation radius, redistribute the angular momentum in the outer layers of
model D2L in a much more pronounced way than for model D2H.
Ω 0 t=
D2L D2H
Fig. 9. Resolution comparison between models D2L and D2H once the instability
has saturated. Only slices of the stars in the equatorial plane are shown.
5 Summary and outlook
We have presented AMR high-resolution simulations of the low T/|W | bar-
mode instability of extremely differentially rotating neutron stars. Our main
motivation has been to revisit the simulations by [15] on such instability,
assessing how sensitive the onset and development of the instability is to
numerical issues such as grid resolution. We have addressed the importance
of a correct treatment of delicate numerical aspects which may spoil three-
dimensional simulations in (Cartesian) grid-based codes, always hampered by
insufficient resolution, namely the handling of the low-density atmosphere sur-
rounding the star, the correction for the center of mass displacement, and the
mass and momentum conservation properties of the numerical scheme. Our
simulations have revealed the complex morphological features involved in the
nonlinear dynamics of the instability. We have found that in the nonlinear
phase of the evolution, the excitation of Kelvin-Helmholtz-like fluid modes
outside the corotation radii of the stellar models leads to the saturation of the
bar-mode deformation. While the overall trends reported in the investigation
of [15] are confirmed by our work, the resolution used to perform the simu-
lations may play a key role on the long-term behaviour of the instability and
on the nonlinear dynamics of rotating stars, which has only become apparent
for some specific models of our sample (namely model D2). This, in turn, has
implications on the attainable amplitudes of the associated gravitational wave
signals.
The work reported in this paper is a first step in our ongoing efforts of study-
ing the dynamical bar-mode instability within the magnetized core collapse
scenario.
Acknowledgements
The authors thank Harry Dimmelmeier, Nick Stergioulas, and Anna Wats
for useful comments. Research supported by the Spanish Ministerio de Edu-
cación y Ciencia (MEC; grants AYA2004-08067-C03-01, AYA2003-08739-C02-
02, AYA2006-02570). VQ is a Ramón y Cajal Fellow of the Spanish MEC.
Computations performed at the Servei d’Informática de la Universitat de
València (CERCA-CESAR).
References
[1] N. Stergioulas, Liv. Rev. Relativ. 6 (2003) 3
[2] J. E. Tohline, R. H. Durisen, & M. McCollough, ApJ 298 (1985) 220
[3] J. L. Houser, J. M. Centrella, & S. Smith, Phys. Rev. Lett. 72 (1994) 1314
[4] K. C. B. New, J. M. Centrella, & J. E. Tohline, Phys. Rev. D 62 (2000) 064019
[5] M. Shibata, T. W. Baumgarte, & S. L. Shapiro, ApJ 542 (2000) 453
[6] M. Rampp, E. Müller, & M. Ruffert, A&A 332 (1998) 969
[7] M. Shibata, & Y. Sekiguchi, Phys. Rev. D 71 (2005) 024014
[8] M. Saijo, Phys. Rev. D 71 (2005) 104038
[9] C. D. Ott, S. Ou, J. E. Tohline, & A. Burrows, ApJ 625 (2005) L119
[10] H. C. Spruit & E. S. Phinney, Nature 393 (1998) 139
[11] A. Heger, S. E. Woosley, & H. C. Spruit, ApJ 626 (2005) 350
[12] C. D. Ott, A. Burrows, T. A. Thompson, E. Livne, & R. Walder, ApJS 164
(2006) 130
[13] S. E. Woosley, & A. Heger, ApJ 637 (2006) 914
[14] J. M. Centrella, K. C. B. New, L. L. Lowe, & J. D. Brown, ApJ 550 (2001)
[15] M. Shibata, S. Karino, & Y. Eriguchi, MNRAS 334 (2002) L27
[16] M. Shibata, S. Karino, & Y. Eriguchi, MNRAS 343 (2003) 619
[17] V. Quilis, MNRAS 352 (2004) 1426
[18] M. J. Berger, P. Colella, J. Comp. Phys. 82 (1989) 64
[19] J. A. Font, M. Miller, W.-M. Suen, & M. Tobias, Phys. Rev. D 61 (2000) 0044011
[20] M. D. Duez, P. Marronetti, S. L. Shapiro, & T. W. Baumgarte, Phys. Rev. D
67 (2003) 024004
[21] Y. Eriguchi, & E. Müller, A&A, 147 (1984) 161
[22] A. L. Watts, N. Andersson, & D. I. Jones, ApJ 618 (2005) L37
[23] A. L. Watts, N. Andersson, H. Beyer, & B. F. Schutz, MNRAS 342 (2003) 1156
[24] A. L. Watts, N. Andersson, & R. L. Williams, MNRAS 350 (2004) 927
[25] M. Saijo & S. Yoshida, MNRAS 368 (2006) 1429
	Introduction
	Mathematical framework
	Numerical approach
	Results
	Initial data
	Stability analysis
	Gravitational waves
	Morphology
	Summary and outlook
	References

  It has been recently argued through numerical work that rotating stars with a
high degree of differential rotation are dynamically unstable against bar-mode
deformation, even for values of the ratio of rotational kinetic energy to
gravitational potential energy as low as O(0.01). This may have implications
for gravitational wave astronomy in high-frequency sources such as core
collapse supernovae. In this paper we present high-resolution simulations,
performed with an adaptive mesh refinement hydrodynamics code, of such low
T/|W| bar-mode instability. The complex morphological features involved in the
nonlinear dynamics of the instability are revealed in our simulations, which
show that the excitation of Kelvin-Helmholtz-like fluid modes outside the
corotation radius of the star leads to the saturation of the bar-mode
deformation. While the overall trends reported in an earlier investigation are
confirmed by our work, we also find that numerical resolution plays an
important role during the long-term, nonlinear behaviour of the instability,
which has implications on the dynamics of rotating stars and on the attainable
amplitudes of the associated gravitational wave signals.

Introduction
Neutron stars following a core collapse supernova are rotating at birth and
can be subject to various nonaxisymmetric instabilities (see e.g. [1] for a re-
view). Among those, if the rotation rate is high enough so that the ratio of
rotational kinetic energy T to gravitational potential energy W , β ≡ T/|W |,
exceeds the critical value βd ∼ 0.27, inferred from studies with incompressible
Maclaurin spheroids, the star is subject to a dynamical bar-mode (l = m = 2
Preprint submitted to Elsevier 30 July 2021
http://arxiv.org/abs/0704.0356v1
f -mode) instability driven by hydrodynamics and gravity. Its study is highly
motivated nowadays as such an instability bears important implications in
the prospects of detection of gravitational radiation from newly-born rapidly
rotating neutron stars.
Simulations of the dynamical bar-mode instability are available in the litera-
ture, both using simplified models based on equilibrium stellar configurations
perturbed with suitable eigenfunctions [2,3,4,5], and more involved models for
the core collapse scenario [6,7,8,9], and in either case both in Newtonian grav-
ity and general relativity. Due to its superior simplicity the former approach
has received much more attention, notwithstanding that the conclusions drawn
from perturbed stellar models may not be straightforwardly extended to the
collapse scenario.
Newtonian simulations of triaxial instabilities following core collapse were first
performed by [6]. These showed that the bar-mode instability sets in when
β ≫ 0.27 and when the progenitor rotates rapidly and highly differentially.
Such conditions are met when the (artificial) depletion of internal energy to
trigger the collapse is large enough to produce a very compact core for which a
significant spun-up can be achieved. More recently, three-dimensional simula-
tions of the core collapse of rotating polytropes in general relativity have been
performed by [7]. These authors studied the evolution of the bar-mode insta-
bility starting with axisymmetric core collapse initial models which reached
values of β ∼ 0.27 during the infall phase. These simulations showed that the
maximum value of β achieved during collapse and bounce depends strongly
on the velocity profile, the total mass of the initial core, and on the equa-
tion of state. In agreement with the findings from the Newtonian simula-
tions of [6], the bar-mode instability sets in if the progenitor rotates rapidly
(0.01 ≤ β ≤ 0.02) and has a high degree of differential rotation. In addition,
the artificial depletion of pressure and internal energy to trigger the collapse,
leading to a compact core which subsequently spins up, also plays a key role
in general relativity for a noticeable growth of the bar-mode instability.
Whether the requirements inferred from numerical simulations are at all met
by the collapse progenitors remains unclear. As shown by [10] magnetic torques
can spin down the core of the progenitor, which leads to slowly rotating neu-
tron stars at birth (∼ 10 − 15ms). The most recent, state-of-the-art compu-
tations of the evolution of massive stars, which include angular momentum
redistribution by magnetic torques and spin estimates of neutron stars at
birth [11,12], lead to core collapse progenitors which do not seem to rotate
fast enough to guarantee the unambiguous growth of the canonical bar-mode
instability. Rapidly-rotating cores might be produced by an appropriate mix-
ture of high progenitor mass (M > 25M⊙) and low metallicity (N. Stergioulas,
private communication). In such case the progenitor could by-pass the Red
Supergiant phase in which the differential rotation of the core produces a
magnetic field by dynamo action which couples the core to the outer layers
of the star, transporting angular momentum outwards and spinning down the
core. According to [13] about 1% of all stars with M > 10M⊙ will produce
rapidly-rotating cores.
On the other hand, Newtonian simulations of the bar-mode instability from
perturbed equilibrium models of rotating stars have shown that βd ∼ 0.27
independent of the stiffness of the equation of state provided the star is not
strongly differentially rotating. The relativistic simulations of [5] yielded a
value of β ∼ 0.24 − 0.25 for the onset of the instability, while the dynamics
of the process closely resembles that found in Newtonian theory, i.e. unstable
models with large enough β develop spiral arms following the formation of
bars, ejecting mass and redistributing the angular momentum. As the degree
of differential rotation becomes higher Newtonian simulations have also shown
that βd can be as low as 0.14 [14]. More recently [15,16] have reported that
rotating stars with an extreme degree of differential rotation are dynamically
unstable against bar-mode deformation even for values of β of O(0.01).
Given its recent discovery and its potential astrophysical implications for post-
bounce core collapse dynamics and gravitational wave astronomy, we present
in this paper high resolution simulations of such low T/|W | bar-mode in-
stabilities. This work is further motivated in the light of the few numerical
simulations available in the literature. Our main goal is to revisit the simula-
tions by [15] on the low T/|W | bar-mode instability, and particularly to check
how sensitive the onset and development of the instability is to numerical
issues such as grid resolution. To this aim we perform Newtonian hydrody-
namical simulations of a subset of models analyzed by [15] using an adaptive
mesh refinement (AMR) code [17] which allows us to perform such three di-
mensional simulations with the highest resolution ever used. Our simulations
reveal the complex morphological features involved in the nonlinear dynamics
of the instability, where the excitation of Kelvin-Helmholtz-like fluid modes
influences the saturation of the bar-mode deformation. We advance that while
the overall trends found by [15] are confirmed by our work, the resolution
employed in the simulations does play a key role for the long-term behaviour
of the instability and for the nonlinear dynamics of rotating stars, which has
implications on the attainable amplitudes of the associated gravitational wave
signals. We note that we plan to upgrade the existing AMR code to account
for the effects of magnetic fields in order to attempt the current study in a
more realistic setup. The present work is a step towards that goal.
The paper is organized as follows: Section 2 gives a brief overview of the
equations to solve. Their solution is outlined in Section 3 which also contains
the bare details of the AMR code. The results of the simulations are discussed
in Section 4. Finally Section 5 presents our conclusions.
2 Mathematical framework
The evolution of a self-gravitating ideal fluid in the Newtonian limit is de-
scribed by the hydrodynamics equations and Poisson’s equation:
+∇ · (ρv) = 0 (1)
+ (v · ∇)v = −1
∇p−∇φ (2)
+∇ · [(E + p)v] = −ρv∇φ (3)
∇2φ = 4πGρ (4)
where x, v = dx
= (vx, vy, vz), and φ(t,x) are, respectively, the Eulerian
coordinates, the velocity, and the Newtonian gravitational potential. The total
energy density, E = ρǫ + 1
ρv2 , is defined as the sum of the thermal energy,
ρǫ, where ρ is the mass density and ǫ is the specific internal energy, and the
kinetic energy (where v2 = v2x + v
y + v
z). Pressure gradients and gravitational
forces are the responsible for the evolution. An equation of state p = p(ρ, ǫ)
closes the system. We use an ideal gas equation of state p = (Γ − 1)ρǫ with
Γ = 2.
The hydrodynamics equations, Eqs. (1–3), can be rewritten in flux-conservative
form:
∂f(u)
∂g(u)
∂h(u)
= s(u) (5)
where u is the vector of unknowns (conserved variables):
u = [ρ, ρvx, ρvy, ρvz, E] . (6)
The three flux functions Fα ≡ {f , g,h} in the spatial directions x, y, z, respec-
tively, are defined by
f(u) =
ρvx, ρv
x + p, ρvxvy, ρvxvz, (E + p)vx
g(u)=
ρvy, ρvxvy, ρv
y + p, ρvyvz, (E + p)vy
h(u) =
ρvz, ρvxvz, ρvyvz, ρv
z + p, (E + p)vz
and the source terms s are given by
Table 1
Overview of the initial models and results of the simulations. The rows report
the name of the model, the ratio of equatorial-to-polar radii (re/rp), the degree of
differential rotation (Â), the ratio of kinetic to potential energy (T/|W |), the size of
the computational grid (L) and the location of the corotation radius (rc) for the two
resolutions used: high (AMR H) and low (AMR L). In models R1H and R2H the
corotation radius lies outside the star. The real (frequency) and imaginary (growth
rate) parts of the bar-mode σ2 are shown, for the low and high resolution simulation
in comparison with the numerical results and linear analysis by [15]. Note that for
model D3 no linear analysis results are available.
Model D1 D2 D3 R1 R2
re/rp 0.805 0.605 0.305 0.305 0.255
Â 0.3 0.3 0.3 1.0 1.0
T/|W | 0.039 0.085 0.149 0.253 0.275
L/re 4.06 3.73 3.21 4.25 4.03
rc/re AMR L 0.38 0.47 0.58
AMR H 0.36 0.48 0.56 - -
Re(σ2)/Ω0 AMR L 0.76 0.58 0.41
AMR H 0.81 0.55 0.43 - 0.82
Shibata 0.80 0.60 0.45 0.92 0.75
linear 0.80 0.58 - 0.92 0.75
Im(σ2)/Ω0 AMR L 0.0042 0.0154 0.0200
AMR H 0.0089 0.0190 0.0240 0.0005 0.1960
Shibata 0.009-0.013 0.019-0.021 0.013 <0.002 0.23
linear 0.015 0.021 - <0.002 0.20
s(u) =
, ρvx
− ρvy
− ρvz
. (10)
System (5) is a three-dimensional hyperbolic system of conservation laws with
sources s(u).
3 Numerical approach
For our study of the low T/|W | bar-mode instability we perform high-resolution
simulations of rotating neutron stars using a Newtonian AMR hydrodynamics
code called MASCLET [17]. The implementation of the AMR technique in the
code follows the procedure developed by [18]. The hydrodynamics equations
are solved using a high-resolution shock-capturing scheme based upon Roe’s
Riemann solver and second-order cell reconstruction procedures, while Pois-
son’s equation for the gravitational field is solved using multigrid techniques.
The accuracy and performance of the MASCLET code has been assessed in a
number of tests [17]. We note that the code was originally designed for cos-
mological applications, and here it is applied to simulations of self-gravitating
stellar objects for the first time.
The simulations are performed with two different grid resolutions. The low
resolution grid consists of a box of size L with 1283 zones, yielding a fixed
resolution of L/128. We note that the effective resolution of our coarse grid
is comparable to that used by [15]. Correspondingly, the high resolution grid
consists of a base coarse grid of 1283 cells, and one level of refinement composed
of patches with maximum size of 643 cells (323 coarse cells). This yields a grid
resolution on the finest grid of L/256. This resolution is enough to resolve the
structures simulated, and hence no deeper refinement levels are needed. The
patches are dynamically allocated covering those regions of the star where the
highest resolution is required (highest densities). Typically only one patch is
needed for spheroidal models, and 4-8 in models with toroidal topology. The
use of AMR techniques in our high resolution simulations, allows us to save
about a factor 4 in CPU time and memory with respect to a unigrid simulation
with 2563 cells. No symmetries are imposed in the simulations. To the best of
our knowledge, in the investigations of the bar-mode instability performed by
previous groups, grid resolutions as high as the ones we use here were never
employed.
As customary in grid-based codes [19,20] the vacuum surrounding the star
is filled with a tenuous numerical atmosphere with density ρ/ρmax ≈ 10−12
and zero velocities, ρmax being the maximun density. Every grid cell with
ρ/ρmax < 10
−6 is reset to the atmosphere values. A correct treatment of the
atmosphere is essential for an accurate description of the stellar dynamics and
correct computation of the growth rates of unstable modes. We have checked
that values for the atmosphere higher than those we chose or a free evolution of
the atmosphere altogether, lead to remarkable changes in the mode behaviour,
growth rates, and frequencies. We have also checked that lower values for the
atmosphere do not produce those changes, which ensures that our evolutions
are not affected by the atmosphere values used in the simulations.
4 Results
4.1 Initial data
Differentially rotating stellar models in equilibrium are built according to the
method of [21], and used as initial data for the AMR evolution code. The stars
obey a polytropic equation of state P = KρΓ with index Γ = 2. As [15] the
profile of the angular velocity Ω is given by
(̟/re)2 + Â2
, (11)
where re is the equatorial radius of the star, Ω0 is the central angular ve-
locity, ̟ is the distance to the rotation axis, and Â parametrizes the degree
of differential rotation, from Â ≪ 1 for highly differentially rotating stars to
Â → ∞ for rigidly rotating stars. For comparison purposes these parameters
are chosen as in some of the models of [15], and are summarized in Table 1.
Models labelled D rotate with a high degree of differential rotation, as Â = 0.3,
and may therefore be subject to the low T/|W | bar-mode instability. We also
consider models almost rigidly rotating, labelled R, prone to experience the
“classical” bar-mode instability. Labels L and H in the models refer to low
and high resolution respectively.
Following [15] we perturb the initial density profile ρ(0) according to
ρ = ρ(0)
1 + δ
x2 − y2
, (12)
the perturbation of the pressure given by the equation of state accordingly.
A perturbation amplitude δ = 0.1 is used in all our simulations. As we show
below this form of the perturbation excites the l = m = 2 bar-mode. In
addition, grid discretization can leak small amounts of energy to all other
possible modes, which could in principle grow provided they were unstable
and the simulations were carried on for sufficiently long times.
4.2 Stability analysis
To compare with [15] we calculate the distortion parameters η+ and η× (and
η = (η2+ + η
)1/2) defined as
0 0.5 1 1.5 2
Ω / Ω0
0 0.21 0.43 0.65 0.86 1.1 1.3 1.5 1.7 1.9
Fig. 1. Power spectra of Am from m = 1 to m = 8 for model D3H.
Ixx − Iyy
Ixx + Iyy
, η× ≡
Ixx + Iyy
, (13)
where Iij(i, j = x, y, z) is the mass-quadrupole moment
Iij =
dx3ρ xixj . (14)
For the study of the growth rate and interaction of the different angular modes
within the star is useful to calculate the global quantity
dx3ρ(x) e−imϕ, (15)
and Am ≡ Am/A0. We follow the time evolution of modes with m ranging
from 1 to 8. Since our initial equilibrium models are axisymmetric and have
equatorial plane symmetry, all Am are zero initially, but once perturbed all
0 20 40 60 80 100
Fig. 2. Evolution of η for models R1H (upper panel) and R2H (lower panel). Expo-
nential fits to the peaks in the growing phase are overplotted as solid lines.
initial models exhibit a dominant m = 2 component. Assuming that the modes
behave as e−i(σmt−mϕ), the real part of σm can be obtained by Fourier trans-
forming Am. In particular Re(σ2), the bar-mode frequency, can be extracted
from either A2 or η as both represent the same mode. This is the dominant
mode in all our simulations and its frequency and growth rate are given in
Table 1. The latter corresponds to the imaginary part of σ2, which is calcu-
lated fitting an exponential to the peak values of η in the growing phase of
the evolution until the modes saturate. Other modes are also identified in the
simulations for values of Am with lower amplitudes. We have checked that
these modes are harmonics of the l = m = 2 mode so that they follow to good
accuracy the relation σm = mσp, σp being the pattern frequency, calculated
as σp = σ2/2. This is shown for model D3H in Fig. 1 which displays the spec-
trum of Am from m = 1 to m = 8 (in arbitrary units). The vertical dashed
lines in this figure indicate the location of the integer multiples of the pattern
frequency σp, their values indicated on the axis at the top of the figure. Each
spectrum for each mode is normalized to its own maximum for plotting pur-
poses. Note that the lower the mode amplitude the noisier the spectrum and
the less accurate the relation σm = mσp.
For the models of our sample subject to the “clasical” bar-mode deformation
(R1H and R2H), our simulations yield a value of β between 0.253 and 0.275,
in good agreement with the critical value for the onset of the dynamical bar-
mode instability. Model R1H is stable and model R2H is unstable. The growth
rates and frequencies reported in Table 1 agree with those of [15]. Note that
0 100 200 300
Fig. 3. Evolution of η for models D1 (upper panel), D2 (central panel) and D3 (lower
panel). Dashed lines correspond to low resolution and solid lines to high resolution.
Exponential fits to the peaks in the growing phase are overplotted as solid lines.
for model R1H, which is stable, the frequency for the m = 2 mode cannot be
computed. The time evolution of η for these two models is displayed in Fig. 2.
For the unstable model R2H, our simulations show the formation of a bar
which saturates for values of η+ and η× close to 1, i. e. in the full nonlinear
regime.
Fig. 3 shows the time evolution of η for models D in our sample, prone to suffer
the low T/|W | bar-mode instability. Solid lines correspond to high resolution
simulations and dashed lines to low resolution. For all three models the pattern
frequencies σp are such that there exists a corotation radius inside the star,
i.e. a radius at which the bar-mode rotates with the same angular velocity as
the fluid. The location of the corotation radius for all models of our sample
is reported in Table 1. As recently discussed by [22] the existence of such
corotation radius is a potential requirement for the ocurrence of the instability.
As becomes clear from Fig. 3, all models are unstable but grid resolution has
an important effect on the saturation of the instability once the nonlinear
phase has been reached, as well as in the long-term dynamics of the stars.
In the linear phase of models D1H and D2H, the growth rates and frequencies
agree with the results of [15] in both, the numerical simulations and the linear
analysis (see Table 1). In the linear phase of model D3H, our frequencies are
similar to the numerical results of [15], although our growth rates are about
a factor two larger. We emphasize that no results are reported in the linear
analysis for this model in the work of [15], and therefore this discrepancy can
be an effect of the resolution used or of the characteristics of each numerical
code. Increasing resolution leads to similar results in the frequencies but to
higher growth rates.
In the nonlinear phase, models D1 and D3 behave similarly for the two resolu-
tions used (see Fig. 3), and also similarly to the results by [15] (compare with
Fig. 3 of that paper). For model D2 we observe a radical change of behavior in
the nonlinear phase of the mode evolution depending on the grid resolution.
This has implications on the long-term dynamics of the star and, in particu-
lar, on the attainable amplitudes of the gravitational radiation emitted, as we
discuss below.
It is worth mentioning the possibility that the unstable mode at the start of
model D2H might excite some other mode in the corotation band, which could
not otherwise be excited for lower grid resolution. As discussed by [23,24] in
their study of differentially rotating shells, there are many zero-step modes in
the band, so that the whole continuous spectrum could potentially be excited.
In such case these modes would have very slow power-law growth.
For all our models we have checked mass conservation along the evolution. The
worst results are obtained for model D3H, for which mass is conserved within
2.5% error when the instability saturates. At the end of the simulation (after
48 orbital periods and 25000 iterations in the coarsest grid) the error has grown
to only 6%. For all other models mass conservation is even more accurate. Note
that these errors are within the round-off error of the code, and it is not related
to the conservation properties of the numerical scheme itself. For a regular grid
with 1283 cells and a simulation employing 25000 iterations, the accumulated
round-off error (binomial distribution) using single-precision arithmetics, is
about
1283 × 25000 × 10−8 = 0.0023 = 0.23%. Correspondingly, for a 2563
grid (with twice the number of iterations for the simulation) the error is about
0.9%. Taking into account that this error affects the nonlinear evolution of the
system, it is not surprising to have an error at the level of a few percent by
the end of our high resolution simulations, for all conserved quantities.
Figure 4 shows the evolution of Am for model D3 and for m ranging from
1 to 8 for our two resolutions. According to this figure, the only two modes
1e-06
0.0001
0 100 200
1e-06
0.0001
m=1, 3, 5, 7
m=1, 3, 5, 6, 7, 8
Fig. 4. Evolution of |Am| for model D3 with low resolution (top) and high resolution
(bottom). The m = 2 mode is represented with thick solid line, m = 4 with thin
solid line, m = 6 with dashed line, m = 8 with dot-dashed line, and all other odd
m with dotted lines.
relevant for the dynamics of the star are m = 2 and m = 4. All other modes
have smaller amplitudes and play no role in the dynamics. Note that for odd
m modes, the value of the integrated quantity Am, if close to zero, is extremely
sensitive to very small numerical asymmetries, which are induced by the patch
creation scheme of our AMR code. This explains the resolution differences in
the initial values for odd m modes in Fig. 4 (at t = 0 they start off at 10−8
level for the low resolution simulation), although they saturate at the same
value irrespective of the resolution.
An important diagnosis for the accuracy of the results is the location of the
center of mass during an evolution. The round-off error of the numerical code
imposes controlled errors in mass and linear momentum, which results in tiny
displacements of the center of mass. However small (one numerical cell in
our runs) this unphysical displacement may hinder the correct analysis of the
0 100 200 300
1e-06
0.0001
Fig. 5. Effects of the artificial displacement of the center of mass (of only one
numerical cell) on the time evolution of |A1| for model D2H. The thin solid line
shows a fictitiuos evolution resulting from the numerical artifact originated by the
center of mass displacement.
mode growth rates. For this reason all integrated quantities shown in Fig. 4
are computed after correcting for the displacement of the center of mass,
xnew = xold − xCM, in a post-processing stage of the data analysis. Were this
not done, a one-armed m = 1 mode would grow much faster than it should to
bring up fictitious features in the plots. This is shown for model D2H in Fig. 5.
The thick solid line in this figure corresponds to the evolution of the m = 1
mode taking into account the correction for the center of mass displacement,
while the thin solid line is the corresponding evolution of this mode without
the correction.
4.3 Gravitational waves
The growth and saturation of the instability is also imprinted on the gravita-
tional waves emitted. The gravitational waveforms h+ and h× for models D1,
D2, and D3, computed using the standard quadrupole formula, are shown in
Fig. 6. For a source of mass M located at a distance R those waveforms can
be calculated from the dimensionless waveform amplitudes a+ and a× as
h+,× = a+,×
sin2 θ
, (16)
using G = c = 1 units. The resulting chirp-like signal in all the models, partic-
ularly apparent for model D2L, indicates the presence of a bipolar distribution
of mass within the star (see Sec. 4.4).
-0.05
0 100 200 300
Fig. 6. Gravitational waves for models D1 to D3 extracted using the standard
quadrupole formula. Thick (thin) solid lines correspond to low (high) resolution.
Only the dimensionless waveform amplitude a+ is plotted.
As mentioned before, the effects of grid resolution on the evolution of the
nonlinear phase of the bar-mode are imprinted on the gravitational waveforms.
Thick solid lines in Fig. 6 are the waveforms which correspond to the low-
resolution models, and thin solid lines to the high-resolution counterparts.
The evolution of η for model D3, displayed in Fig. 3, shows little deviations
with grid resolution, and this translates into very similar gravitational wave
patterns (bottom panel of Fig. 6), the differences becoming more noticeable in
the nonlinear phase following saturation (Ω0t ≥ 75). For model D1 (top panel),
the differences also become more apparent at later times during the evolution,
in good agreement with the dissimilar behaviour of the matter dynamics in
this model, as encoded in the evolution of η in Fig. 3. As happens for model
D3 the first few cycles of the gravitational waveform, when the mode is still
in the linear phase, are accurately captured for both resolutions.
The major dependence of the waveform on the grid resolution is found for
model D2. Again, the linear phase for the growth of the bar deformation is
Ω 0 t=
70.40.0
Fig. 7. Snapshots of the density, vorticity, and specific angular momentum, for
model D3H, at three representative instants of the evolution. All snapshots show
slices of the stars in the equatorial plane. Quantities are normalized as follows:
max, rew
s , and l
ϕ/(rev
s ), where v
s is the initial velocity at the surface
of the star.
accurately captured irrespective of the resolution (and agrees with the per-
turbative results of [15]). This is signalled in the perfect overlapping of both
gravitational waveforms during the first three cycles (see the middle panel of
Fig. 6). However, the different nonlinear dynamics of the bar-mode deforma-
tion for this model, shown in the middle panel of Fig. 3, is severely imprinted
on the gravitational waveform. Model D2H emits gravitational waves which
have roughly one order of magnitude smaller amplitude than those computed
for the corresponding low resolution model.
4.4 Morphology
We next describe the morphological features encountered during the evolution
of some representative models. Fig. 7 shows three snaphsots of the evolution
of model D3H for the density (top), the azimuthal component of the vorticity,
~wϕ = (∇ × ~v)ϕ (middle), and the specific angular momentum, ~l = ~r × ~v
(bottom). From left to right the snapshots correspond to the initial time (Ω0t =
0), a time when the bar-mode instability is growing (Ω0t = 33.6), and the time
when the instability saturates (Ω0t = 70.4). Only the equatorial plane of the
stars is shown in all these plots. Animations of all simulations performed are
available at www.uv.es/∼cerdupa/bars/. We note that our AMR code is able
to dynamically place patches (e. g. between 4 and 8 in the D3H model) and
evolves the system with continuous matching between patches, as exemplified
in Fig. 7.
The evolution of model D3H shows that as the m = 2 mode grows the star
develops an ellipsoidal shape which remains spinning beyond saturation. Since
the low β m = 2 mode saturates at lower values (η ∼ 0.1) than the classical
bar-mode instability (η ∼ 1), no clear bars are visible in the density plot. At
late times (Ω0t > 100) a “boxy” structure becomes apparent as the m = 4
mode has grown to almost similar amplitude as the m = 2 mode (see anima-
tions and Fig. 4). No other global features can be seen, consistent with the
fact that |Am| ≪ 1 for all modes other than m = 2 and 4. The vorticity plot
shows that the m = 2 mode at Ω0t = 33.6 adopts the form of a two-armed
spiral winding up around the central parts of the star. As the mode begins to
saturate (Ω0t = 70.4) the spirals break apart into the outer layers in a turbu-
lent flow reminiscent of the (shear) Kelvin-Helmholtz instability, and shock as
they reach the atmosphere. These trends are also visible in the specific angular
momentum plot.
The presence of a corotation radius, at r/re = 0.56 for model D3H, seems
to play a role in the growth and saturation of the instability, in agreement
with the recent findings of [25]. As the bar-mode grows, pressure waves carry
angular momentum outside the corotation radius, which is deposited in the
outer layers of the star. This excites Kelvin-Helmholtz-like instabilities in the
fluid that break the mode outside the corotation radius. When this happens the
m = 2 instability stops growing and no more angular momentum is extracted.
Figure 8 shows late-time snapshots of the equatorial plane distribution of the
density perturbation, i.e. (ρ−ρ(0))/ρ(0)max, for models D2H and D3H. The times
are chosen well inside the nonlinear and saturation phase of the instability.
This figure helps to interpret the mode dynamics and its saturation along the
lines mentioned before: During the evolution the density perturbations are
shed in waves from the center towards the outer layers of the star. At late
times, when the instability saturates, such shedding stops, and the density
Fig. 8. Snapshots of the density perturbation at the equatorial plane for models D2H
and D3H. The white solid curves indicate the location of the corotation radius. The
white dashed boxes indicate the location of the patches for model D3H.
perturbation reaches the largest values outside the corotation radius (depicted
with white solid lines in Fig. 8), for either model.
We note in passing that the corotation radius in all our high resolution models
lies well inside the outer boundary of the finest box set up by the AMR
refinement pattern. (see, e.g. the white dashed boxes depicted in the right
panel of Fig. 8 indicating the location of the AMR patches for model D3H)
This rules out the possibility of a numerical artifact resulting from the patch
creation scheme of our AMR code being the cause for the different long-term
evolution between low and high resolution models, particularly noticeable for
model D2 in Fig. 3.
Finally, Fig. 9 shows a comparison between models D2L and D2H at Ω0t = 101
(i.e. well within the nonlinear phase), to highlight the effects of the numer-
ical resolution on the morphology. From top to bottom this panel shows a
schlieren plot (|∇ log ρ|) , ~wϕ, and ~l. The resolution differences in the evolu-
tion of model D2 become apparent from this figure. In particular, the “boxy”
structure becomes much more clearly visible in the low resolution simulation
(D2L), indicating an excessive growth rate of the m = 4 mode. The presence
of pressure waves is emphasized in the schlieren plot, very accurately captured
in model D2H. Those waves, once the flow is driven to turbulence past the
corotation radius, redistribute the angular momentum in the outer layers of
model D2L in a much more pronounced way than for model D2H.
Ω 0 t=
D2L D2H
Fig. 9. Resolution comparison between models D2L and D2H once the instability
has saturated. Only slices of the stars in the equatorial plane are shown.
5 Summary and outlook
We have presented AMR high-resolution simulations of the low T/|W | bar-
mode instability of extremely differentially rotating neutron stars. Our main
motivation has been to revisit the simulations by [15] on such instability,
assessing how sensitive the onset and development of the instability is to
numerical issues such as grid resolution. We have addressed the importance
of a correct treatment of delicate numerical aspects which may spoil three-
dimensional simulations in (Cartesian) grid-based codes, always hampered by
insufficient resolution, namely the handling of the low-density atmosphere sur-
rounding the star, the correction for the center of mass displacement, and the
mass and momentum conservation properties of the numerical scheme. Our
simulations have revealed the complex morphological features involved in the
nonlinear dynamics of the instability. We have found that in the nonlinear
phase of the evolution, the excitation of Kelvin-Helmholtz-like fluid modes
outside the corotation radii of the stellar models leads to the saturation of the
bar-mode deformation. While the overall trends reported in the investigation
of [15] are confirmed by our work, the resolution used to perform the simu-
lations may play a key role on the long-term behaviour of the instability and
on the nonlinear dynamics of rotating stars, which has only become apparent
for some specific models of our sample (namely model D2). This, in turn, has
implications on the attainable amplitudes of the associated gravitational wave
signals.
The work reported in this paper is a first step in our ongoing efforts of study-
ing the dynamical bar-mode instability within the magnetized core collapse
scenario.
Acknowledgements
The authors thank Harry Dimmelmeier, Nick Stergioulas, and Anna Wats
for useful comments. Research supported by the Spanish Ministerio de Edu-
cación y Ciencia (MEC; grants AYA2004-08067-C03-01, AYA2003-08739-C02-
02, AYA2006-02570). VQ is a Ramón y Cajal Fellow of the Spanish MEC.
Computations performed at the Servei d’Informática de la Universitat de
València (CERCA-CESAR).
References
[1] N. Stergioulas, Liv. Rev. Relativ. 6 (2003) 3
[2] J. E. Tohline, R. H. Durisen, & M. McCollough, ApJ 298 (1985) 220
[3] J. L. Houser, J. M. Centrella, & S. Smith, Phys. Rev. Lett. 72 (1994) 1314
[4] K. C. B. New, J. M. Centrella, & J. E. Tohline, Phys. Rev. D 62 (2000) 064019
[5] M. Shibata, T. W. Baumgarte, & S. L. Shapiro, ApJ 542 (2000) 453
[6] M. Rampp, E. Müller, & M. Ruffert, A&A 332 (1998) 969
[7] M. Shibata, & Y. Sekiguchi, Phys. Rev. D 71 (2005) 024014
[8] M. Saijo, Phys. Rev. D 71 (2005) 104038
[9] C. D. Ott, S. Ou, J. E. Tohline, & A. Burrows, ApJ 625 (2005) L119
[10] H. C. Spruit & E. S. Phinney, Nature 393 (1998) 139
[11] A. Heger, S. E. Woosley, & H. C. Spruit, ApJ 626 (2005) 350
[12] C. D. Ott, A. Burrows, T. A. Thompson, E. Livne, & R. Walder, ApJS 164
(2006) 130
[13] S. E. Woosley, & A. Heger, ApJ 637 (2006) 914
[14] J. M. Centrella, K. C. B. New, L. L. Lowe, & J. D. Brown, ApJ 550 (2001)
[15] M. Shibata, S. Karino, & Y. Eriguchi, MNRAS 334 (2002) L27
[16] M. Shibata, S. Karino, & Y. Eriguchi, MNRAS 343 (2003) 619
[17] V. Quilis, MNRAS 352 (2004) 1426
[18] M. J. Berger, P. Colella, J. Comp. Phys. 82 (1989) 64
[19] J. A. Font, M. Miller, W.-M. Suen, & M. Tobias, Phys. Rev. D 61 (2000) 0044011
[20] M. D. Duez, P. Marronetti, S. L. Shapiro, & T. W. Baumgarte, Phys. Rev. D
67 (2003) 024004
[21] Y. Eriguchi, & E. Müller, A&A, 147 (1984) 161
[22] A. L. Watts, N. Andersson, & D. I. Jones, ApJ 618 (2005) L37
[23] A. L. Watts, N. Andersson, H. Beyer, & B. F. Schutz, MNRAS 342 (2003) 1156
[24] A. L. Watts, N. Andersson, & R. L. Williams, MNRAS 350 (2004) 927
[25] M. Saijo & S. Yoshida, MNRAS 368 (2006) 1429
	Introduction
	Mathematical framework
	Numerical approach
	Results
	Initial data
	Stability analysis
	Gravitational waves
	Morphology
	Summary and outlook
	References
