Binary systems

6.3 Binary systems

As seen in Section 6.2, not many groups using spectral methods are able to solve all the three-dimensional Einstein equations in a stable way. When dealing with black holes, the situation is even worse. Only very recently, the Caltech/Cornell group succeeded in following 16 orbits, merger and ring-down of an equal-mass nonspinning black-hole–binary system [186

]. Moreover, we can report on three recent partial results in the field using spectral methods, dealing with each type of binary system (neutron stars and/or black holes) and leaving space for future study in this rapidly evolving field. We note, of course, that successful numerical evolutions of such systems have been performed with other numerical methods, which we very briefly summarize here. The first successful fully-relativistic simulation of neutron star binaries was obtained by Shibata et al. [197 , 193 ] and now, more groups are also able to study such systems: the Louisiana State University (LSU) group [4 ] and the Albert Einstein Institute (AEI, Golm) group [16 ]. We should also mention here the more microphysically-detailed simulations by Oechslin and Janka [162 ], although with the conformally-flat approximation, and those of Liu et al. [141 ] evolving magnetized neutron star binaries. Shibata and Uryū [198 , 199 ] have successfully evolved black-hole–neutron-star binaries using the puncture technique for the modeling of the black hole. As far as black hole binary systems are concerned, after many years of hard work and codes evolving the binary system for a restricted time, a first stable simulation up to the merger phase has been performed by Pretorius [176

], who used general harmonic coordinates together with constraint-damping terms and a compactification of spatial infinity. He also used the excision technique for a region surrounding the singularity inside the horizon. This first success was followed a few moths later by the Texas/Brownsville group [55 ] and the NASA/Goddard group [18 ], using very different techniques, namely BSSN with moving punctures and “1+log” slicing together with a “ $Γ$ -driver” shift condition. These techniques have rapidly become standards for many other groups, which are now able to stably evolve black hole binaries, such as the AEI/LSU collaboration [174 ], the group in Jena [93 ], that at Pennsylvania State University [114 ] and at Florida Atlantic University [220 ]. The results have reached a high level of confidence and it is possible to compare gravitational waveforms obtained from numerical evolution with post-Newtonian template families [166 ].

6.3.1 Neutron star binaries

Numerical simulations of the final stage of inspiral and merger of neutron star binaries have been performed by Faber et al. [75 ], who used spectral methods in spherical coordinates (based on Lorene library [99 ]) to solve Einstein’s equations in the conformally-flat approximation (see Sections 5 and 6.1.1). The hydrodynamic evolution has been computed using a Lagrangian smoothed particle hydrodynamics (SPH) code. As for the initial conditions described in Section 5.5, the equations for the gravitational field reduce, in the case of the conformally-flat approximation, to a set of five nonlinear coupled elliptic (Poisson-type) PDEs. The considered fields (lapse, shift and conformal factor) are “split” into two parts, each component being associated with one of the stars in the binary. Although this splitting is not unique, the result obtained is independent from it, because the equations with the complete fields are solved iteratively, for each timestep. Boundary conditions are imposed on each solution of the field equations at radial infinity thanks to a multidomain decomposition and a $u = 1∕r$ compactification in the last domain. The authors used $∼$ 10⁵ SPH particles for each run, with an estimated accuracy level of 1 – 2%. Most of the CPU time was spent in calculating the values of quantities known from their spectral representation, at SPH particle positions. Another difficulty has been the determination of the domain boundary containing each neutron star, avoiding any Gibbs phenomena. Because the conformally-flat approximation discards gravitational waves, the dissipative effects of gravitational radiation back reaction were added by hand. The authors used the slow-motion approximation [227 ] to induce a shrinking of the binary systems, and the gravitational waves were calculated with the lowest-order quadrupole formulae. The code has passed many tests and, in particular, it has evolved several quasiequilibrium binary configurations without adding the radiation reaction force with resulting orbits that were nearly circular (change in separation lower than 4%). The code was thus able to follow irrotational neutron star binaries, including radiation reaction effects, up to the merger and the formation of a differentially rotating remnant, which is stable against immediate gravitational collapse for reasonably stiff equations of state. All the results agreed pretty well with previous relativistic calculations.

6.3.2 Black-hole–neutron-star binaries

A similar combination of numerical techniques has been used by Faber et al. [74 ] to compute the dynamic evolution of merging black-hole–neutron-star binaries. In addition to the conformally-flat approximation and similar to Taniguchi et al. [210 ], Faber et al. [74 ] considered only the case of an extremely large mass ratio between the black hole and the neutron star, thus holding the black hole position fixed and restricting the spectral computational grid to the neighborhood of the neutron star. The metric describing the space surrounding the black hole was thus, supposed to keep the form of a Schwarzschild black hole in isotropic coordinates. The neutron star was restricted to low compactness (only a few percent) in order to have systems that disrupt well outside the last stable orbit. The system was considered to be in corotation and, just as for neutron star binaries, the gravitational radiation reaction was added by hand. As stated above, the numerical methods used SPH discretization to treat dynamic evolution of matter, and the spectral library Lorene [99 ] to solve the Einstein field Poisson-like equations in the conformally-flat approximation. But here, the spectral domains associated with the neutron star did not extend to radial infinity (no compactified domain) and approximate boundary conditions were imposed, using a multipole expansion of the fields. The main reason being that the black hole central singularity could not be well described on the neutron star grid.

Faber et al. [74 ] have studied the evolution of neutron-star–black-hole binaries with different polytropic indices for the neutron star matter equation of state, the initial data being obtained as solutions of the conformal thin-sandwich decomposition of Einstein’s equations. They found that, at least for some systems, the mass transfer from the neutron star to the black hole plays an important role in the dynamics of the system. For most of these systems, the onset of tidal disruption occurred outside the last stable orbit, contrary to what had been previously proposed in analytical models. Moreover, they have not found any evidence for significant shocks within the body of the neutron star. This star possibly expanded during the mass loss, eventually losing mass outward and inward provided that it was not too far within the last stable orbit. Although the major part of released matter remained bound to the black hole, a significant fraction could be ejected with sufficient velocity to become unbound from the binary system.

6.3.3 Black hole binaries

Encouraging results concerning black-hole–binary simulations with spectral methods have been first obtained by Scheel et al. [189 ]. They have used two coordinate frames to describe the motion of black holes in the spectral grid. Indeed, when using excision techniques (punctures are not regular enough to be well represented by spectral methods), excision boundaries are fixed on the numerical grid. This can cause severe problems when, due to the movement of the black hole, the excision surface becomes timelike and the whole evolution problem becomes ill-posed in the absence of boundary conditions. One solution seems to be the use of comoving coordinates, but the authors report that the GH formulation they use appear to be unstable with this setting. They, therefore, consider a first system of inertial coordinates (with respect to spatial infinity) to define the tensor components in the triad associated with these coordinates, and a second system of comoving (in some sense) coordinates. In the case of their black-hole–binary tests [189 ], they define the comoving coordinates dynamically, with a feedback control system that adjusts the moving coordinate frame to control the location of each apparent horizon center.

The spectral code uses 44 domains of different types (spherical and cylindrical shells, rectangular blocks) to describe the binary system. Most of the numerical strategy to integrate Einstein’s equations is taken from the tests on the GH formulation of Lindblom et al. [137 ] and has already been detailed in Section 6.2.1. The important technical ingredient detailed by Scheel et al. [189 ] is the particular filtering of tensor fields in terms of spherical harmonics. The dual-coordinate-frame representation can mix the tensor’s spherical harmonic components. So, in their filtering of the highest-order tensor spherical-harmonic coefficients, Scheel et al. [189 ] had to take into account this mixing by transforming spatial tensors into a rotating-frame tensor spherical-harmonic basis before filtering and then transforming back to an inertial frame basis. This method allowed them to evolve black-hole–binary spacetimes for more than four orbits, until $t ≳ 600MADM$ .

However, a central problem has been the capability of the code to follow the merger phase, and even though the code was able to compute the inspiral quite accurately, it used to fail just before the black holes merged. The problem was that, when the black holes came close to each other, their horizons became extremely distorted and strong gradients would develop in the dynamic fields. This has been explained as a gauge effect, coming from the incapacity of the gauge condition to react and change the geometry when the two black holes begin to interact strongly, and can be seen as a coordinate singularity developing close to the merger. Nevertheless, a cure has been found, as explained in Scheel et al. [186 ]. The original gauge is kept until some given time and then smoothly changed to a new one, based on the gauge treatment by Pretorius [177 , 176 ] (for the lapse): the gauge source function is evolved through a damped, driven wave equation, which drives the lapse toward unity and the shift vector toward zero near the horizons. Thus, the horizons expand in coordinate space and the dynamic fields are smoothed out near the horizons, preventing gauge singularities from developing. With this transition of the gauge condition, the evolution of the black holes can be tracked until the formation of a common horizon encompassing both black holes. Then, the evolution of this single-distorted dynamic black hole is achieved by first interpolating all variables onto a new computational domain containing only one excised region, then by choosing a new comoving coordinate system, and finally by again modifying the gauge condition to drive the shift vector to a time-independent state.

These new gauge conditions have allowed Scheel et al. [186 ] to follow the inspiral during 16 orbits, and the merger and ring-down phase of an equal-mass nonspinning black-hole–binary system. They were able to compute the mass and the spin of the final black hole with very high accuracy (10^–5 and 10^–4 relative accuracy for the mass and spin, respectively), and to extract the physical waveform accurately to 0.01 radians in phase. This is the first spectral numerical simulation of the full evolution of a black-hole–binary system.