Hydrodynamical evolutions of polytropic models of spherical neutron stars can be used as test-bed computations for multi-dimensional codes. Representative examples are the simulations by [111], with pseudo-spectral methods, and by [244] with HRSC schemes. These investigations adopted radial-gauge polar-slicing coordinates in which the general relativistic equations are expressed in a simple way that resembles Newtonian hydrodynamics. Gourgoulhon [111] used a numerical code to detect, dynamically, the zero value of the fundamental mode of a neutron star against radial oscillations. Romero et al. [244] highlighted the accuracy of HRSC schemes by finding, numerically, a change in the stability behavior of two slightly different initial neutron star models: For a given EOS, a model with mass is stable and a model of is unstable. More recently, in [274] a method based on the nonlinear evolution of deviations from a background stationary-equilibrium star was applied to study nonlinear radial oscillations of a neutron star. The accuracy of the approach permitted a detailed investigation of nonlinear features such as quadratic and higher order mode-coupling and nonlinear transfer of energy.

Axisymmetric pulsations of rotating neutron stars can be excited in several scenarios, such as core collapse, crust and core quakes, and binary mergers, and they could become detectable either via gravitational waves or high-energy radiation. An observational detection of such pulsations would yield valuable information about the EOS of relativistic stars. As a first step towards the study of pulsations of rapidly rotating relativistic stars, Font, Stergioulas, and Kokkotas [97] developed an axisymmetric numerical code that integrates the equations of general relativistic hydrodynamics in a fixed background spacetime. The finite difference code is based on a state-of-the-art approximate Riemann solver [70] and incorporates different second- and third-order TVD and ENO numerical schemes. This code is capable of accurately evolving rapidly rotating stars for many rotational periods, even for stars at the mass-shedding limit. The test simulations reported in [97] showed that, for non-rotating stars, small amplitude oscillations have frequencies that agree to better than 1% with linear normal-mode frequencies (radial and non-radial) in the so-called Cowling approximation (i.e., when the evolution of the spacetime variables is neglected). Axisymmetric modes of pulsating non-rotating stars are computed in [269], both in Cowling and fully coupled evolutions. Contrary to the 2+1 approach followed by [97], the code used in [269] evolves the relativistic stars on null spacetime foliations (see Section 2.2.2).

Until very recently (see below), the quasi-radial modes of rotating relativistic stars had been studied only under simplifying assumptions, such as in the slow-rotation approximation or in the relativistic Cowling approximation. An example of the latter can be found in [88], where a comprehensive study of all low-order axisymmetric modes of uniformly and rapidly rotating relativistic stars was presented, using the code developed by [97]. The frequencies of quasi-radial and non-radial modes with spherical harmonic indices and 3 were computed through Fourier transforms of the time evolution of selected fluid variables. This was done for a sequence of appropriately perturbed stationary rotating stars, from the non-rotating limit to the mass-shedding limit. The frequencies of the axisymmetric modes are affected significantly by rotation only when the rotation rate exceeds about 50% of the maximum allowed. As expected, at large rotation rates, apparent mode crossings between different modes appear.

In [89], the first mode frequencies of uniformly rotating stars in full
general relativity and rapid rotation were obtained, using the
three-dimensional code
`GR_ASTRO`

described in Section
3.3.2
. Such frequencies were computed both in fixed spacetime
evolutions (Cowling approximation) and in coupled hydrodynamical
and spacetime evolutions. The simulations used a sequence of
(perturbed)
*N*
=1,
*K*
=100 () polytropes of central density
, in which the rotation rate
varies from zero to 97% of the maximum allowed rotational
frequency,
. The Cowling runs allowed a comparison with earlier results
reported by [88], obtaining agreement at the 0.5% level. The fundamental
mode-frequencies and their first overtones obtained in fully
coupled evolutions show a dependence on the increased rotation
which is similar to the one observed for the corresponding
frequencies in the Cowling approximation [88].

Relativistic hydrodynamical simulations of nonlinear
*r*
-modes are presented in [279] (see also [155] for Newtonian simulations). The gravitational radiation
reaction-driven instability of the
*r*
-modes might have important astrophysical implications, provided
that the instability does not saturate at low amplitudes by
nonlinear effects or by dissipative mechanisms. Using a version
of the
`GR_ASTRO`

code, Stergioulas and Font [279] found evidence that the maximum
*r*
-mode amplitude in isentropic stars is of order unity. The
dissipative mechanisms proposed by different authors to limit the
mode amplitude include shear and bulk viscosity, energy loss to a
magnetic field driven by differential rotation, shock waves, or
the slow leak of the
*r*
-mode energy into some short wavelength oscillation modes
(see [16] and references therein). The latter mechanism would
dramatically limit the
*r*
-mode amplitude to a small value, much smaller than those found
in the simulations of [279,
155] (see [278] for a complete list of references on the subject). Energy leak
of the
*r*
-mode into other fluid modes has been recently considered
by [113] through Newtonian hydrodynamical simulations, finding a
catastrophic decay of the amplitude only once it has grown to a
value larger than that reported by [16].

The bar mode instability in differentially rotating stars in general relativity has been analyzed by [261] by means of 3+1 hydrodynamical simulations. Using the code discussed in Section 3.3.1, Shibata et al. [261] found that the critical ratio of rotational kinetic energy to gravitational binding energy for compact stars with is , slightly below the Newtonian value for incompressible Maclaurin spheroids. All unstable stars are found to form bars on a dynamical timescale.

The accurate simulation of a binary neutron star coalescence is, however, one of the most challenging tasks in numerical relativity. These scenarios involve strong gravitational fields, matter motion with (ultra-)relativistic speeds, and relativistic shock waves. The numerical difficulties are exacerbated by the intrinsic multi-dimensional character of the problem and by the inherent complexities in Einstein's theory of gravity, such as coordinate degrees of freedom and the possible formation of curvature singularities (e.g., collapse of matter configurations to black holes). It is thus not surprising that most of the (few) available simulations have been attempted in the Newtonian (and post-Newtonian) framework (see [235] for a review). Many of these studies employ Lagrangian particle methods such as SPH, and only a few have considered (less viscous) high-order finite difference methods such as PPM [246].

Concerning relativistic simulations, Wilson's formulation of the hydrodynamic equations (see Section 2.1.2) was used in Refs. [303, 304, 169]. Such investigations assumed a conformally flat 3-metric, which reduces the (hyperbolic) gravitational field equations to a coupled set of elliptic (Poisson-like) equations for the lapse function, the shift vector, and the conformal factor. These early simulations revealed the unexpected appearance of a ``binary-induced collapse instability'' of the neutron stars, prior to the eventual collapse of the final merged object. This effect was reduced, but not eliminated fully, in revised simulations [169], after Flanagan [85] pointed out an error in the momentum constraint equation as implemented by Wilson and coworkers [303, 304]. A summary of this controversy can be found in [235]. Subsequent numerical simulations with the full set of Einstein equations (see below) did not find this effect.

Nakamura and coworkers have been pursuing a programme to simulate neutron star binary coalescence in general relativity since the late 1980's (see, e.g., [196]). The group developed a three-dimensional code that solves the full set of Einstein equations and self-gravitating matter fields [214]. The equations are finite-differenced in a uniform Cartesian grid using van Leer's scheme [290] with TVD flux limiters. Shock waves are spread out using a tensor artificial viscosity algorithm. The hydrodynamic equations follow Wilson's Eulerian formulation and the ADM formalism is adopted for the Einstein equations. This code has been tested by the study of the gravitational collapse of a rotating polytrope to a black hole (comparing to the axisymmetric computation of Stark and Piran [276]). Further work to achieve long term stability in simulations of neutron star binary coalescence is under way [214]. We note that the hydrodynamics part of this code is at the basis of Shibata's code (Section 3.3.1), which has successfully been applied to simulate the binary coalescence problem (see below).

The head-on collision of two neutron stars (a limiting case of
a coalescence event) was considered by Miller et al. [183], who performed time-dependent relativistic simulations using
the code described in Section
3.3.2
. These simulations analyzed whether the collapse of the final
object occurs in prompt timescales (a few milliseconds) or
delayed (after neutrino cooling) timescales (a few seconds).
In [254] it was argued that in a head-on collision event, sufficient
thermal pressure is generated to support the remnant in
quasi-static equilibrium against (prompt) collapse prior to slow
cooling via neutrino emission (delayed collapse). In [183], prompt collapse to a black hole was found in the head-on
collision of two
neutron stars modeled by a polytropic EOS with
and
. The stars, initially separated by a proper distance of
, were boosted toward one another at a speed of
(the Newtonian infall velocity). The simulation employed a
Cartesian grid of
points. The time evolution of this simulation can be followed in
the movie in Figure
11
. This animation simultaneously shows the rest-mass density and
the internal energy evolution during the on-axis collision. The
formation of the black hole in prompt timescales is signalled by
the sudden appearance of the apparent horizon at
(*t*
=63.194 in code units). The violet dotted circles indicate the
trapped photons. The animation also shows a moderately
relativistic shock wave (Lorentz factor
) appearing at
(code units; yellow-white colors), which eventually is followed
by two opposite moving shocks (along the infalling
*z*
direction) that propagate along the atmosphere surrounding the
black hole.

The most advanced simulations of neutron star coalescence in full general relativity are those performed by Shibata and Uryu [258, 266, 267]. Their numerical code, briefly described in Section 3.3.1, allows the long-term simulation of the coalescences of both irrotational and corotational binaries, from the innermost stable circular orbit up to the formation and ringdown of the final collapsed object (either a black hole or a stable neutron star). Their code also includes an apparent horizon finder, and can extract the gravitational waveforms emitted in the collisions. Shibata and Uryu have performed simulations for a large sample of parameters of the binary system, such as the compactness of the (equal mass) neutron stars (), the adiabatic index of the -law EOS (), and the maximum density, rest mass, gravitational mass, and total angular momentum. The initial data correspond to quasi-equilibrium states, either corotational or irrotational, the latter being more realistic from considerations of viscous versus gravitational radiation timescales. These initial data are obtained by solving the Einstein constraint equations and the equations for the gauge variables under the assumption of a conformally flat 3-metric and the existence of a helical Killing vector (see [267] for a detailed explanation). The binaries are chosen at the innermost orbits for which the Lagrange points appear at the inner edge of the neutron stars, and the plunge is induced by reducing the initial angular momentum by .

The comprehensive parameter space survey carried out by [258, 266, 267] shows that the final outcome of the merger depends sensitively on the initial compactness of the neutron stars before plunge. Hence, depending on the stiffness of the EOS, controlled through the value of , if the total rest mass of the system is times larger than the maximum rest mass of a spherical star in isolation, the end product is a black hole. Otherwise, a marginally-stable massive neutron star forms. In the latter case the star is supported against self-gravity by rapid differential rotation. The star could eventually collapse to a black hole once sufficient angular momentum has dissipated via neutrino emission or gravitational radiation. The different outcome of the merger is reflected in the gravitational waveforms [267]. Therefore, future detection of high-frequency gravitational waves could help to constrain the maximum allowed mass of neutron stars. In addition, for prompt black hole formation, a disk orbiting around the black hole forms, with a mass less than 1% the total rest mass of the system. Disk formation during binary neutron star coalescence, a fundamental issue for cosmological models of short duration GRBs, is enhanced for unequal mass neutron stars, in which the less massive one is tidally disrupted during the evolution (Shibata, private communication).

A representative example of one of the models simulated by
Shibata and Uryu is shown in Figure
12
. This figure is taken from [267]. The compactness of each star in isolation is
*M*
/
*R*
=0.14 and
. Additional properties of the initial model can be found in
Table 1 of [267]. The figure shows nine snapshots of density isocontours and the
velocity field in the equatorial plane (*z*
=0) of the computational domain. At the end of the simulation a
black hole has formed, as indicated by the thick solid circle in
the final snapshot, representing the apparent horizon. The
formation timescale of the black hole is larger the smaller the
initial compactness of each star. The snapshots depicted in
Figure
12
show that once the stars have merged, the object starts
oscillating quasi-radially before the complete collapse takes
place, the lapse function approaching zero
non-monotonically [267]. The collapse toward a black hole sets in after the angular
momentum of the merged object is dissipated through gravitational
radiation. Animations of various simulations (including this
example) can be found at Shibata's website [257].

To close this section we mention the work of Duez et al. [72] where, through analytic modelling of the inspiral of corotational and irrotational relativistic binary neutron stars as a sequence of quasi-equilibrium configurations, the gravitational wave-train from the last phase (a few hundred cycles) of the inspiral is computed. These authors further show a practical procedure to construct the entire wave-train through coalescence by matching the late inspiral waveform to the one obtained by fully relativistic hydrodynamical simulations as those discussed in the above paragraphs [266]. Detailed theoretical waveforms of the inspiral and plunge similar to those reported by [72] are crucial to enhance the chances of experimental detection in conjunction with matched-filtering techniques.

Alternative: single figures.

Numerical Hydrodynamics in General Relativity
José A. Font
http://www.livingreviews.org/lrr-2003-4
© Max-Planck-Gesellschaft. ISSN 1433-8351 Problems/Comments to livrev@aei-potsdam.mpg.de |