Dailey, Schnetter, and Afshordi (Class. Quantum Grav. 42, 025002, 2025) implemented the first black-hole excision evolutions of the full nonlinear Einstein equations in the generalized harmonic formulation on a strictly rectangular Cartesian grid, using finite-difference summation-by-parts (SBP) operators with simultaneous-approximation-term (SAT) boundary penalties and an embedded-boundary treatment in which the excision surface need not coincide with grid points. Their goal is a complete, provably stable initial boundary value problem (IBVP) that can place a reflecting boundary just outside an apparent horizon to simulate black-hole echoes from putative quantum-gravity structure, and compact Cauchy domains for Cauchy-characteristic matching (CCM). The paper itself identifies what remains unsolved in its future-work section and stability discussion: (1) a mathematical framework, based on quasi-local conservation laws, that defines what reflecting gravitational radiation means in 3D and maps prescribed boundary data to changes in quasi-local mass, linear momentum, and angular momentum; (2) embedded-boundary SBP operators of higher than second order with genuine multidimensional SBP stability proofs (existing proofs are effectively one-dimensional, and the second-order D2−1D_{2-1} operator — a specific second-order-accurate SBP finite-difference stencil — was chosen because its norm stays positive definite); (3) dynamical (moving) inner boundaries with robust initialization of grid points entering the domain, needed for binary-black-hole excision; (4) non-spherical/non-convex inner and non-cubical outer boundaries; (5) reducing or eliminating the numerical dissipation currently required for stability; and (6) coupling the scheme to CCM for binary-black-hole simulations.

Background

Numerical relativity evolves Einstein's field equations as an initial boundary value problem (IBVP): given spacetime data on an initial spatial slice and conditions on the boundaries of a finite computational domain, one computes the future spacetime. In the generalized harmonic (GH) formulation, the metric gμνg_{\mu\nu} is evolved subject to gauge constraints Cμ=0\mathcal{C}_{\mu} = 0 that fix the coordinates through a prescribed wave equation, which casts the equations as a quasilinear wave system for which well-posed IBVP theory is comparatively mature. Black-hole simulations excise the interior of each hole (a region inside the horizon, causally disconnected from the exterior) so that singularities never enter the grid; the excision surface is an inner boundary on which the physical fields are purely outflowing and no boundary conditions may be over-imposed. An apparent horizon is a marginally trapped surface of the slice whose presence signals that an event horizon lies outside it. To model black-hole echoes — delayed repeated pulses that would follow a gravitational wave reflecting off structure replacing the classical horizon, as predicted by several quantum-gravity scenarios — one must place a boundary at or just outside the apparent horizon and impose boundary conditions that are genuinely reflecting for gravitational radiation, then evolve the full nonlinear system stably. Dailey, Afshordi, and Schnetter first did this in spherical symmetry (Class. Quantum Grav. 40, 195007, 2023), defining reflecting boundary conditions through conservation laws. Summation-by-parts (SBP) finite-difference operators mimic integration by parts at the discrete level: together with boundary conditions weakly imposed through simultaneous-approximation-term (SAT) penalty terms, they yield a discrete energy estimate closely tied to the well-posedness of the continuous problem, hence a route to provable numerical stability. Mattsson and Almquist (J. Comput. Phys. 334, 255, 2017) constructed embedded-boundary SBP operators, in which interpolation/extrapolation closures let an arbitrarily shaped boundary cut between grid points of a Cartesian mesh; however, their SBP proof is one-dimensional, with the multidimensional operator assembled as a tensor product that treats each grid line as its own 1D problem — which, as the source paper states plainly, falls short of an SBP proof for multidimensional problems, so numerical dissipation is added to ensure stability in practice. Sharan, Brady, and Livescu (J. Comput. Phys. 464, 111341, 2022) report a dimensionally-split Cartesian embedded-boundary method for non-dissipative schemes, and Gabbard and van Rees (J. Comput. Phys. 507, 112979, 2024) built high-order finite differences with moving immersed boundaries for continuum mechanics — adjacent machinery whose adaptation to the Einstein equations is unexplored. Separately, Cauchy-characteristic matching (CCM) joins a standard Cauchy evolution inside a worldtube to a characteristic (null-coordinate) evolution outside it, in principle extracting waveforms at future null infinity with no artificial outer boundary; Ma et al. demonstrated the first fully relativistic 3D CCM (Phys. Rev. D 109, 124027, 2024) and then nine CCM binary-black-hole merger simulations with late-time tails (Phys. Rev. D 112, 024003, 2025), but in the spectral SpECTRE code (a community spectral numerical-relativity code) with boundary-conforming spherical Cauchy domains. Hemberger et al. (Class. Quantum Grav. 30, 115001, 2013) had earlier achieved dynamical excision in spectral binary-black-hole evolutions using dual coordinate frames — a co-rotating coordinate system that carries the grid along with the holes — rather than embedded boundaries. The source paper's contribution is the marriage of all these threads in a finite-difference framework: a boundary-condition decomposition for the GH system that controls all ten incoming metric degrees of freedom (constraint-controlling, physical-curvature, and gauge-controlling operators), a proposal to use quasi-local conservation laws — balance laws relating the change of surface integrals over a closed 2-surface to fluxes through it, depending only on the geometry near the surface — to define reflection of gravitational radiation and to interpret boundary data as prescribed changes of quasi-local mass, linear momentum, and angular momentum, and the first nonlinear GR evolutions with embedded excision boundaries on a Cartesian grid. What the paper does not supply, and states as future work, defines the open problem addressed in this task.

Problem Statement

Construct a complete, provably stable embedded-boundary initial boundary value problem for the full nonlinear 3D Einstein equations in the generalized harmonic formulation on Cartesian grids, in the sense of Dailey–Schnetter–Afshordi (2025), by resolving the following jointly serving components, each of which the source states as future work. (1) Develop the proposed quasi-local-conservation-law construction into a proper mathematical framework: a definition, on a closed 2-surface in a dynamical spacetime, of what boundary data mean by "reflecting" gravitational radiation, together with an interpretation that maps a prescribed set of boundary values to the induced changes in the quasi-local mass, linear momentum, and angular momentum of the enclosed region, usable to derive boundary conditions for a desired boundary behavior in 3D. (2) Provide embedded-boundary SBP finite-difference operators of higher than second order for the GH system whose stability is established by a genuine multidimensional SBP argument rather than a line-by-line 1D proof tensorized across the grid (higher-order upwind embedded-boundary operators now exist for the scalar wave equation — Dailey 2026 — but neither the multidimensional proof nor the nonlinear GH application does), and determine whether the added numerical dissipation (coefficient ε=0.1\varepsilon=0.1 in the source) can be reduced or eliminated while retaining stability for nonlinear evolutions. (3) Extend the framework to dynamical (moving) inner boundaries tied to black-hole evolution, including a robust, accuracy-preserving way to initialize grid points that enter the computational domain as the embedded boundary moves across the grid, as required for binary-black-hole excision with finite differences. (4) Demonstrate the framework for non-spherical and non-convex inner boundaries and non-cubical outer boundaries, beyond the spherical-excision/cubical-domain tests performed to date. (5) Integrate the resulting SBP-SAT embedded-boundary GH scheme with Cauchy-characteristic matching in binary-black-hole simulations, so that compact Cauchy domains — with the outer boundary close to the strong-field region and the inner excision boundaries at or just outside the apparent horizons, imposing the quasi-local reflecting conditions of (1) — deliver gravitational waveforms extracted at infinity. The common objective is a single provably stable 3D IBVP framework in which arbitrarily shaped, moving, embedded boundaries carry physically interpretable reflecting conditions, enabling first-principles nonlinear numerical-relativity predictions of black-hole echoes and efficient CCM-based compact domains. The five components jointly serve this single objective, and the task is resolved by establishing its framework-and-stability core: component (1) — the quasi-local-conservation-law definition of reflecting gravitational radiation in 3D together with the interpretation of prescribed boundary data as induced changes of quasi-local mass, linear momentum, and angular momentum — together with, from component (2), a genuine multidimensional SBP stability argument for embedded-boundary operators of higher than second order, where a bounded sub-case of component (2) resolves that sub-case. Components (3)–(5) and the nonlinear-GH application of component (2) are extensions of the same framework, and the component-(2) determination of whether the numerical dissipation can be reduced or eliminated is independent of the operators-and-proof requirement.

References

  1. Conner Dailey, Erik Schnetter, Niayesh Afshordi, "Formulating the complete initial boundary value problem in numerical relativity to model black hole echoes", Class. Quantum Grav. 42, 025002 (2025), DOI 10.1088/1361-6382/ad9701, arXiv:2409.17970, https://arxiv.org/abs/2409.17970
  2. Conner Dailey, "Upwind embedded boundary SBP operators: New high order numerical schemes for arbitrarily shaped domains with Cartesian grids", arXiv:2607.14403 (2026), https://arxiv.org/abs/2607.14403
  3. Conner Dailey, Niayesh Afshordi, Erik Schnetter, "Reflecting boundary conditions in numerical relativity as a model for black hole echoes", Class. Quantum Grav. 40, 195007 (2023), DOI 10.1088/1361-6382/acde2f, arXiv:2301.05778, https://arxiv.org/abs/2301.05778
  4. Sizheng Ma, Mark A. Scheel, Jordan Moxon, Kyle C. Nelli, Nils Deppe, Lawrence E. Kidder, William Throwe, Nils L. Vu, "Merging black holes with Cauchy-characteristic matching: Computation of late-time tails", Phys. Rev. D 112, 024003 (2025), DOI 10.1103/jd26-8q5w, arXiv:2412.06906, https://arxiv.org/abs/2412.06906
  5. Sizheng Ma, Jordan Moxon, Mark A. Scheel, Kyle C. Nelli, Nils Deppe, Marceline S. Bonilla, Lawrence E. Kidder, Prayush Kumar, Geoffrey Lovelace, William Throwe, Nils L. Vu, "Fully relativistic three-dimensional Cauchy-characteristic matching for physical degrees of freedom", Phys. Rev. D 109, 124027 (2024), DOI 10.1103/PhysRevD.109.124027, arXiv:2308.10361, https://arxiv.org/abs/2308.10361
  6. Emanuele Berti et al., "Black hole spectroscopy: from theory to experiment", Class. Quantum Grav. 43, 123001 (2026), DOI 10.1088/1361-6382/ae59e2, arXiv:2505.23895, https://arxiv.org/abs/2505.23895
  7. Ken Mattsson, Martin Almquist, "A high-order accurate embedded boundary method for first order hyperbolic equations", J. Comput. Phys. 334, 255-279 (2017), DOI 10.1016/j.jcp.2016.12.034, https://doi.org/10.1016/j.jcp.2016.12.034
  8. Ken Mattsson, Martin Almquist, "A solution to the stability issues with block norm summation by parts operators", J. Comput. Phys. 253, 418-442 (2013), DOI 10.1016/j.jcp.2013.07.013, https://doi.org/10.1016/j.jcp.2013.07.013
  9. Nek Sharan, Peter T. Brady, Daniel Livescu, "High-order dimensionally-split Cartesian embedded boundary method for non-dissipative schemes", J. Comput. Phys. 464, 111341 (2022), DOI 10.1016/j.jcp.2022.111341, https://doi.org/10.1016/j.jcp.2022.111341
  10. James Gabbard, Wim M. van Rees, "A high-order finite difference method for moving immersed domain boundaries and material interfaces", J. Comput. Phys. 507, 112979 (2024), DOI 10.1016/j.jcp.2024.112979, https://doi.org/10.1016/j.jcp.2024.112979
  11. Daniel A. Hemberger, Mark A. Scheel, Lawrence E. Kidder, Béla Szilágyi, Geoffrey Lovelace, Nicholas W. Taylor, Saul A. Teukolsky, "Dynamical excision boundaries in spectral evolutions of binary black hole spacetimes", Class. Quantum Grav. 30, 115001 (2013), DOI 10.1088/0264-9381/30/11/115001, https://doi.org/10.1088/0264-9381/30/11/115001
  12. Marco Melis, Fabrizio Corelli, Robin Croft, Paolo Pani, "Black hole spectroscopy and nonlinear echoes in Einstein-Maxwell-scalar theory", Phys. Rev. D 111, 064072 (2025), DOI 10.1103/physrevd.111.064072, https://doi.org/10.1103/physrevd.111.064072
  13. Christina G. Taylor, Jesse Chan, "An Entropy Stable High-Order Discontinuous Galerkin Method on Cut Meshes", arXiv:2412.13002 (2024), https://arxiv.org/abs/2412.13002
  14. Ken Mattsson, David Niemelä, Andrew R. Winters, "Boundary-optimized closures for diagonal-norm upwind SBP operators", arXiv:2602.05727 (2026), https://arxiv.org/abs/2602.05727