Accurate black hole evolutions by fourth-order numerical relativity

We present techniques for successfully performing numerical relativity simulations of binary black holes with fourth-order accuracy. Our simulations are based on a new coding framework which currently supports higher-order finite differencing for the Baumgarte-Shapiro-Shibata-Nakamura formulation of Einstein’s equations, but which is designed to be readily applicable to a broad class of formulations. We apply our techniques to a standard set of numerical relativity test problems, demonstrating the fourth-order accuracy of the solutions. Finally we apply our approach to binary black hole head-on collisions, calculating the waveforms of gravitational radiation generated and demonstrating significant improvements in waveform accuracy over second-order methods with typically achievable numerical resolution.

Figure 1

The $L_2$ norm of $\delta {\gamma }_{xx}={\gamma }_{xx}-{\gamma }_{xx}^{\mathrm{analytic}}$, rescaled by ${\rho }^4/16$, for the one-dimensional gauge wave test with $A=0.01$. Note the near perfect overlap for 630 crossing times and that the larger $\rho$ runs are convergent longer. The runs are acceptably accurate when the (nonrescaled) norm of the error is smaller than ${10}^{-4}$ (i.e. 1% of $A$).

Figure 2

The $L_2$ norm of $\mathcal{H}$, rescaled by ${\rho }^4/16$, for the one-dimensional gauge wave test with $A=0.01$. Note the near perfect overlap for 800 crossing times and that the larger $\rho$ runs are convergent longer. The runs are acceptably accurate when the (nonrescaled) norm of the Hamiltonian constraint is smaller than ${10}^{-4}$ (i.e. 1% of $A$).

Figure 3

The $L_2$ norm of $\mathcal{H}$, rescaled by ${\rho }^4/16$, for the one-dimensional gauge wave test with $A=0.1$. Note the near perfect overlap for 80 crossing times and that the larger $\rho$ runs are convergent longer. The runs are acceptably accurate when the (nonrescaled) norm of the Hamiltonian constraint is smaller than ${10}^{-3}$ (i.e. 1% of $A$).

Figure 4

The $L_2$ norm of $\mathcal{H}$, rescaled by ${\rho }^4/16$, for the two-dimensional gauge wave test with $A=0.1$. Note the near perfect overlap for 55 crossing times and that the larger $\rho$ runs are convergent longer. The runs are acceptably accurate when the (nonrescaled) norm of the Hamiltonian constraint is smaller than ${10}^{-3}$ (i.e. 1% of $A$).

Figure 5

The $L_2$ norm of $\mathcal{H}$, rescaled by ${\rho }^4/16$, for the one-dimensional Gowdy wave test. Note the good agreement between the curves for 1000 crossing times and that the evolution is backwards in time.

Figure 6

The $L_2$ norm of the error in ${\gamma }_{zz}$, rescaled by the $L_2$ norm of ${\gamma }_{zz}$ and by ${\rho }^4/16$, for the one-dimensional Gowdy wave test. Note the good agreement between the curves for 1000 crossing times and that the relative error in ${\gamma }_{zz}$ is less than $2·{10}^{-4}$ at 1000 crossing times. The evolution is backwards in time.

Figure 7

The $L_2$ norm of $\mathcal{H}$, rescaled by ${\rho }^4/16$, for the two-dimensional Gowdy wave test. The test was limited to 50 crossing times due to limited resources. The early time lack of convergence is due to low amplitude high frequency noise in the $\rho =8$ results. Although the two higher resolution curves do not lie on top of $\rho =2$ curve, the overlap of the two high resolution curves indicates that the Hamiltonian constraint is converging to fourth order. Note that the evolution is backwards in time.

Figure 8

The $L_2$ norm of the error in ${\gamma }_{yy}$, rescaled by ${\rho }^4/16$, for the two-dimensional Gowdy wave test. The test was limited to 50 crossing times due to limited resources. The excellent overlap of all three curves indicates that the code is converging to fourth order. Note that the evolution is backwards in time.

Figure 9

The top plot shows the $(\ell =2,m=0)$ mode of ${\psi }_4$ at $r=5M$ for the second-order evolution of the MWBL data for grid sizes $h=1.1515/(9\rho )$, as well as the Richardson extrapolated value. The bottom plot shows the differences $({\psi }_4{|}_{\rho =1}-{\psi }_4{|}_{\rho =2})$ and $4({\psi }_4{|}_{\rho =2}-{\psi }_4{|}_{\rho =4})$ between these waveforms for this mode. Note that the latter difference has been rescaled by a factor of 4. Second-order convergence is demonstrated by the reasonable overlap between these two differences.

Figure 10

Convergence rate of the $(\ell =2,m=0)$ mode of ${\psi }_4$ at $r=5M$ for a purely second-order evolution of MWBL head-on collision data.

Figure 11

The $L_2$ norm of the Hamiltonian constraint violation versus time for the second-order accurate head-on collision runs. Note that the constraint violation increases with resolution.

Figure 12

Hamiltonian constraint violation (absolute value) versus $z$ along the $z$ axis at time $t=80M$ for the $\rho =4$ second-order head-on collision run. Note that the large constraint violation inside the apparent horizon ($r=2.8M$) has not leaked out into the exterior spacetime.

Figure 13

The top plot shows the $(\ell =2,m=0)$ mode of ${\psi }_4$ at $r=5M$, produced using fourth-order evolution with LOR, for 3 resolutions, along with the Richardson extrapolation curve of these data. The $\rho =2$ and $\rho =4$ curves are indistinguishable from the extrapolated curve. The bottom plot shows rescaled differences ${\psi }_4{|}_{\rho =1}-{\psi }_4{|}_{\rho =2}$ and $16({\psi }_4{|}_{\rho =2}-{\psi }_4{|}_{\rho =4})$ (note the factor of 16). The convergence rate for the fourth-order runs using LOR inside the apparent horizons is not strictly fourth order due to the phase discrepancy between the two curves, but average waveform differences between resolutions is consistent with fourth-order accuracy. Note that the amplitude of the dotted curve is less than the amplitude of the solid curve indicating better than fourth-order reduction in the amplitude of the errors.

Figure 14

A comparison of the $(\ell =2,m=0)$ mode of ${\psi }_4$ at $r=5M$ for second and fourth-order evolutions with LOR. Note that the $\rho =1$ fourth-order waveform is better than the $\rho =2$ second-order waveform. Similarly the $\rho =2$ fourth-order waveform is better than the $\rho =4$ second-order waveforms. The Richardson extrapolated waveform (based on the second-order waveforms) is indistinguishable from the $\rho =2$ and $\rho =4$ fourth-order waveforms.

Figure 15

A comparison of the $(\ell =2,m=0)$ mode of ${\psi }_4$ at $r=5M$ for second- and fourth-order evolutions with LOR. Note that the $\rho =2$ fourth-order waveform is better than the $\rho =4$ second-order waveforms. The second-order $\rho =4$ curve lies above the extrapolated curve while the fourth-order $\rho =1$ curve lies below.

Figure 16

The $(\ell =4,m=0)$ mode of ${\psi }_4$ for $\rho =4$ as well as the differences between the $\rho =1$ and $\rho =2$ waveforms and the $\rho =2$ and $\rho =4$ waveforms (the latter rescaled by 16). The convergence rate of the $(\ell =4,m=0)$ mode of ${\psi }_4$ at $r=5M$ for the fourth-order runs with LOR is better than fourth order. Note that the error in the $\rho =1$ waveform (as evident by the difference between the $\rho =1$ and $\rho =2$ waveforms) is approximately 50% of the amplitude of the waveform at $20M$ and larger than the waveform beyond $30M$.

Figure 17

The $(\ell =4,m=0)$ mode of ${\psi }_4$ at $r=5M$ for fourth-order evolutions with LOR as well as second-order evolutions. Note that the $\rho =1$ fourth-order waveform has an order 100% error and that the $\rho =2$ fourth-order waveform has a substantially smaller error than the $\rho =4$ second-order waveform.

Figure 18

The absolute value of the Hamiltonian constraint along the $z$ axis at $t=47.2M$ for the second- and fourth-order runs. The points near the puncture have been removed (the violation at the puncture was ${10}^{117}$) from the fourth-order data. Note that the constraint violation is ${10}^4$ times bigger inside the horizon for the fourth-order run, and that the extreme constraint violation does not leak out of the horizon (located at $z=2.8M$).

Figure 19

The unstable mode in $\mathcal{H}$ for the $\rho =4$ fourth-order runs.

Figure 20

The $L_2$ norm, restricted to the region $3M<r_{\mathrm{physical}}<26M$, of all constraint violations. The boundaries were located at $63M$. All norms have been multiplied by ${10}^5$ and the $\rho =2$ norms have been rescaled by an additional factor of 16. In each panel the solid curves correspond to $\rho =1$ and the dashed curves to $\rho =2$. The reasonable overlap of the dashed and solid curves indicate that the constraint violations are converging to zero to fourth order. The constraints no longer converge to zero after $t\sim 43M$ due to boundary effects.

Figure 21

The $L_2$ norm of the Hamiltonian constraint violation for the $\rho =4$ fourth-order evolution with LOR for the standard and modified $1+\log $ lapses. Note the $e^{t^2}$ blowup in the run using the standard $1+\log $ lapse and that the run using the modified $1+\log $ lapse crashes at $26.5M$.

Figure 22

The horizon mass versus time for the fourth-order runs with LOR. The thick solid line and the dotted line show the horizon mass versus time for the $\rho =2$ resolution with the physical boundary at $63M$ and $26M$, respectively. Note that the sharp increase in mass at later times is due to boundary effects. The dashed line shows the horizon mass for the $\rho =4$ run with the boundaries at $26M$. The error in the horizon mass at $t=63.5M$ is 2.6% for the $\rho =2$ run with boundary at $26M$.

Figure 23

The Hamiltonian constraint along the $z$ axis at $t=40.8M$ for the fourth-order (with LOR) runs with grid spacing $h=1.1515/18$ (i.e. $\rho =2$). $z_{\mathrm{fish}\mathrm{eye}}=12.3M$ corresponds to $26M$ in physical coordinates and $z_{\mathrm{fish}\mathrm{eye}}=24.6M$ corresponds to $63M$ in physical coordinates. Boundary contamination for the larger run has just reached $z_{\mathrm{fish}\mathrm{eye}}=12M$. Note that in the range $5M<z_{\mathrm{fish}\mathrm{eye}}<12M$ the Hamiltonian constraint is 100 times smaller when the boundary is pushed out to $63M$.

Figure 24

The $(\ell =2,m=0)$ mode of ${\psi }_4$ (observer at $r=5M$) for the fourth-order (with LOR) $\rho =2$ run with the physical boundary at $26M$ (standard) and $63M$. The effect of the boundary is not significant until $t=55M$.

Figure 25

The $(\ell =4,m=0)$ mode of ${\psi }_4$ (observer at $r=5M$) for the fourth-order (with LOR) $\rho =2$ run with the physical boundary at $26M$ (standard) and $63M$. The effect of the boundary is significant at $t>33M$.

Figure 26

Convergence plot for the mixed fourth-ordered centered differencing with second-order upwinded differencing runs. Note that the difference ${\psi }_4{|}_{\rho =2}-{\psi }_4{|}_{\rho =4}$ has been multiplied by 4, indicating approximate second-order convergence.

Figure 27

A comparison of the $(\ell =2,m=0)$ mode of ${\psi }_4$ at $r=5M$ produced using the standard second-order evolution and a mixed fourth-order with second-order upwinding evolution. Note that the $\rho =2$ waveform from the mixed fourth-order/second-order upwinding is of similar quality to the $\rho =4$ second-order waveform. The extrapolated values used in this plot are based on the waveforms from the fourth-order runs with LOR.