5 Additional Information4 Simulations of Astrophysical Phenomena4.2 Accretion onto black holes

4.3 Hydrodynamical evolution of neutron stars 

The numerical evolution of neutron stars in general relativity is nowadays a very active field of research [206, 130, 75Jump To The Next Citation Point In The Article, 140Jump To The Next Citation Point In The Article, 202Jump To The Next Citation Point In The Article, 203Jump To The Next Citation Point In The Article, 205Jump To The Next Citation Point In The Article, 76]. The numerical investigation of many interesting astrophysical processes involving neutron stars, such as the rotational evolution of proto-neutron stars or the gravitational radiation from unstable pulsation modes or, more importantly, from the catastrophic coalescence and merger of neutron star-neutron star or black hole-neutron star binaries, requires the ability to perform accurate long-term hydrodynamical evolutions employing relativistic gravity.

Evolutions of polytropic models of spherical neutron stars (i.e. with EOS tex2html_wrap_inline4091, 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 [87Jump To The Next Citation Point In The Article], employing pseudo-spectral methods, and by [184Jump 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 which resembles Newtonian hydrodynamics. Gourgoulhon [87] used a numerical code to detect, dynamically, the zero value of the fundamental mode of a neutron star against radial oscillations. Romero et al. [184] 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 tex2html_wrap_inline4095 is stable and a model of tex2html_wrap_inline4097 is unstable.

Three-dimensional hydrodynamical evolutions of relativistic, self-gravitating neutron stars have been considered in [75Jump To The Next Citation Point In The Article, 202Jump To The Next Citation Point In The Article, 8Jump To The Next Citation Point In The Article, 203Jump To The Next Citation Point In The Article, 205Jump To The Next Citation Point In The Article]. In [75Jump To The Next Citation Point In The Article] 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 [183] and Marquina's flux formula [57]. The Einstein equations were formulated using two different approaches: (i) the standard ADM formalism and (ii) a hyperbolic formulation developed in [35]. 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 [75Jump To The Next Citation Point In The Article] 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 [140Jump To The Next Citation Point In The Article] and is briefly described below.

The simulations presented in [202Jump To The Next Citation Point In The Article, 203Jump To The Next Citation Point In The Article, 205Jump To The Next Citation Point In The Article, 8Jump To The Next Citation Point In The Article] are performed using a reformulation of the ADM equations, first proposed by Shibata and Nakamura [204Jump To The Next Citation Point In The Article], and recently taken up by Baumgarte and Shapiro [23Jump To The Next Citation Point In The Article]. 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 [23Jump To The Next Citation Point In The Article]. This formulation has proven to be much more long-term stable (though somehow less accurate, see [8Jump To The Next Citation Point In The Article]) than standard ADM. The tests analyzed include evolutions of weak gravitational waves [23Jump To The Next Citation Point In The Article], self-gravitating matter configurations with prescribed (analytic) time evolution [22], spherical dust collapse [202Jump To The Next Citation Point In The Article], strong gravitational waves (which even collapse to black holes) [8Jump To The Next Citation Point In The Article], black holes [8Jump To The Next Citation Point In The Article], and boson stars [8]. Hydrodynamical evolutions of neutron stars are extensively considered in [202Jump To The Next Citation Point In The Article, 203, 205Jump To The Next Citation Point In The Article].

  

Click on thumbnail to view image

Figure 10: Schematic illustration showing the evolution of a neutron star binary (left side) together with the numerical approach (right side) best suited for an accurate description of each portion of the evolution. General relativistic hydrodynamic simulations are unavoidable during the most dynamical phase from the innermost stable circular orbit to the final merge.

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 tex2html_wrap_inline4101 and tex2html_wrap_inline4103 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 [186]. The interested reader is referred to the recent review by Rasio and Shapiro [180Jump To The Next Citation Point In The Article] and an upcoming Living Reviews article by Swesty [213], 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 [230Jump To The Next Citation Point In The Article, 231Jump To The Next Citation Point In The Article] 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 [180] for an updated discussion). In particular, a radial stability analysis carried out by [21] 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 [69Jump To The Next Citation Point In The Article] 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 [131], which incorporate the correction identified by Flanagan [69], 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., [150]). A three-dimensional code employing the full set of Einstein equations and self-gravitating matter fields has been developed [164Jump To The Next Citation Point In The Article]. In this code the complete set of equations, spacetime and hydrodynamics, are finite-differenced on a uniform Cartesian grid using van Leer's scheme [219Jump To The Next Citation Point In The Article] 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 [211]) and applied to the coalescence of a binary neutron star system. Further work to achieve long term stability is currently under way [164Jump To The Next Citation Point In The Article].

The most advanced simulations of neutron star coalescence in full general relativity are those reported recently by Shibata [202Jump To The Next Citation Point In The Article, 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 [164]: The hydrodynamic equations are formulated following Wilson's approach and they are solved using van Leer's [219] 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 [204] and recently slightly modified by [23]. In [202] Shibata computed the merger of two models of corotating binary neutron stars of tex2html_wrap_inline3493 in contact and in approximate quasi-equilibrium orbits. The central density of each star in the two models is tex2html_wrap_inline4115 (mildly relativistic) and tex2html_wrap_inline4117 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 tex2html_wrap_inline4123 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 tex2html_wrap_inline4125, which is nearly the maximum allowed density along the sequence of stable neutron stars of K =10 and tex2html_wrap_inline3493 . 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 (tex2html_wrap_inline4131, J being the angular momentum and tex2html_wrap_inline4135 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.

  

Click on thumbnail to view image

Figure 11: Animation of a head-on collision simulation of two tex2html_wrap_inline3541 neutron stars obtained with a relativistic code [75Jump To The Next Citation Point In The Article, 140Jump 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_inline4139 ( t=63.194 in code units), which is indicated by violet dotted circles representing the trapped photons. See  [241] for download options of higher quality versions of this movie.

Recently, Miller et al. [140Jump To The Next Citation Point In The Article] have studied the head-on collision of two neutron stars by means of time-dependent relativistic simulations using the code of Font et al. [75]. 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 [195] 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 [140] prompt collapse to a black hole was found in the head-on collision of two tex2html_wrap_inline3541 neutron stars modeled by a polytropic EOS with tex2html_wrap_inline4145 and tex2html_wrap_inline4147 . The stars are initially separated by a proper distance of tex2html_wrap_inline4149 and are boosted towards one another at a speed of tex2html_wrap_inline4151 (the Newtonian infall velocity). The simulation employed a Cartesian grid of tex2html_wrap_inline4153 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 tex2html_wrap_inline4155 (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 tex2html_wrap_inline4159 (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.



5 Additional Information4 Simulations of Astrophysical Phenomena4.2 Accretion onto black holes

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