Accuracy Issues for Numerical Waveforms

We study the convergence properties of our implementation of the ‘moving punctures’ approach at very high resolutions for an equal-mass, non-spinning, black-hole binary. We find convergence of the Hamiltonian constraint on the horizons and the L2 norm of the Hamiltonian constraint in the bulk for sixth and eighth-order finite difference implementations. The momentum constraint is more sensitive, and its L2 norm shows clear convergence for a system with consistent sixth-order finite differencing, while the momentum and BSSN constraints on the horizons show convergence for both sixth and eighth-order systems. We analyze the gravitational waveform error from the late-inspiral, merger, and ringdown. We find that using several lower-order techniques for increasing the speed of numerical relativity simulations actually lead to apparently non-convergent errors. Even when using standard high-accuracy techniques, rather than seeing clean convergence, where the waveform phase is a monotonic function of grid resolution, we find that the phase tends to oscillate with resolution, possibly due to stochastic errors induced by grid refinement boundaries. Our results seem to indicate that one can obtain gravitational waveform phases to within 0.05 rad. (and possibly as small as 0.015 rad.), while the amplitude error can be reduced to 0.1%. We then compare with the waveforms obtained using the cZ4 formalism. We find that the cZ4 waveforms have larger truncation errors for a given resolution, but the Richardson extrapolation phase of the cZ4 and BSSN waveforms agree to within 0.01 rad., even during the ringdown.

💡 Research Summary

The paper presents a thorough convergence study of the moving‑puncture approach applied to an equal‑mass, non‑spinning binary‑black‑hole system, focusing on the accuracy of the resulting gravitational waveforms. Using sixth‑ and eighth‑order finite‑difference schemes, the authors examine the Hamiltonian and momentum constraints both on the horizons and in the bulk (via L₂ norms). The Hamiltonian constraint exhibits clear sixth‑ and eighth‑order convergence for both implementations, confirming that the discretisation error decreases as expected when the grid spacing is halved. The momentum constraint is more delicate: its bulk L₂ norm converges only for a consistently sixth‑order scheme, while on the horizons convergence is observed for both sixth‑ and eighth‑order runs, indicating that boundary treatment and time‑integration choices have a pronounced impact on this constraint.

The core of the investigation concerns the gravitational‑wave phase and amplitude errors across the late‑inspiral, merger, and ringdown phases. The authors find that several low‑order “speed‑up” techniques (e.g., reduced reconstruction order, larger Courant factors) introduce non‑convergent errors that mask the true convergence behavior. Even with high‑accuracy settings—high‑order differencing, small time steps, and fine adaptive‑mesh‑refinement (AMR) grids—the waveform phase does not decrease monotonically with resolution. Instead, the phase oscillates as the resolution changes, a phenomenon the authors attribute to stochastic errors generated at AMR refinement boundaries. These errors appear to inject small, random phase shifts that accumulate during the evolution, preventing a clean monotonic convergence pattern.

Quantitatively, the study demonstrates that, after careful error control, the phase error can be reduced to about 0.05 rad (and in the most favorable cases to ≈0.015 rad), while the amplitude error can be limited to roughly 0.1 %. These figures represent a significant improvement over typical NR simulations and are approaching the accuracy requirements for next‑generation gravitational‑wave data analysis.

To assess the influence of the evolution formalism, the authors also perform simulations using the conformal Z4 (cZ4) system and compare them with the standard BSSN formulation. At a given resolution, cZ4 waveforms exhibit slightly larger truncation errors than BSSN, reflecting the different constraint‑damping mechanisms. However, when Richardson extrapolation is applied to both sets of data, the extrapolated phases agree to within 0.01 rad throughout the entire signal, including the ringdown. This agreement indicates that, despite differing numerical error magnitudes, both formalisms converge to the same physical solution when properly extrapolated.

Overall, the paper highlights two key messages for the NR community. First, achieving sub‑0.1 rad phase accuracy requires not only high‑order spatial discretisation but also meticulous handling of AMR boundaries, as stochastic boundary‑induced errors can dominate the error budget. Second, the choice between BSSN and cZ4 does not materially affect the final extrapolated waveform, allowing researchers to select the formulation that best suits their stability and implementation preferences. The results provide a benchmark for future high‑precision simulations and suggest that further reductions in phase error will likely depend on improved AMR strategies or global spectral methods that eliminate refinement‑boundary artifacts.