Numerical-Relativity Waveforms

• Numerical-relativity waveforms are computed by directly solving Einstein’s equations for dynamic systems like merging black holes and neutron stars.
• They employ advanced methods such as BSSN and generalized harmonic formulations combined with extraction techniques like extrapolation and CCE to minimize errors.
• Hybridization with analytic models and inclusion of higher harmonics enhance the precision and applicability of these waveforms for gravitational-wave detection and tests of general relativity.

Numerical-relativity waveforms are gravitational-wave signals computed through direct numerical integration of Einstein’s equations for strongly gravitating, dynamical systems—such as merging black hole and neutron star binaries. These waveforms, often used as the “ground truth” in gravitational-wave astronomy, are critical for constructing templates, calibrating analytic and semi-analytic models, and for direct parameter estimation and testing of general relativity in the strong-field regime. Their development requires careful treatment of initial data, evolution, waveform extraction, and post-processing, as well as rigorous error budgeting to ensure sufficient accuracy for advanced ground- and space-based detectors.

1. Waveform Generation and Physical Content

The majority of numerical-relativity (NR) simulations of compact binaries are based on formulations such as the BSSN (Baumgarte–Shapiro–Shibata–Nakamura) or generalized harmonic (GH) systems, typically implemented in codes such as SpEC, Maya, BAM, or the Einstein Toolkit (Hinder et al., 2013, Ferguson et al., 2023). Initial conditions are set via puncture or excision methods and often employ 3.5PN-corrected orbital parameters for minimal eccentricity (Habib et al., 2020). For black hole binaries, the relevant intrinsic parameters span mass ratio, spin vectors, and initial orbital eccentricity; for neutron stars, additional parameters such as the equation of state and tidal deformability are included (Bernuzzi et al., 2011, Kiuchi et al., 2017).

The direct output from NR simulations is usually the Weyl scalar $\psi_4$, a complex combination of the Riemann curvature encoding outgoing gravitational radiation. The measured or detector-relevant waveform, the strain $h = h_+ - i h_\times$, is related to $\psi_4$ by

$$h(t) = \int_{-\infty}^t dt' \int_{-\infty}^{t'} dt''\, \psi_4(t'')$$

necessitating a double time integration (see Section 3).

Waveforms are decomposed as

$$h(t,\theta,\phi) = \sum_{\ell \ge 2} \sum_{m=-\ell}^{\ell} h^{\ell m}(t)\, {}_{-2}Y_{\ell m}(\theta,\phi),$$

where ${}_{-2}Y_{\ell m}$ are spin-weight $-2$ spherical harmonics. Contemporary catalogs include not only the dominant $(\ell,|m|)=(2,2)$ mode but also higher-order modes up to $\ell=4$ or greater (Joshi et al., 2022, Ferguson et al., 2023).

2. Post-Processing: Integration and Extraction to Infinity

Waveform extraction is performed on spheres of finite coordinate radii and then extrapolated to null infinity, usually via fitting amplitude and phase as functions of $1/r$ (Chu et al., 2015). Two principal techniques are:

• Extrapolation: The waveform at radius $r$ is fit as

$$h_{\ell m}(t_\mathrm{ret}, r) = \sum_{k=0}^{n} h_{\ell m}^{(k)}(t_\mathrm{ret}) \left(\frac{1}{r}\right)^k,$$

with $h^{(0)}_{\ell m}(t)$ representing the asymptotic waveform (Chu et al., 2015).

• Cauchy-Characteristic Extraction (CCE): The data from a finite-radius 'worldtube' is supplied as inner boundary data for a characteristic evolution to $\mathcal{I}^+$ (future null infinity), giving $h(u,\theta,\phi)$ directly in the Bondi-Sachs framework (Chu et al., 2015, Mitman et al., 2021).

Because simulations output $\psi_4$, robust integration is essential. Standard time- or frequency-domain techniques introduce secular nonlinear drifts due to spectral leakage and amplification of low-frequency noise (Reisswig et al., 2010). The Fixed Frequency Integration (FFI) prescription

$$\widetilde{h}(\omega) = \frac{\widetilde{\psi}_4(\omega)}{(-\omega^2)}$$

with a low-frequency cutoff

$$\widetilde{h}(\omega) = \frac{\widetilde{\psi}_4(\omega)}{-\omega_0^2} \quad \text{for } |\omega| < \omega_0$$

mitigates these artifacts effectively.

The effect of the extraction method is non-negligible: errors from waveform extraction, finite-length windowing, and numerical truncation all contribute comparably, with extraction sometimes dominating at a noise-weighted mismatch level of $\sim 3 \times 10^{-4}$ (Chu et al., 2015).

3. Frame and Gauge Alignment

Physical comparison and hybridization of waveforms require consistent choice of asymptotic frames. In asymptotically flat spacetimes, the Bondi-Metzner-Sachs (BMS) group governs the coordinate freedom at $\mathcal{I}^+$ (Mitman et al., 2021). The BMS frame includes translations, boosts, rotations, and supertranslations (angle-dependent time shifts), all of which can affect mode mixing and physical interpretation.

Traditional Newtonian center-of-mass corrections do not generally suffice, especially in high mass-ratio or precessing systems. Instead, asymptotic Poincaré charges calculated from the waveform data provide accurate identification of the center-of-mass and boost frame, using

$G_{(\Psi)} = \frac{E_{(\Psi)}}{P^t}$

where $E_{(\Psi)}$ and $P^t$ are energy-moment and energy charges derived from the radiative multipoles (Mitman et al., 2021).

After fixing the Poincaré frame, remaining supertranslation freedom is resolved by mapping to a post-Newtonian Bondi frame through minimization of the $L^2$ norm of the difference between NR and 3PN PN strain waveforms over the early inspiral:

$$\|f(u,\theta,\phi)\|^2 = \int_{u_1}^{u_1+T} \int_{S^2} f(u,\theta,\phi)\, \overline{f(u,\theta,\phi)}\, d\Omega\, du$$

with $$f(u,\theta,\phi) = h^{\mathrm{NR}}(u,\theta,\phi) - h^{\mathrm{PN}}(u,\theta,\phi)$$ (Mitman et al., 2021). This ensures consistent hybridization, dramatically reducing spurious mode mixing and physical inaccuracies in the hybrid waveform.

4. Characterization of Precision and Application to Detectors

Waveform accuracy is quantified by the overlap (or mismatch) computed in a detector’s noise-weighted inner product:

$$\langle h_1, h_2 \rangle = 4 \, \mathrm{Re} \int_{f_{\mathrm{low}}}^{f_{\mathrm{high}}} \frac{\tilde{h}_1(f)\, \tilde{h}_2^*(f)}{S_n(f)}\,df$$

with $S_n(f)$ the noise power spectral density (Hinder et al., 2013, Wang et al., 27 Jan 2024). Mismatch $\mathcal{M} = 1 - FF$ with $$FF=\max_{t,\phi} \langle h_1, h_2\rangle/(\|h_1\|\|h_2\|)$$ is required to be typically $<1\%$ for detection-level analyses (Wang et al., 27 Jan 2024).

Studies find that state-of-the-art NR waveforms uniformly achieve accuracy $\gtrsim 99\%$ (mismatch $<~1\%$) when compared over the frequency bands relevant for ground-based (LIGO) and space-based (LISA, Taiji, Tianqin) detectors. The limiting accuracy is set by numerical truncation, extraction, and finite-length errors; improvement to $<3 \times 10^{-4}$ in mismatch is routinely reached for state-of-the-art aligned-spin binaries (Chu et al., 2015, Wang et al., 27 Jan 2024).

5. Hybridization, Parameter-Space Coverage, and Model Calibration

Direct NR simulations are computationally expensive and typically cover only the late inspiral through ringdown (spanning tens to hundreds of cycles). For complete inspiral-merger-ringdown (IMR) signals required by matched-filter searches, NR is smoothly hybridized with post-Newtonian or EOB models (Ajith et al., 2012, Hannam et al., 2010, Kumar et al., 2013, Szilagyi et al., 2015). The hybridization procedure involves:

• Matching the phases and amplitudes (via least-squares minimization) over a window in time or frequency,
• Blending the waveforms using a window function (e.g., cosine-squared, Planck-taper),
• Rigorous verification through overlaps and direct time/frequency-domain inspection.

Hybrid waveform length requirements are driven by post-Newtonian error budgets; for example, to achieve $\mathcal{M}<3\%$ in nonspinning BBH, 5–10 pre-merger GW cycles of NR data suffice for $q\leq4$ and moderate spins, but up to 20 or more cycles may be required for high spins or more extreme mass ratios (Hannam et al., 2010).

Templates built either directly from NR (using stochastic placement in parameter space) or via hybridization with analytic models produce robust template banks effectual for detection and parameter estimation across the mass and spin parameter space (Kumar et al., 2013). Error analysis and cross-comparisons using metrics such as the match integral and root-mean-square phase errors underpin the validation and calibration of semi-analytic waveform models.

6. Higher-Order Modes, Eccentricity, and Advanced Source Modeling

Until recently, most NR waveform catalogs focused on nonprecessing and quasicircular systems, reporting only the dominant $(2,2)$ mode. With the improved reach of gravitational wave detectors and the increasing importance of parameter estimation, NR simulations now routinely include higher modes up to $(\ell,|m|)=4$, eccentric and precessing binaries, and extreme mass ratios (Joshi et al., 2022, Ferguson et al., 2023). Inclusion of higher harmonics enhances the signal-to-noise by 3.5–35%, significantly affecting the detectability and parameter recovery of eccentric and high-mass ratio mergers (Joshi et al., 2022).

Recent comparisons of PN and NR waveforms in eccentric binaries highlight good agreement in the inspiral for $22$ modes but increasing phase differences and amplitude residuals as merger approaches, or as eccentricity and mass ratio increase; these findings constrain the parameter-regions in which analytic models remain valid (Wang et al., 26 Sep 2024).

Advanced catalogs (e.g., the second MAYA catalog) systematically capture high mass ratios, precession, and eccentricity, enabling robust calibration and validation of next-generation waveform models intended for LIGO-Virgo-KAGRA, LISA, and future observatories (Ferguson et al., 2023).

7. Future Directions and Best Practices

Continued improvements in initial data (with accurate eccentricity reduction (Habib et al., 2020)), error assessment, and extraction methods expand the range and reliability of NR waveforms. Developments in BMS frame fixing (Mitman et al., 2021) are establishing best practices needed for surrogate modeling and accurate hybridization. Increasingly, NR waveforms inform effective-one-body and phenomenological models via systematic injection of NR results, even for dynamical captures and highly eccentric orbits (Andrade et al., 2023).

Public waveform catalogs with machine-readable formats (e.g., mayawaves Python library (Ferguson et al., 2023)), validated convergence, and clear documentation are accelerating integration with data analysis pipelines and fostering reproducibility. As detector sensitivities increase and parameter space coverage broadens, NR waveforms will play a central role in gravitational-wave astrophysics and tests of general relativity.

---

This article has synthesized key aspects including waveform generation and extraction, frame and gauge ambiguities, methods for post-processing and hybridization, practical accuracy metrics for detector applications, and the incorporation of higher modes and eccentricity. Each section is grounded in recent and foundational literature, reflecting the rigorous standards and evolving best practices in the field.