There are two distinct schemes used in all binary merger calculations performed to date, the BSSN (Baumgarte–Shapiro–Shibata–Nakamura) [277, 27] and Generalized Harmonic formalisms. For general reviews of these formalisms, as well as other developments in numerical relativity, we refer readers to two recent books on numerical relativity [4, 30]. Here we only briefly summarize these two schemes.
The BSSN formalism was adapted from the 3+1 ADM approach, with quantities defined as in Eqs. 7 and 8. While the original ADM scheme proved to be numerically unstable, defining the auxiliary quantities and treating these expressions as independent variables stabilized the system and allowed for long-term evolutions. While slight variants exist, the 19 evolved variables are typically either the conformal factor or its logarithm , the conformal 3-metric , the trace of the extrinsic curvature, the trace free extrinsic curvature and the conformal connection functions . The evolution equations themselves are given in Appendix A.
The BSSN scheme was used in the binary merger calculations of the KT collaboration [287, 288, 285, 286], the first completely successful NS-NS calculations ever performed in full GR. Ever since the UTB/RIT  and Goddard  groups showed simultaneously that a careful choice of gauge allows BHs to be evolved in BSSN schemes without the need to excise the singularity, these “puncture gauges” have gained widespread hold, and have been used to evolve NS-NS binaries (and in some cases, BH-NS binaries) by the KT collaboration , UIUC , and Whisky .
The generalized harmonic formalism, developed over about two decades from initial theoretical suggestions up to its current numerical implementation [112, 115, 232, 130, 231] was used to perform the first calculations of merging BH-BH binaries by Pretorius , and has since been applied to NS-NS binaries by the HAD collaboration  and to BH-NS mergers by HAD , SXS [85, 84, 108], and the Princeton group [294, 88]. It involves introducing a set of auxiliary quantities denoted representing the action of the wave operator on the spacetime coordinates themselves[171, 6] evolve equations for the spacetime metric along with its spatial derivative, and projected time derivative , subject to consistency constraints on the spatial derivatives. The first BH-NS merger calculations in the GH formalism used a first-order reduction  of the Einstein equations and specified the source functions to damp to zero exponentially in time, while the first binary NS merger work  used a similar first-order reduction and chose harmonic coordinates with .
In both formalisms, most groups employ grid-based finite differencing to evaluate spatial derivatives. While finite differencing operators may be easily written down to arbitrary orders of accuracy, there is a trade-off between the computational efficiency achievable by using smaller, second-order stencils and the higher accuracy that can be attained using larger, higher-order ones. For the moment, many groups are now moving to at least fourth-order accurate differencing techniques, with a high likelihood that at least the field sector of NS-NS merger calculations will soon be performed at comparable order to BH-BH calculations, at either sixth  or eighth-order  accuracy, if not higher. The main limitation to date involves the complexity of shock capturing using higher-order schemes, as we discuss in Section 5.2.3 below.
Imposing physically realistic and accurate boundary conditions remains a difficult task for numerical codes. Ideally, one wishes to impose boundary conditions at large distances that preserve the GR constraints and yield a well-posed initial-boundary value problem. On physical grounds, the boundary should only permit outgoing waves, preventing the reflection of spurious waves back into the numerical grid. Otherwise, reflections may be avoided by placing the outer boundaries so far away from the matter that they remain causally disconnected from the merging objects for the full duration of the simulation. Building upon previous work (see, e.g., [113, 151, 12, 150, 256, 242], Winicour  derived boundary conditions that satisfy all of the desired conditions for the generalized harmonic formalism. No such conditions have been derived for the full BSSN formalism, though progress has been made (see, e.g., [43, 129]) so that we may now construct well-posed boundary conditions in the weak-gravity linearized limit of BSSN  and for related first-order gravitational formalisms .
While unigrid schemes are extremely convenient, they tend to be extremely inefficient, since one must resolve small-scale features in the central regions of a merger but also extend grids out to the point where the GW signal has taken on its roughly asymptotic form. Thus, nearly every code incorporates some means of focusing resolution on the high-density regions via some form of mesh refinement. A simple approach, for instance, is to use “fisheye” coordinates that represent a continuous radial deformation of a grid , a technique that had previously been used successfully, e.g., for BH-BH mergers [21, 62].
While fixed mesh refinement offers the chance for greater computational efficiency and accuracy, much current work focuses on adaptive mesh refinement, in which the level of refinement of the grid is allowed to evolve dynamically to react to the evolving binary configuration. Several codes, both public and private, now implement this functionality. The publicly available Carpet computational toolkit [263, 262], which is distributed to the community as part of the open source Einstein Toolkit [90, 175] uses a “moving boxes” approach, and has been designed to be compatible with the widely used and publicly available Cactus framework. It has been successfully implemented by the Whisky group  to perform NS-NS mergers, by UIUC for their BH-NS mergers , and a host of other groups for BH-BH mergers and additional problems. The KT code SACRA also implements an adaptive mesh refinement (AMR) scheme for NS-NS and BH-NS mergers , as does the most recent version of the HAD collaboration’s code , which is based on the publicly available infrastructure of the same name [169, 8], and the BAM code employed by the Jena group [308, 41, 122]. The Princeton group also has an AMR code, which has been used to perform BH-NS mergers to date [294, 88, 89]
One drawback of employing rectangular grids is that memory costs scale like , where is the number of grid cells across a side, and total computational effort like once one accounts for the reduction in the timestep due to the Courant–Friedrich–Levy criterion. Since one does not necessarily need high angular resolution at large radii, there is great current interest in developing schemes that use some form of spheroidal grid, for which the memory scaling is merely . A group at LSU has implemented a multi-patch method , in which space is broken up into a number of non-overlapping domains in such a way that the six outermost regions (projections of the faces of a cube onto spheres of constant radius), maintain constant angular resolution and thus produce linear dependence of the total number of grid points on the number of radial points. To date, it has been used primarily for vacuum spacetimes  and hydrodynamics on a fixed background. The SXS collaboration, begun at Caltech and Cornell and now including members at CITA and Washington State, has used a spectral evolution code with multiple domains to evolve BH-NS binaries, which achieves the same scaling by expanding the metric fields in modes rather than in position space . Their first published results on NS-NS binaries are currently in preparation (see  for a brief summary of work to date).
While all of these grid techniques produce tremendous advantages in computational efficiency, each required careful study since deformations of a grid or the introduction of multiple domains can introduce inaccuracies and non-conservative effects. As an example, in AMR schemes, one must deal with the same reflection issues at refinement boundaries that are present at the physical boundaries of the grid, as discussed above, though the interior nature of the boundaries allows for additional techniques, such as tapered grid boundaries , to be used to minimize reflections there. The study of how to minimize spurious effects in these schemes continues, and will represent an important topic for years to come, especially as evolution schemes become more complicated by including more physical effects.
Fluids cannot be treated in the same way as the spacetime metric because finite differencing operators do not return meaningful results when placed across discontinuities induced by shocks. Instead, GR(M)HD calculations must include some form of shock-capturing that accounts for these jumps. These are typically implemented in conservative schemes, in which the fluid variables are transformed from the standard “primitive” set , which includes the fluid density, pressure, and velocity (and magnetic field in MHD evolutions), into a new set so that the evolution equations may be written in the formflux functions and source terms can be expressed in terms of the primitive variables but not their derivatives. These schemes allow one to evolve the resulting MHD set of equations so that numerical fluxes are conserved to numerical precision across cell walls as the fluid evolves in time. One widely used scheme, often referred to as the Valencia formulation , is described in Appendix B..
There are important mathematical reasons for casting the GRHD/GRMHD system in conservative form, primarily since the mathematical techniques describing Godunov methods may be called into play . In such methods, we assume that the evolution of the quantities may be expressed in the form.
First one reconstructs the primitives from the conserved variables on both sides of an interface, using interpolation schemes designed to respect the presence of shocks. Simple schemes involving limiters yield second-order accuracy in general but first-order accuracy at shocks, while higher-order methods such as PPM (piecewise parabolic method) and essentially non-oscillatory (ENO) schemes such as CENO (central ENO) and WENO (weighted ENO) yield higher accuracy but at much higher computational cost. Once primitives are interpolated to the cell interfaces, flux terms are evaluated there and one solves the Riemann problem describing the evolution of two distinct fluid configurations placed on either side of a membrane (see  for a description). While complete solutions of the Riemann problem are painstaking to evolve, a number of approximation techniques exist and do not reduce the order of accuracy of the code. Finally, one must take the conservative solution, now advanced forward in time, and recover the underlying primitive variables, a process that requires solving as many as eight simultaneous equations in the case of GRMHD or five for GRHD systems. A number of methods to do this have been widely studied , and simplifying techniques are known for specific cases (for the case of polytropic EOSs in GRHD evolutions, one need only invert a single non-analytic expression and the remaining variables can then be derived analytically).
The inclusion of magnetic fields in hydrodynamic calculations adds another layer of complexity beyond shock capturing. Magnetic fields must be evolved in such a way that they remain divergence-free, much in the same way that relativistic evolutions must satisfy the Hamiltonian and momentum constraints. Brute force attempts to subtract away any spurious divergence often lead to instabilities, so more intricate schemes have been developed. “Divergence cleaning” schemes typically introduce a new field representing the magnetic field divergence and use parabolic/hyperbolic equations to damp the divergence away while moving it off the computational domain; the approach is relatively simple to implement but prone to small-scale numerical errors . “Constrained transport” schemes stagger the grids on which different physical terms are calculated to enforce the constraints (see, e.g.,  for a particular implementation), and have been applied widely to many different physical configurations. Recently, a new scheme in which the vector potential is used rather than the magnetic field was introduced by Etienne, Liu, and Shapiro [93, 95], and found to yield successful results for a variety of physical configurations including NSs and BHs.
Living Rev. Relativity 15, (2012), 8
This work is licensed under a Creative Commons License.