High order numerical discretizations of the Einstein-Euler equations in the Generalized Harmonic formulation

TL;DR

This paper develops and validates two high-order, well-balanced numerical schemes for the Generalized Harmonic (GH) formulation of the Einstein–Euler system: a CWENO finite-difference method on Cartesian grids and an ADER-DG method on 2D unstructured polygonal meshes. Both schemes are designed to preserve exact equilibria by subtracting known stationary states and employing path-conservative formulations for nonconservative products, enabling robust long-time evolutions in vacuum spacetimes and when coupling to matter. The authors demonstrate the methods’ accuracy and stability across a comprehensive set of tests, including robust stability, linearized and gauge gravitational waves, single and rotating black holes in 2D/3D, spherical accretion, and a non-rotating neutron star in equilibrium, using up to high polynomial degrees and refined meshes. The results establish a solid foundation for extending to fully 3D simulations with moving meshes and topology changes, paving the way for complex astrophysical source modeling with DG schemes on unstructured meshes and future gravitational wave analyses.

Abstract

We propose two new alternative numerical schemes to solve the coupled Einstein-Euler equations in the Generalized Harmonic formulation. The first one is a finite difference (FD) Central Weighted Essentially Non-Oscillatory (CWENO) scheme on a traditional Cartesian mesh, while the second one is an ADER (Arbitrary high order Derivatives) discontinuous Galerkin (DG) scheme on 2D unstructured polygonal meshes. The latter, in particular, represents a preliminary step in view of a full 3D numerical relativity calculation on moving meshes. Both schemes are equipped with a well-balancing (WB) property, which allows to preserve the equilibrium of a priori known stationary solutions exactly at the discrete level. We validate our numerical approaches by successfully reproducing standard vacuum test cases, such as the robust stability, the linearized wave, and the gauge wave tests, as well as achieving long-term stable evolutions of stationary black holes, including Kerr black holes with extreme spin. Concerning the coupling with matter, modeled by the relativistic Euler equations, we perform a classical test of spherical accretion onto a Schwarzschild black hole, as well as an evolution of a perturbed non-rotating neutron star, demonstrating the capability of our schemes to operate also on the full Einstein-Euler system. Altogether, these results provide a solid foundation for addressing more complex and challenging simulations of astrophysical sources through DG schemes on unstructured 3D meshes.

Figures (13)

• Figure 1: An example of a rectangular 2D unstructured mesh covered by non-overlapping polygons $P_{i}^{n}$. Above, a zoomed-in portion highlights the polygonal structure.
• Figure 2: Time evolution of the constraints for the robust stability test case with a random initial perturbation of amplitude $\pm 10^{-7}/\varrho^2$ applied to all variables, performed on a sequence of successively refined meshes on the unit square in 2D. Top panel from left to right: fourth order ADER-DG scheme with $10$ ($\varrho = 1$), $20$ ($\varrho = 2$), and $40$ ($\varrho = 4$) elements in the $x$ direction. Bottom panel from left to right: fifth order CWENO-FD scheme on $50 \times 50$ ($\varrho = 1$), $100\times100$ ($\varrho = 2$), and $200\times200$ ($\varrho = 4$) grid points.
• Figure 3: Solution of the linearized wave test using the fourth order ADER-DG scheme on an unstructured polygonal mesh. Left panel: time evolution of the constraints. Right panel: below, 2D color contour of $\Pi_{22}$; above, comparison of the exact and numerical solution at final time $t = 1000$ on a one dimensional cut at $y = 0$.
• Figure 4: Solution of the gauge wave test with amplitude $A = 0.1$ using the fifth order CWENO-FD scheme on a Cartesian mesh. Left panel: time evolution of the constraints. Right panel: comparison of the exact and numerical solution for the variable $g_{00}$ at final time $t = 1000$ on a one dimensional cut at $y = 0$.
• Figure 5: Solution of the gauge wave test with high amplitude $A = 0.5$ using the fourth order ADER-DG scheme on an unstructured polygonal mesh. Left panel: time evolution of the constraints. Right panel: below, 2D color contour of $g_{00}$; above, comparison of the exact and numerical solution at final time $t = 1000$ on a one dimensional cut at $y = 0$.