5 Additional Information4 Hydrodynamical Simulations in Relativistic 4.2 Accretion onto black holes

4.3 Hydrodynamical evolution of neutron stars 

The numerical investigation of many interesting astrophysical processes involving neutron stars, such as the rotational evolution of proto-neutron stars (which can be affected by a dynamical bar mode instability and by the Chandrasekhar-Friedman-Schutz instability) or the gravitational radiation from unstable pulsation modes or, more importantly, from the catastrophic coalescence and merger of neutron star compact binaries, requires the ability of accurate, long-term hydrodynamical evolutions employing relativistic gravity. These scenarios are receiving increasing attention in recent years [263, 168, 96Jump To The Next Citation Point In The Article, 183Jump To The Next Citation Point In The Article, 258Jump To The Next Citation Point In The Article, 262, 266Jump To The Next Citation Point In The Article, 97Jump To The Next Citation Point In The Article, 89Jump To The Next Citation Point In The Article, 260Jump To The Next Citation Point In The Article].

4.3.1 Pulsations of relativistic stars 

The use of relativistic hydrodynamical codes to study the stability properties of neutron stars and to compute mode frequencies of oscillations of such objects has increased in recent years (see, e.g., the Living Reviews article by Stergioulas [278Jump To The Next Citation Point In The Article] and references therein). An obvious advantage of the hydrodynamical approach is that it includes, by construction, nonlinear effects, which are important in situations where the linearized equations (commonly employed in calculations of mode-frequencies of pulsating stars) break down. In addition, due to the current status of hydrodynamics codes, the computation of mode frequencies in rapidly rotating relativistic stars might be easier to achieve through nonlinear numerical evolutions than by using perturbative computations (see, e.g., the results in [89Jump To The Next Citation Point In The Article, 260]).

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 [111Jump To The Next Citation Point In The Article], with pseudo-spectral methods, and by [244Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline5448 is stable and a model of tex2html_wrap_inline5450 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 [97Jump To The Next Citation Point In The Article] 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 [97Jump To The Next Citation Point In The Article] 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 [269Jump To The Next Citation Point In The Article], both in Cowling and fully coupled evolutions. Contrary to the 2+1 approach followed by [97Jump To The Next Citation Point In The Article], 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 [88Jump To The Next Citation Point In The Article], 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 tex2html_wrap_inline5452 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 (tex2html_wrap_inline5282) polytropes of central density tex2html_wrap_inline5462, in which the rotation rate tex2html_wrap_inline4980 varies from zero to 97% of the maximum allowed rotational frequency, tex2html_wrap_inline5466 . The Cowling runs allowed a comparison with earlier results reported by [88Jump To The Next Citation Point In The Article], 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 [279Jump To The Next Citation Point In The Article] (see also [155Jump To The Next Citation Point In The Article] 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 [279Jump To The Next Citation Point In The Article] 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 [16Jump To The Next Citation Point In The Article] 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 [261Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline5480 is tex2html_wrap_inline5482, slightly below the Newtonian value tex2html_wrap_inline5484 for incompressible Maclaurin spheroids. All unstable stars are found to form bars on a dynamical timescale.

4.3.2 Binary neutron star coalescence 

Many of the current efforts in code development in relativistic astrophysics aim to simulate the coalescence of compact binaries. Neutron star binaries are among the most promising sources of gravitational radiation to be detected by the various ground-based interferometers worldwide. The computation of the gravitational waveform during the most dynamical phase of the coalescence and plunge depends crucially on hydrodynamical, finite-size effects. This phase begins once the stars, initially in quasi-equilibrium orbits of gradually smaller orbital radius (due to the emission of gravitational waves) reach the so-called innermost stable circular orbit. About tex2html_wrap_inline5486 years after formation of the binary system, the gravitational wave frequency enters the LIGO/VIRGO high frequency band. The final plunge of the two objects takes place on a dynamical timescale of a few ms. During the last 15 minutes or so until the stars finally merge, the frequency is inside the LIGO/VIRGO sensitivity range. About 16,000 cycles of waveform oscillation can be monitored, while the frequency gradually shifts from tex2html_wrap_inline5488 to tex2html_wrap_inline5490 . A perturbative treatment of the gravitational radiation in the quadrupole approximation is valid as long as tex2html_wrap_inline5492 and tex2html_wrap_inline5494 simultaneously, M being the total mass of the binary, R the neutron star radius, and r the separation of the two stars. As the stars approach each other and merge, both inequalities are less valid and therefore, fully relativistic hydrodynamical calculations become necessary.

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 [235Jump To The Next Citation Point In The Article] 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. [303Jump To The Next Citation Point In The Article, 304Jump To The Next Citation Point In The Article, 169Jump To The Next Citation Point In The Article]. 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.


Click on thumbnail to view movie

Figure 11: Movie showing the animation of a head-on collision simulation of two tex2html_wrap_inline4658 neutron stars obtained with a relativistic code [96, 183Jump To The Next Citation Point In The Article]. The movie shows the evolution of the density and internal energy. The formation of the black hole in prompt timescales is demonstrated by the sudden appearance of the apparent horizon at tex2html_wrap_inline4660 ( t=63.194 in code units), which is indicated by violet dotted circles representing the trapped photons. See [28] for download options of higher quality versions of this movie.

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 [214Jump To The Next Citation Point In The Article]. 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. [183Jump To The Next Citation Point In The Article], 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 tex2html_wrap_inline4658 neutron stars modeled by a polytropic EOS with tex2html_wrap_inline5234 and tex2html_wrap_inline5512 . The stars, initially separated by a proper distance of tex2html_wrap_inline5514, were boosted toward one another at a speed of tex2html_wrap_inline5516 (the Newtonian infall velocity). The simulation employed a Cartesian grid of tex2html_wrap_inline5518 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 tex2html_wrap_inline4660 (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 tex2html_wrap_inline5524) appearing at tex2html_wrap_inline5526 (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 [258Jump To The Next Citation Point In The Article, 266Jump To The Next Citation Point In The Article, 267Jump To The Next Citation Point In The Article]. 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 (tex2html_wrap_inline5530), the adiabatic index of the tex2html_wrap_inline4726 -law EOS (tex2html_wrap_inline5534), 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 [267Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline5536 .

The comprehensive parameter space survey carried out by [258, 266Jump To The Next Citation Point In The Article, 267Jump To The Next Citation Point In The Article] 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 tex2html_wrap_inline4726, if the total rest mass of the system is tex2html_wrap_inline5540 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 [267Jump To The Next Citation Point In The Article]. 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 [267Jump To The Next Citation Point In The Article]. The compactness of each star in isolation is M / R =0.14 and tex2html_wrap_inline5544 . Additional properties of the initial model can be found in Table 1 of [267Jump To The Next Citation Point In The Article]. 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 [267Jump To The Next Citation Point In The Article]. 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. [72Jump To The Next Citation Point In The Article] 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.


Click on thumbnail to view image

Alternative: single figures.


Click on thumbnail to view imageClick on thumbnail to view imageClick on thumbnail to view image
Click on thumbnail to view imageClick on thumbnail to view imageClick on thumbnail to view image
Click on thumbnail to view imageClick on thumbnail to view imageClick on thumbnail to view image

Figure 12: Snapshots of the density contours in the equatorial plane for a binary neutron star coalescence that leads to a rotating black hole (see [267Jump To The Next Citation Point In The Article] for the characteristics of the model). Vectors indicate the local velocity field tex2html_wrap_inline4664 . tex2html_wrap_inline4666 denotes the orbital period of the initial configuration. The thick solid circle in the final panel indicates the apparent horizon. The figure is taken from [267] (used with permission).

5 Additional Information4 Hydrodynamical Simulations in Relativistic 4.2 Accretion onto black holes

image Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to livrev@aei-potsdam.mpg.de