Evolutions of polytropic models of spherical neutron stars (i.e. with EOS , K being the polytropic constant) using relativistic hydrodynamics can be used as test-bed computations for multidimensional codes. One-dimensional hydrodynamical studies of relativistic stars have been performed by , employing pseudo-spectral methods, and by  with HRSC schemes. These investigations adopted radial gauge polar slicing coordinates in which the general relativistic equations are expressed in a simple way which resembles Newtonian hydrodynamics. Gourgoulhon  used a numerical code to detect, dynamically, the zero value of the fundamental mode of a neutron star against radial oscillations. Romero et al.  highlighted the accuracy of HRSC schemes by finding, numerically, a change in the stability behavior of two slightly different initial neutron star models: A model with mass is stable and a model of is unstable.
Three-dimensional hydrodynamical evolutions of relativistic, self-gravitating neutron stars have been considered in [75, 202, 8, 203, 205]. In  a new efficient parallel general relativistic hydrodynamical code was presented and thoroughly tested. In this code the Valencia hydrodynamical (conservative) formulation was adopted and a variety of state-of-the-art Riemann solvers were implemented, including Roe's solver  and Marquina's flux formula . The Einstein equations were formulated using two different approaches: (i) the standard ADM formalism and (ii) a hyperbolic formulation developed in . The code was subjected to a series of convergence tests, with (twelve) different combinations of the spacetime and hydrodynamics finite differencing schemes, demonstrating the consistency of the discrete equations with the differential equations. The simulations performed in  include, among others, the evolution of equilibrium configurations of compact stars (solutions to the TOV equations), and the evolution of relativistically boosted TOV stars (v =0.87 c) transversing diagonally across the computational domain. In the long term plan this code is designed to simulate the inspiral coalescence of two neutron stars in general relativity. The more ``academic'' scenario of a head-on collision has been already analyzed in  and is briefly described below.
The simulations presented in [202, 203, 205, 8] are performed using a reformulation of the ADM equations, first proposed by Shibata and Nakamura , and recently taken up by Baumgarte and Shapiro . This formulation is based upon a conformal decomposition of the three metric and extrinsic curvature. Additionally, three ``conformal connection functions'' are introduced so that the principal part of the conformal Ricci tensor is an elliptic operator acting on the components of the conformal three metric. In this way the evolution equations reduce to a coupled set of non-linear, inhomogeneous wave equations for the conformal three metric, which are coupled to evolution equations for the (gauge) connection functions. Further details can be found in . This formulation has proven to be much more long-term stable (though somehow less accurate, see ) than standard ADM. The tests analyzed include evolutions of weak gravitational waves , self-gravitating matter configurations with prescribed (analytic) time evolution , spherical dust collapse , strong gravitational waves (which even collapse to black holes) , black holes , and boson stars . Hydrodynamical evolutions of neutron stars are extensively considered in [202, 203, 205].
Nowadays, most of the current efforts in developing codes in relativistic astrophysics are strongly motivated by the simulation of the coalescence of compact binaries. These scenarios are considerered the most promising sources of gravitational radiation to be detected by the planned laser interferometers going online worldwide in the next few years. The computation of the gravitational waveform during the most dynamical phase of the coalescence stage 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 innermost stable circular orbit (see the schematic plot in Fig. 10). From here on, the final merger of the two objects takes place in a dynamical timescale and lasts for a few milliseconds. A treatment of the gravitational radiation as a perturbation in the quadrupole approximation would be valid as long as and simultaneously, M being the mass of the binary, R the neutron star radius and r the separation of the two stars. As the stars gradually approach each other and merge, both inequalities are less valid and fully relativistic 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/or strong shock waves. The difficulties of a successful numerical integration are exacerbated by the intrinsic multidimensional character of the problem and by the inherent complexities in Einstein's theory of gravity, e.g. coordinate degrees of freedom and the possible formation of curvature singularities (e.g. collapse of matter configurations to black holes).
For these reasons it is not surprising that most of the available simulations have been attempted in the Newtonian (and post-Newtonian) framework. Most 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 . The interested reader is referred to the recent review by Rasio and Shapiro  and an upcoming Living Reviews article by Swesty , for detailed descriptions of the current status of Newtonian simulations.
Concerning relativistic simulations Wilson's hydrodynamical formulation has been applied to the study of neutron star binary coalescence in [230, 231] under the assumption of a conformally flat spacetime, which leads to a considerable simplification of the gravitational field equations. In this case they reduce to a coupled set of elliptic (Poisson-like) equations for the lapse function, shift vector and conformal factor. Their simulations revealed, unexpectedly, the appearance of a ``binary-induced collapse instability'' of the neutron stars, prior to the eventual collapse of the final merged object, with the central density of each star increasing by an amount proportional to 1/ r . The numerical results of Wilson and collaborators have received considerable attention in the literature, and their unexpected outcome has been strongly criticized by many authors on theoretical grounds (see references in  for an updated discussion). In particular, a radial stability analysis carried out by  showed that fully relativistic, corotating binary configurations are stable against collapse to black holes all the way down to the innermost stable circular orbit. More recently Flanagan  has pointed out the use of an incorrect form of the momentum constraint equation in the simulations performed by Wilson and collaborators [230, 231], which gives rise to a first post-Newtonian-order error in the scheme, showing analytically that this error can cause the observed increase of the central densities obtained in the simulations. However, in revised hydrodynamical simulations performed by Mathews and Wilson , which incorporate the correction identified by Flanagan , it was found that the compression effect is not completely eliminated, although its magnitude significantly diminishes at a given angular momentum. Reliable numerical simulations with the full set of Einstein equations will ultimately clarify these results.
Nakamura and co-workers have been pursueing a programme to simulate neutron star binary coalescence in general relativity since the late 80's (see, e.g., ). A three-dimensional code employing the full set of Einstein equations and self-gravitating matter fields has been developed . In this code the complete set of equations, spacetime and hydrodynamics, are finite-differenced on a uniform Cartesian grid using van Leer's scheme  with TVD flux limiters. Shock waves are spread out using a tensor artificial viscosity algorithm. The hydrodynamic equations follow Wilson's 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 ) and applied to the coalescence of a binary neutron star system. Further work to achieve long term stability is currently under way .
The most advanced simulations of neutron star coalescence in full general relativity are those reported recently by Shibata [202, 205]. His code is able to simulate a coalescence event for a long time from the innermost circular orbit up to the formation of a final merged object (either a black hole or a neutron star). The code shares many features with that of Oohara and Nakamura : The hydrodynamic equations are formulated following Wilson's approach and they are solved using van Leer's  second order finite difference scheme with artificial viscosity. The most important difference concerns the reformulation of the ADM Einstein equations into a conformal traceless system, as mentioned previously. This formulation was originally introduced by Shibata and Nakamura  and recently slightly modified by . In  Shibata computed the merger of two models of corotating binary neutron stars of in contact and in approximate quasi-equilibrium orbits. The central density of each star in the two models is (mildly relativistic) and in geometrized units. The quasi-equilibrium models are constructed assuming a polytropic EOS with K =10. For both initial models, the neutron stars begin to merge forming spiral arms at half the orbital period P of the quasi-equilibrium states, and by the final object is a rapidly and differentially rotating highly flattened neutron star. For the more relativistic model the central density of the merged object is , which is nearly the maximum allowed density along the sequence of stable neutron stars of K =10 and . Hence, a black hole could be formed directly in the merger of initially more massive neutron stars. For the cases considered by Shibata, the new star is strongly supported by rapid rotation (, J being the angular momentum and the gravitational mass) and could eventually collapse to a black hole once sufficient angular momentum has dissipated through, e.g., neutrino emission or gravitational radiation.
Recently, Miller et al.  have studied the head-on collision of two neutron stars by means of time-dependent relativistic simulations using the code of Font et al. . These simulations are aimed at investigating whether the collapse of the final object occurs in prompt timescales (i.e. a few milliseconds) or delayed (after neutrino cooling) timescales (i.e. a few seconds). In  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). Nevertheless, in  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 are initially separated by a proper distance of and are boosted towards 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 Quicktime movie in Fig. 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 demonstrated 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 of about 1.2) appearing at (code units; yellow-white colors) which eventually is followed by two opposite moving shocks (along the infalling z direction) which propagate along the atmosphere surrounding the black hole.
|Numerical Hydrodynamics in General Relativity
José A. Font
© Max-Planck-Gesellschaft. ISSN 1433-8351
Problems/Comments to firstname.lastname@example.org