Black-hole–binary systems

5.6 Black-hole–binary systems

5.6.1 Digging the holes

Though the computation of black hole binaries in circular orbit has a lot of common features with the neutron star case, there are also some differences that need to be addressed. In at least one aspect, black holes are much simpler objects because they are a solution of Einstein’s equations without matter. So the whole issue of investigating various equations of state is irrelevant and there is no need to solve any equation for the matter. However, there is a price to pay and one must find a way to impose the presence of black holes in the spacetime. Two ideas have been proposed.

In the puncture method, the full spacetime contains three asymptotically-flat regions. One is located at $r = ∞$ and the other two at two points, $M 1$ an $M 2$ , which are called the punctures. The presence of flat regions near the punctures is enforced by demanding that some quantities, like the conformal factor, diverge at those points (typically as $1∕r$ ). The discontinuities are taken out analytically and the equations are solved numerically for the regular parts in the whole space. This idea dates back to the work of Brill and Lindquist [50 ], at least in the case of black holes initially at rest.The puncture approach has been successfully applied to the computation of quasicircular orbits by means of spectral methods in [8 ].

The apparent horizon method relies on initial works by Misner [151 ] and Lindquist [140 ]. In this case, the space has only two asymptotically-flat regions. One can show that this is equivalent to solving Einstein’s equations outside two spheres on which boundary conditions must be imposed. The boundary conditions are based on the concept of trapped surfaces and apparent horizons. The physical state of the black holes are precisely encoded in the boundary conditions.

5.6.2 First configurations

The first configurations of black hole binaries computed by means of spectral methods can be found in [110 ]. The formalism and various hypotheses are given in the companion paper [98 ]. The assumptions are very similar to those used for neutron star binaries (see Section 5.5.1). Helical symmetry is enforced and conformal flatness assumed. The holes are described by the apparent horizon technique. However, the boundary conditions used have been shown to be only approximately valid, up to a rather good accuracy. This effect is discussed in the original paper [110 ] and further explored by Cook in [62 ]. The numerical techniques are very similar to the ones employed for neutron-star–binary configurations (see Section 5.5.2). Two sets of spherical domains are used, one for each black hole. Boundary conditions are imposed on the surface between the nucleus and the first shell. Both sets extend up to infinity using a compactification in $1∕r$ .

For the first time, good agreement was found between numerical results and post-Newtonian ones. A detailed comparison can be found in [67 ]. In particular, the location of the energy minimum is shown to coincide to within a few percent. This improvement with respect to previous numerical work is mainly due to a difference in the physical hypothesis (i.e. the use of helical symmetry). One important product of [110 ] is the use of a new criterion to determine the appropriate value of the orbital angular velocity $Ω$ . Indeed, for neutron stars, this is done by demanding that the fluid of both stars be in equilibrium [100 ]. This, of course, is not applicable for black holes. Instead, in [98 , 110 ] it is proposed that one find $Ω$ by demanding that the ADM mass and the Komar-like mass coincide. Without going into too much detail, this amounts to demanding that, far from the binary and at first order in $1∕r$ , the metric behave like the Schwarzschild. It is shown in [98 ] that it can be linked to a relativistic virial theorem. Since then it has been shown that this criterion can also be used for neutron stars [215 ] and that it is equivalent to the use of a variational principle called the effective potential method [60 , 22 , 173 , 59 ], where the binding energy is minimized with respect to $Ω$ .

5.6.3 Further investigation

More recently, two other spectral codes have been developed in the context of black hole binaries and successfully applied to address some of the issues raised by the work of [98 , 110 ].

One of these codes comes from the Caltech/Cornell group of Pfeiffer et al. and is described extensively in [172 , 168 ]. The code is multidomain and two main types of domains are used i) square domains in which each Cartesian-like coordinate is expanded in terms of Chebyshev polynomials and ii) spherical domains in which spherical harmonics are used for the angles $(𝜃,φ )$ and Chebyshev polynomials for the radial coordinate. Space can be compactified by standard use of the variable $1∕r$ . The two types of domains can be set up in various manners to accommodate the desired geometry. When using both square and spherical domains, there are regions of space that belong to more than one domain. This is to be contrasted with work by the Meudon group in which different domains are only touching but not overlapping. The algorithm of [172 ] solves differential equations by using a multidimensional collocation method. The size of the resulting system is roughly equal to the number of collocation points. It is then solved iteratively via a Newton–Raphson algorithm with a line search. At each step of the Newton–Raphson method, the linear system is solved by means of an iterative scheme (typically GMRES). This inner iterative solver requires careful preconditioning to work properly. Various tests are passed by the code in [172 ], in which elliptic equations and systems are solved in either spherical or bispherical topologies. In the cases presented, the error decays spectrally.

In [170 ] the code is used to investigate different ways of solving the constraint equations. Three different decompositions are used: the conformal TT one, the physical TT one and the thin-sandwich decomposition. When solving for the constraint equations only, one must also choose some freely specifiable variables, which describe the physical state of the system. In [170 ], these specifiable variables are fixed using a superposition of two Kerr–Schild black holes. The net result of [170 ] is that global quantities, like the total energy, are very sensitive to the choice of decomposition. The variation of total energy can be as large as 5%, which is the order of the energy released by gravitational waves. It is also shown that the choice of extrinsic curvature tensor is more crucial than the one of conformal metric, in accordance with an underlying result of [110 ]. Let us point out that the equations derived form the helical Killing vector approach in [98 , 110 ] are equivalent to the ones obtained by making use of the thin-sandwich decomposition of the constraints. The freely specifiable variables are obtained by both the imposition of the helical Killing symmetry and by solving an additional equation for the lapse function (resulting in the extended thin-sandwich formalism).

In [63 ] the boundary conditions based on the apparent horizon formalism [62 ] are implemented and tested numerically in the case of one and two black holes. In the latter case, the main difference from [110 ] lies in the use of more elaborate and better boundary conditions on the horizons of the black holes. By allowing for a nonvanishing lapse on the horizons, the authors of [63 ] solve the constraint equations exactly. This is to be contrasted with [110 ], where the momentum constraint equation is only solved up to a small correction. Both results show rather good agreement. This is not surprising as the correction used by the Meudon group was known to be small (see Figures 10 and 11 of [110 ]). More results are presented in [59 ], for both corotating and irrotational black holes. An important result of [59 ] is the comparison of the two criteria for determining the orbital angular velocity $Ω$ . They indeed show that the effective potential method first introduced in [60 ] and the method based on the virial theorem proposed in [98 ] are in very good agreement.

By slightly extending the boundary conditions used in [59 ], the authors of [169 ] propose to reduce the eccentricity of the black-hole–binary configurations. This is done by giving the black holes a small radial velocity by modifying the boundary condition on the shift vector. The code and other equations are the same as in [59 ]. Time evolution of the obtained initial data does show that this technique can reduce the eccentricity of the binary. However, the effect on the emitted gravitational wave is small and probably unimportant.

Another application of the Caltech/Cornell solver can be found in [143 ], where the emphasis is put on nearly maximum spinning black holes. Initial data are constructed for both single black holes and binaries. Three families of initial data are investigated. Using a formalism based on the Kerr–Schild spacetime, the authors are able to reach spins as large as $a = 0.9998$ . Such nearly-maximum spinning black holes may be relevant from the astrophysical point of view. Evolutions of these data are also discussed there.

The other spectral code used to compute configurations of black hole binaries comes from Ansorg [6 ]. It shares a lot of features with the code developed by the same author in the context of rotating stars [9 , 10 ] already discussed in Section 5.2.7. Space is decomposed into two domains. One of them lies just outside the horizons of the black holes and bispherical-like coordinates are used. The other domain extends up to infinity and an appropriate mapping is used in order to avoid the singularity of the bispherical coordinates at spatial infinity (see Section IV of [6 ]). The angle of the bispherical coordinates (i.e. the angle around the x-axis joining the two black holes) is expanded in terms of a Fourier series and the two other coordinates in terms of Chebyshev polynomials. As in [13 , 172 ], the partial differential equations are solved using a collocation method and the resulting system is solved by the Newton–Raphson method. Once again the linear system coming from the Jacobian is solved by an iterative scheme with preconditioning. The code is used to compute essentially the same configuration as those shown in [59 ]. An interesting point of [6 ] is the detailed investigation of the convergence of the results with increased resolution. As can be seen in Figure 4 of [6 ], the error initially decreases exponentially, but, for high number of points, it seems that the error follows only a power law. This is an indication that some non- $𝒞∞$ fields must be present. It is conjectured in [6 ] that this comes from logarithm terms that cannot be dealt with properly with a compactification in $1∕r$ . The same kind of effect is investigated in some detail in [109 ], where some criteria for the appearance of such terms are discussed.

A code very similar to the one used in [6 ] has also been used to compute spacetimes with black holes using the puncture approach [8 ]. Given that the black holes are no longer described by their horizons, one does not need to impose inner boundary conditions. The absence of this requirement enables the author of [8 ] to use a single domain to describe the whole space, from the puncture up to infinity. The other features of the spectral solver are the same as in [6 ]. This method has been successfully applied to the computation of black-hole–binary configurations in the puncture framework. The authors have, in particular, investigated high mass ratios between the bodies and compared their results with the ones given in the test-mass limit around a Schwarzschild black hole. The discrepancy is found to be on the order of 50% for the total energy. It is believed that this comes from the fact that the mass of each puncture cannot be directly related to the local black hole mass (see discussion in Section VII of [8 ]).

Finally, let us mention that the algorithms developed by Ansorg in [9 , 10 , 8 , 6 ] have all been unified in [7 ] to accommodate any type of binaries. Various domain decompositions are exhibited that can be used to represent neutron stars, excised black holes or puncture black holes, with the compactification of space. The algorithms are shown to be applicable to limiting cases such as large mass ratios.