General relativistic hydrodynamics in curvilinear coordinates

In this paper we report on what we believe is the first successful implementation of relativistic hydrodynamics, coupled to dynamical spacetimes, in spherical polar coordinates without symmetry assumptions. We employ a high-resolution shock-capturing scheme, which requires that the equations be cast in a flux-conservative form. One example of such a form is the “Valencia” formulation, which has been adopted in numerous applications, in particular in Cartesian coordinates. Here we generalize this formulation to allow for a reference-metric approach, which provides a natural framework for calculations in curvilinear coordinates. In spherical polar coordinates, for example, it allows for an analytical treatment of the singular $r$ and $\sin \theta$ terms that appear in the equations. We experiment with different versions of our generalized Valencia formulation in numerical implementations of relativistic hydrodynamics for both fixed and dynamical spacetimes. We consider a number of different tests—nonrotating and rotating relativistic stars, as well as gravitational collapse to a black hole—to demonstrate that our formulation provides a promising approach to performing fully relativistic astrophysics simulations in spherical polar coordinates.

Figure 1

Time evolution of the difference $|{\rho }_c(t)-{\rho }_c(0)|$ for the spherical relativistic star in the Cowling approximation, for both the full (lower panel) and partial approaches (upper panel) using three different resolutions in the radial direction such that the grid spacing varies as $\mathrm{\Delta }r=0.2$, 0.1, 0.05. The full approach produces noisier results initially and leads to a larger drift in the long term evolution of the central rest-mass density than the partial approach. Overall, the error decreases with increasing resolution in both approaches.

Figure 2

We show in the upper panel the time evolution of the normalized central density for the spherical relativistic star using a grid spacing of $\mathrm{\Delta }r=0.05$ for both the full (dashed line) and partial approaches (solid line). As we also saw for coarser grids, the drift in the time evolution of the central density for the full approach is larger than for the partial approach. The middle panel displays the time evolution of the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ computed inside the star for different resolutions for the partial approach. Finally, we show in the lower panel the convergence rate of the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ at $t=5$ ms is approximately 2.03 for the partial approach.

Figure 3

Upper panel: Snapshots of the rest-mass density $\rho$ at the initial time $t=0$ and at a later time $t=5$ ms for the evolution of a uniformly rotating star in the Cowling approximation. We show profiles along one ray very close to the equator, and another close to the pole. Both profiles remain very similar to their initial data throughout the evolution. Middle panel: the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ for two simulations with $N_r=100$, $N_{\theta }=8$, and $N_{\phi }=2,8$. Lower panel: the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ for three simulations performed with grids consisting of (100,8,2), (150,12,2), and (200,16,2) points, respectively.

Figure 4

We show in the upper panel the time evolution of the difference $|{\rho }_c(t)-{\rho }_c(0)|$ for the TOV model in a dynamical spacetime, using three different resolutions in the radial direction such that the grid spacing varies as $\mathrm{\Delta }r=0.2$, 0.1, 0.05. The middle panel graphs the time evolution of the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ computed inside the star, and the lower panel shows that the convergence rate of the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ at $t=15\mathrm{ms}$ is approximately 2.04.

Figure 5

Dynamical spacetime evolution of a rotating relativistic star. The figure shows the time evolution of the L1-norm, computed inside the star, of the rest-mass density $\parallel \rho (t)-\rho (0){\parallel }_1$ (upper panel) and the lapse function $\parallel \alpha (t)-\alpha (0){\parallel }_1$ (middle panel). For all simulations we impose the outer boundary at $r=60$ in our code units. Error originating from the outer boundaries reaches the center at around $t=0.3$ ms, and triggers the oscillations visible in the graph. As expected, the amplitude of the initial oscillation does not decrease with increasing resolution; however, for lower resolutions the amplitude continues to increase, while for higher resolutions it does not. We therefore measure the convergence rate before outer boundary and surface effects become the main source of error. We show (lower panel) that the convergence rate of the L1-norm $\parallel \rho (t)-\rho (0){\parallel }_1$ at $t=0.002\mathrm{ms}$ (for comparison, this corresponds to about one-tenth of the time needed by light to travel from the surface to the center, and approximately 1000 time steps) is approximately 1.98.

Figure 6

Oppenheimer-Snyder collapse of a dust cloud to a black hole. In the upper panel we show the time evolution of the central rest-mass density up to the approximate time of black hole formation. The solid (blue) line is the analytical solution for the central rest-mass density as a function of the proper time, and the (red) crosses are the numerical solution for the same quantity at a coordinate location $r=0.1M$ (which avoids numerical artifacts that are the result of the larger truncation errors caused the coordinate singularity at $r=0$). The bottom panel shows the values of the lapse at the center, ${\alpha }_c$ (blue-solid), together with its lower limit, ${\alpha }_{\text{LL}}$ (red-dashed), given by eq. (71).

Figure 7

Collapse of a marginally stable spherical star to a black hole. In the upper panel we show the time evolution of the normalized central density (measured at a coordinate radius $r=0.075$), and in the lower panel the apparent-horizon mass (solid line) in units of the ADM mass of the system (dashed line).

Figure 8

Radial profile of the conformal factor $\psi$ at time $t=300$ for the collapse of a marginally stable star to a black hole. The (red) crosses mark our numerical results, while the (blue) line is the analytical solution for a maximally sliced trumpet solution (see [57]).