3D simulations of Einstein's equations: Symmetric hyperbolicity, live gauges, and dynamic control of the constraints

We present three-dimensional simulations of Einstein's equations implementing a symmetric hyperbolic system of equations with dynamical lapse. The numerical implementation makes use of techniques that guarantee linear numerical stability for the associated initial-boundary value problem. The code is first tested with a gauge wave solution, where rather larger amplitudes and for significantly longer times are obtained with respect to other state of the art implementations. Additionally, by minimizing a suitably defined energy for the constraints in terms of free constraint-functions in the formulation one can dynamically single out preferred values of these functions for the problem at hand. We apply the technique to fully three-dimensional simulations of a stationary black hole spacetime with excision of the singularity, considerably extending the lifetime of the simulations.

Figure 1

This figure shows the logarithm of the energy for the constraints, in simulations of the gauge wave with amplitude $A=0.5$. Three different resolutions are used, the coarsest run exhibits a slow growth which is negligible at higher resolutions.

Figure 2

The logarithm of the relative error in the maximum value attained by $g_{xx}$ (compared to its exact value) versus number of crossing times for the simulations of Fig. 1.

Figure 3

Convergence of the constraints energy to zero as a function of the number of crossing times, for the simulations of Fig. 1, 2. This figure shows the convergence factors obtained by dividing the energy found at different resolutions in consecutive pairs. The oscillations observed are a product of phase velocities differing at different resolutions.

Figure 4

Logarithm of the energy for the constraints for the $A=1$ gauge wave simulations. The coarsest resolution exhibits a marked growth after 600 crossing times, though with higher resolutions this effect becomes negligible.

Figure 5

The logarithm of the relative error in the maximum attained by $g_{xx}$ for the simulations of Fig. 4. The errors exhibit a slow growth which diminishes with resolution.

Figure 6

Convergence of the constraint energy to zero, for the simulations of Figs. 4, 5, obtained by dividing the energy for each resolution in consecutive pairs. As in Fig. 3, after some time the convergence factors deviate from the expected value of two, but this difference diminishes when the two highest resolutions are used to compute the convergence factor.

Figure 7

Snapshots at 280 and 560 crossing times displaying the exact solution and the numerically calculated ones for the middle ($△=0.0125$) and fine ($△=0.00625$) resolutions. Convergence to the exact solution is evident as resolution is increased.

Figure 8

Convergent factor calculated with the numerical solutions for the middle and fine resolutions. At late times accumulation of errors makes the obtained value deviate from its expected one.

Figure 9

$L_2$ norm of error in the numerical solution. Although a faster than linear growth is observed for the middle resolution, this effect converges away as a finer resolution is considered.

Figure 10

Black hole simulations, with fixed outer boundaries (at 5  M), and different position for the inner boundary (IB). The resolution is $\Delta =\mathrm{M}/5$.

Figure 11

Two-constraint-function family, with fixed values $\gamma =0=\eta$, inner and outer boundaries at 1.1  M and 5  M, respectively.

Figure 12

Generally $\chi$ would need to change sign in order to achieve certain control on the constraints. Therefore, it would not have a continuous limit when $△t\to 0$. For this reason, the two-constraint-function formulation is used throughout this paper.

Figure 13

Energy for the constraints and its time derivative, $\dot{\mathcal{N}}$, computed through the semidiscrete prediction and through numerical differentiation. The remarkable agreement between the two indicates that the semidiscrete analysis used to calculate the constraint-minimization is faithful representation of the fully discrete evolution.

Figure 14

Time derivative of the energy of the constraints computed through the semidiscrete prediction and through numerical differentiation for the gauge wave case. The agreement of the two curves is evident.

Figure 15

The plots in this figure and those in Fig. 16 show the effects of varying the frequency of performing the constraint-minimization, as well as the dependence on $n_a$. The constraint-function $\eta$ for $n_a={10}^4$ is not shown as the run crashes very early and the scale in the time axis for the associated plot is not logarithmic.

Figure 16

This plot shows simulations as those of Fig. 15, except that here the minimization is done every ten iterations. Comparing it with the previous plot, it seems clear that performing the minimization at every iteration seems a better option.

Figure 17

This figure and Fig. 18 demonstrate the influence of the tolerance value $T$ on the results of the constraint-minimization. The two-constraint-function formulation is used, with $\gamma =0$ and $\eta (t)$ chosen dynamically. The results are not sensitive to the chosen tolerance value, provided that one avoids large variations in $\eta (t)$ by using appropriate values of $n_a$. However, notice that “appropriate” values of $n_a$ naturally keep the value of the energy close to its initial, discretization value (here given by $0.99925\times {10}^{-4}$. Therefore, keeping the constraint energy near its initial discretization value appears to give the best results. This figure shows runs with $T={10}^{-2}$ and $T={10}^{-3}$. Figure 18 shows runs with $T={10}^{-4}$ and $T={10}^{-5}$.

Figure 18

This figure compares results for $T={10}^{-4}$ and $T={10}^{-5}$. See Fig. 17 for additional information.

Figure 19

These plots show the associated $\eta (t)$ for the $n_a=1,1000$ cases of the previous plot, in more detail. Notice the scale in the $n_a=1$ case; $\eta$ changes a lot per time step. On the other hand, notice how in the $n_a={10}^4$ case $\eta (t)$ changes slower but still changes, and eventually it needs to take very large values in order to control the constraints.

Figure 20

Dynamic minimization done with boundaries at 5  M, $△=\mathrm{M}/5$, $T={10}^{-3}$, and $n_a={10}^3$.

Figure 21

Same as previous figure, but keeping the constraint-functions constant after 750  M. The figure compares the resulting energy for the constraints with that of the previous figure (shown at late times only, since because of the setup the runs are identical up to $t=750\mathrm{M}$).

Figure 22

Convergence test, with dynamic constraint-functions kept constant after 750  M.

Figure 23

Apparent horizon mass, with dynamic constraint-function values. The run with higher resolution ran out of computing time (in the equivalent plots of the previous figure an apparent horizon was not searched for).

Figure 24

Comparison between dynamic and fixed, fine-tuned constraint-function values given by the asymptotic values to which they approach [cf. Equation (36)].

Figure 25

Dynamic minimization done with boundaries at 10  M, $△=\mathrm{M}/5$, $T={10}^{-4}$, and $n_a={10}^2$.

Figure 26

Running with boundaries at 10  M, using fixed, but fine-tuned constraint-functions, obtained from runs with boundaries at 5  M and 10  M.

Figure 27

Apparent horizon mass for the simulation of Fig. 25 and the simulation of Fig. 26 with constraint-functions given by Eq. (37). A logarithmic scale in time is used in order to show the oscillations in more detail. The oscillations are not caused by time variations in the constraint-functions, since they are also present in the fixed-constraint-functions case.

Figure 28

Dynamic minimization done with boundaries at 15  M, $△=\mathrm{M}/5$, $T={10}^{-5}$, and $n_a={10}^2$. The constraint-functions are constant for $t\ge 450\mathrm{M}$, where they are $\eta =-0.135$, $\gamma =-3.389$; thus, there is no response of the constraint-functions when the code is about to crash.

Figure 29

This figure shows a fully 3D run with fixed and fined-tuned parameters, given by Eq. (37), compared to the same run in octant symmetry, which was stopped at 18 000  M.