Numerical relativity in spherical polar coordinates: Evolution calculations with the BSSN formulation

In the absence of symmetry assumptions most numerical relativity simulations adopt Cartesian coordinates. While Cartesian coordinates have some desirable properties, spherical polar coordinates appear better suited for certain applications, including gravitational collapse and supernova simulations. Development of numerical relativity codes in spherical polar coordinates has been hampered by the need to handle the coordinate singularities at the origin and on the axis, for example by careful regularization of the appropriate variables. Assuming spherical symmetry and adopting a covariant version of the Baumgarte-Shapiro-Shibata-Nakamura equations, Montero and Cordero-Carrión recently demonstrated that such a regularization is not necessary when a partially implicit Runge-Kutta method is used for the time evolution of the gravitational fields. Here we report on an implementation of the Baumgarte-Shapiro-Shibata-Nakamura equations in spherical polar coordinates without any symmetry assumptions. Using a partially implicit Runge-Kutta method we obtain stable simulations in three spatial dimensions without the need to regularize the origin or the axis. We perform and discuss a number of tests to assess the stability, accuracy and convergence of the code, namely weak gravitational waves, “hydro-without-hydro” evolutions of spherical and rotating relativistic stars in equilibrium, and single black holes.

Figure 1

A schematic representation of our cell-centered grid structure in spherical polar coordinates, for one fixed value of $\varphi$. Grid points, marked by the crosses, are placed at the center of grid cells, so that no grid point ends up at the center ($r=0$) or on the axes ($\sin \theta =0$ or $\sin \theta =\pi$). Our interior grid, bordered by solid lines in the figure, covers the region $0\le r\le r_{\max }$ and $0\le \theta \le \pi$ (as well as $0\le \varphi \le 2\pi$). As suggested by the two highlighted stencils, our fourth-order differencing scheme requires two levels of ghost zones outside of the interior grid, indicated by the dotted lines.

Figure 2

Snapshots of the metric coefficient $h_{rr}$ for an axisymmetric $m=0$ small-amplitude gravitational wave at different instances of time. For this simulation we used a grid of size $(160,40,2)$ and imposed the outer boundary at $r_{\max }=8.0$. We show data as a function of $r$ in the (arbitrary) direction $\theta =1.61$ and $\varphi =4.71$. Differences between the numerical results (marked by crosses) and the analytical solution (solid lines) are smaller than the width of the lines in this graph.

Figure 3

The norm of the error in the quantity $h_{rr}$ as a function of time for a small-amplitude, axisymmetric gravitational wave. We show results for simulations with a grid of size $(40N,10N,2)$, for $N=1$, $N=2$ and $N=4$, with the outer boundary imposed at $r=8.0$. At these early times, before the results are affected by the outer boundary, the error appears to be dominated by the fourth-order differencing of the spatial derivatives.

Figure 4

Snapshots of the metric coefficient $h_{rr}$ for a nonaxisymmetric $m=2$ small-amplitude gravitational wave at different instances of time. For this simulation we used a grid of size $(40,32,64)$ and imposed the outer boundary at $r_{\max }=4.0$; we show data as a function of $r$ in the direction $\theta =1.62$ and $\varphi =3.19$. Numerical results are marked by the crosses, while the analytical solution is shown as the solid line.

Figure 5

Snapshots of the conformal exponent $\varphi$ and the lapse $\alpha$ for a rapidly rotating star (see text for details). We show both functions both at the initial time, and at a late time $t=318M$. We also show both functions along rays in two different directions, one very close to the equator, the other pointing close to the pole. Both profiles remain very similar to their initial data throughout the evolution.

Figure 6

The difference between the maximum of the radial component of the shift vector, ${\beta }_{\max }^r$, and its initial value ${\beta }_{\max }^r{|}_{\mathrm{t}=0}$, as a function of time. For these simulations we used a grid of size $(160N,2,2)$ for $N=1$, 2, 4 and 8, and imposed the outer boundary at $r_{\max }=16.0M$. We rescale all differences with $N^2$, so that the convergence of these lines demonstrates second-order convergence in our code.

Figure 7

Initial data and final profiles of the conformal exponent $\varphi$, the lapse function $\alpha$, and the shift ${\beta }^{\widehat{r}}$, showing the coordinate transition from wormhole initial data to time-independent trumpet data. The (blue) long-dashed lines represent the initial data at $t=0$, the (red) dashed lines show our numerical results at time $t=79M$, and the (black) solid lines show the analytical trumpet solution [41]. The initial data appear double valued because we graph this function as a function of the areal radius $R$ (see text for details). For these simulations we adopted a grid size $(10240,2,2)$ with the outer boundary imposed at $r=256M$. (In these graphs we did not include the innermost two grid points, which are affected by the singular behavior of the puncture.)

Figure 8

The maximum of the radial shift ${\beta }^r$ as a function of time. We show results for different grid sizes $(1280N,2,2)$ for $N=1$, 2, 4 and 8, with the outer boundary imposed at $256M$. After a brief transition from the initial data ${\beta }^r=0$, the shift settles down into a new equilibrium. For relatively coarse grid resolutions the shift experiences a slow drift, but this drift disappears as the grid resolution is increased.

Figure 9

Profiles of violations of the Hamiltonian constraint (13) at time $t=79M$. As in Fig. 8 we show results for grid sizes $(1280N,2,2)$ for $N=1$, 2, 4 and 8, with the outer boundary imposed at $256M$. All results are rescaled with $N^2$; the convergence of the resulting lines demonstrates second-order convergence of our code.

Figure 10

Values of the norms of the components of the Ricci tensor, $||R_{ij}||$, for the flat metric (a1) with functions (a6), evaluated using grid sizes $(16N,8N,16N)$ for $N=1$, 2, 4 and 8. All values are rescaled with $N^4$, so that the convergence of these results indicates fourth-order convergence of our implementation.