A well-balanced discontinuous Galerkin method for the first--order Z4 formulation of the Einstein--Euler system

TL;DR

This work targets the stable, high-order numerical solution of the coupled Einstein–Euler equations by employing a first-order hyperbolic Z4 formulation with strong hyperbolicity on general metrics. The authors introduce a novel well-balanced ADER-DG scheme with an a posteriori sub-cell finite volume limiter, augmented by a discrete equilibrium subtraction to preserve equilibria exactly and a robust primitive-variable recovery filter for low-density regimes. They demonstrate the method’s capabilities through a comprehensive suite of tests, including long-time evolutions of stationary black holes (even with extreme spin), a neutron star in vacuum, and the head-on collision of puncture black holes, achieving unprecedented accuracy and robustness in 3D GR settings. The results highlight significant advances in accurate quasi-normal mode studies, neutron star dynamics without artificial atmospheres, and access to previously intractable binary black-hole configurations within the Z4 framework, with potential broad impact on numerical relativity and relativistic astrophysics.

Abstract

In this paper we develop a new well-balanced discontinuous Galerkin (DG) finite element scheme with subcell finite volume (FV) limiter for the numerical solution of the Einstein--Euler equations of general relativity based on a first order hyperbolic reformulation of the Z4 formalism. The first order Z4 system, which is composed of 59 equations, is analyzed and proven to be strongly hyperbolic for a general metric. The well-balancing is achieved for arbitrary but a priori known equilibria by subtracting a discrete version of the equilibrium solution from the discretized time-dependent PDE system. Special care has also been taken in the design of the numerical viscosity so that the well-balancing property is achieved. As for the treatment of low density matter, e.g. when simulating massive compact objects like neutron stars surrounded by vacuum, we have introduced a new filter in the conversion from the conserved to the primitive variables, preventing superluminal velocities when the density drops below a certain threshold, and being potentially also very useful for the numerical investigation of highly rarefied relativistic astrophysical flows. Thanks to these improvements, all standard tests of numerical relativity are successfully reproduced, reaching three achievements: (i) we are able to obtain stable long term simulations of stationary black holes, including Kerr black holes with extreme spin, which after an initial perturbation return perfectly back to the equilibrium solution up to machine precision; (ii) a (standard) TOV star under perturbation is evolved in pure vacuum ($ρ$=$p$=0) up to t=1000 with no need to introduce any artificial atmosphere around the star; and, (iii) we solve the head on collision of two punctures black holes, that was previously considered un--tractable within the Z4 formalism.

Figures (15)

• Figure 1: Lelft panel: plot of the polynomial $f(y)$ built in such a way to have $f(0)=1$, $f(y_0)=0$, $f^\prime(0)=0$, $f^\prime(y_0)=0$. Right panel: plot of the ratio between velocity and momentum $v_i/S_i$ when the filter function is applied. If $y=\rho h W^2 < 10^{-4}$, the filter is activated, and the velocity decreases smoothly to zero.
• Figure 2: Linearized gravitational wave test solved with an ADER-DG scheme of order 6. Left panel: $K_{zz}$ component of the extrinsic curvature at the final time, compared to the exact solution. Right panel: Einstein constraints monitored all along the duration of the simulation.
• Figure 3: Solution of the gauge wave test at $t=1000$ with $A=0.1$ using an ADER DG scheme of order 4. Left panel: profile of the lapse $\alpha$ compared to the exact solution. Right panel: Evolution of the Einstein constraints.
• Figure 4: Robust stability test case with a random initial perturbation of amplitude $10^{-7}/\rho^2$ in all quantities on a sequence of successively refined meshes on the unit square in 2D. The gamma--driver shift condition, $1+\log$ slicing and ADER-DG $P_3$ scheme have been used. Top left: $10\times10$ elements, corresponding to $40\times40$ degrees of freedom ($\varrho=1$). Top right: $20\times20$ elements, corresponding to $80\times80$ degrees of freedom ($\varrho=2$). Bottom left: $40\times40$ elements, corresponding to $160\times160$ degrees of freedom ($\varrho=4$). Bottom right: $80\times80$ elements, corresponding to $320\times320$ degrees of freedom ($\varrho=8$).
• Figure 5: Solution of the spherical accretion of matter onto a non-rotating black hole. The final rest mass density (left panel) and the radial velocity (right panel) at time $t=1000\,M$ are compared to their initial profiles.