5.2 GR numerical techniques

5.2.1 GR formalisms and gauge choice

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. 7View Equation and 8View Equation. While the original ADM scheme proved to be numerically unstable, defining the auxiliary quantities &tidle;Γ i = − &tidle;γi,jj 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 &tidle;γij, the trace K of the extrinsic curvature, the trace free extrinsic curvature Aij and the conformal connection functions &tidle;Γ i. The evolution equations themselves are given in Appendix A.

The BSSN scheme was used in the binary merger calculations of the KT collaboration [287Jump To The Next Citation Point, 288Jump To The Next Citation Point, 285Jump To The Next Citation Point, 286Jump To The Next Citation Point], the first completely successful NS-NS calculations ever performed in full GR. Ever since the UTB/RIT [61Jump To The Next Citation Point] and Goddard [22Jump To The Next Citation Point] 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 [332Jump To The Next Citation Point], UIUC [172Jump To The Next Citation Point], and Whisky [17Jump To The Next Citation Point].

The generalized harmonic formalism, developed over about two decades from initial theoretical suggestions up to its current numerical implementation [112, 115, 232, 130Jump To The Next Citation Point, 231Jump To The Next Citation Point] was used to perform the first calculations of merging BH-BH binaries by Pretorius [231Jump To The Next Citation Point], and has since been applied to NS-NS binaries by the HAD collaboration [6Jump To The Next Citation Point] and to BH-NS mergers by HAD [66Jump To The Next Citation Point], SXS [85Jump To The Next Citation Point, 84Jump To The Next Citation Point, 108Jump To The Next Citation Point], and the Princeton group [294Jump To The Next Citation Point, 88Jump To The Next Citation Point]. It involves introducing a set of auxiliary quantities denoted μ H representing the action of the wave operator on the spacetime coordinates themselves

1 √ --- H μ ≡ Γ μ ≡ gαβΓ μαβ = − √----∂ν( − gg μν) = g αβ∇ α∇ βxμ = □ xμ, (25 ) − g
which are treated as independent gauge variables whose evolution equation must be specified. Current first-order formulations [171Jump To The Next Citation Point, 6Jump To The Next Citation Point] evolve equations for the spacetime metric gμν along with its spatial derivative, Φiμν = ∂igμν and projected time derivative α Πμν = − n ∂αgμν, subject to consistency constraints on the spatial derivatives. The first BH-NS merger calculations in the GH formalism used a first-order reduction [171] of the Einstein equations and specified the source functions to damp to zero exponentially in time, while the first binary NS merger work [6Jump To The Next Citation Point] used a similar first-order reduction and chose harmonic coordinates with μ H = 0.

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 [136] or eighth-order [179] 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 [329] 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, 129Jump To The Next Citation Point]) so that we may now construct well-posed boundary conditions in the weak-gravity linearized limit of BSSN [204] and for related first-order gravitational formalisms [49].

5.2.2 Grid decompositions

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 [172Jump To The Next Citation Point], 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 [90Jump To The Next Citation Point, 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 [17Jump To The Next Citation Point] to perform NS-NS mergers, by UIUC for their BH-NS mergers [94Jump To The Next Citation Point], 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 [332Jump To The Next Citation Point], as does the most recent version of the HAD collaboration’s code [6Jump To The Next Citation Point], which is based on the publicly available infrastructure of the same name [169, 8Jump To The Next Citation Point], and the BAM code employed by the Jena group [308Jump To The Next Citation Point, 41Jump To The Next Citation Point, 122Jump To The Next Citation Point]. The Princeton group also has an AMR code, which has been used to perform BH-NS mergers to date [294Jump To The Next Citation Point, 88Jump To The Next Citation Point, 89]

One drawback of employing rectangular grids is that memory costs scale like 3 N, where N is the number of grid cells across a side, and total computational effort like N 4 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 ∝ N. A group at LSU has implemented a multi-patch method [337], 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 [224] 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 [85Jump To The Next Citation Point]. Their first published results on NS-NS binaries are currently in preparation (see [141] 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 [168], 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.

5.2.3 Hydrodynamics, MHD, and high-resolution shock capturing

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 ⃗P, which includes the fluid density, pressure, and velocity (and magnetic field in MHD evolutions), into a new set ⃗U so that the evolution equations may be written in the form

∂(⃗U ) = ∇ ⋅F⃗ + ⃗S, (26 ) t
where the flux functions ⃗F(P⃗) and source terms S⃗(⃗P ) 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 [186Jump To The Next Citation Point], 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 [121]. In such methods, we assume that the evolution of the quantities ⃗ U may be expressed in the form

⃗ ⃗ Δt-( ⃗ ⃗ ) U (x = xi,t = tn+1) = U (x = xi,t = tn) + Δx F (x = xi−1∕2) − F (x = xi+1∕2) , (27 )
where the points have spatial coordinates xi ≡ x0 + iΔx and the timesteps satisfy tn = t0 + nΔt. The fluxes must be determined by solving the Riemann problem at each cell face (thus the half-integer indices), either exactly or approximately. It can be shown that solutions constructed this way, when convergent, must converge to a solution of the original problem, even when shocks are present [163].

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 [106] 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 [203], 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 [24]. “Constrained transport” schemes stagger the grids on which different physical terms are calculated to enforce the constraints (see, e.g., [310] 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.

  Go to previous page Go up Go to next page