High-resolution shock capturing scheme for ideal hydrodynamics in general relativity optimized for quasistationary solutions

In this article, we present a numerical scheme for relativistic ideal hydrodynamics on arbitrary curved spacetimes in up to three spatial dimensions. The aim of this scheme is an improved treatment of stationary or quasistationary scenarios without loss of general applicability. This is achieved by treating pressure and gravitational forces, which are usually computed with unrelated numerical methods, within a unified framework. We derive a special formulation of the source terms in the hydrodynamic evolution equations, which is then utilized to construct the new scheme. For single star simulations, the optimizations take full effect only for rigidly rotating isentropic stars in a corotating coordinate frame. To test the algorithm, simulations of single polytropic neutron stars (rigidly rotating and nonrotating) have been carried out, demonstrating its usefulness for studies of stellar oscillation modes. We found that most of the numerical error is caused by the treatment of the stellar surface. To study the evolution scheme itself, we introduce a test bed involving an artificial spacetime without a fluid-vacuum boundary. The oscillation modes of the stellar models used here have already been investigated in the literature. A comparison with our results is given, exhibiting good agreement.

Figure 1

Illustration of the modified reconstruction scheme on the example of a polytropic gas at rest in a homogeneous Newtonian gravitational field. Plotted is the specific energy density. (a) Standard piecewise constant/linear reconstruction. (b) $\mathrm{\text{Piecewise equilibrium}}+\mathrm{\text{linear reconstruction}}$ if the true profile is close to equilibrium. (c) Same case if the system is in equilibrium. (d) Same case if the profile is homogeneous, thus far from equilibrium. The fluxes themselves are not captured very well, but the flux differences are still exact (zero) in this particular example. However, the source terms, which are proportional to the height of the sawtooth, are dominant in this case anyway.

Figure 2

Evolution of the radial velocity at $r=7\mathrm{km}$, $\vartheta =\pi /4$ for a polytropic, axisymmetric 2D-simulation of model N with a resolution of ${10}^2$ cells. The corresponding plots for 1D and 3D simulations are qualitatively the same.

Figure 3

Results of a 3D simulation of model N, with a resolution of ${47}^3$ points. The initial data was perturbed with $l=2$, $A=0.001$. Shown are time evolution and corresponding Fourier spectrum of the velocity in $\vartheta$-direction on the space diagonal at $R=7\mathrm{km}$. The vertical lines mark mode frequencies for the same stellar model as reported by 9.

Figure 4

Results obtained with a resolution of only ${10}^3$ points in 3D, for model N and a polytropic EOS. The initial data was perturbed with $l=0$, $A=0.001$. (a) Time evolution of the radial velocity $v_r$ at $R=7\mathrm{km}$ on the space diagonal. (b) Radial velocity divided by an exponential decay ($\tau =1.07292\mathrm{ms}$) fitted to its maxima. Note the higher modes are damped away almost instantaneously.

Figure 5

Radial eigenfunctions (plotted in arbitrary units) for the $H_1$-mode, obtained via mode recycling from a 1D simulation. $r$ is the circumferential radius. Evolution time was 20 ms. Different resolutions are plotted for comparison.

Figure 6

(a) Convergence of the frequency for the lowest (known) mode of the 1D and 2D toy models. Shown is the relative difference to the “true” frequency $f_0$ obtained by fitting convergence law Eq. (80) to the data. The fit results are $f_0=2.181\mathrm{kHz}$, $p=1.74$ in 1D and $f_0=1.257\mathrm{kHz}$, $p=1.93$ in 2D. (b) Convergence of the numerical damping time scale $\tau$ (in units of the oscillation period $T_o$) for the same mode. The convergence orders obtained from the fits (straight lines) to the convergence law $\tau =b{N}^p$ are $p=3.24$ (1D) and $p=3.14$ (2D). Note the lowest resolution ($N=15$) was excluded for all fits.

Figure 7

Convergence behavior for the $F$- and $H_1$-modes of model N. (a) Convergence of the frequencies $f$. $f_0$ is the frequency obtained from a 1D simulation with a resolution of 1133. In detail $f_0=2.6845679\mathrm{kHz}$ for the $F$-mode and $f_0=4.5456764\mathrm{kHz}$ for the $H_1$-mode. (b) Convergence of the numerical damping time scale $\tau$ in units of the respective oscillation period $T_o$. The straight lines are fits of the convergence law $\tau =a\Delta x^p$ to the data. The convergence order is $p=1.62$ ($H_1$-mode, 2D), $p=1.80$ ($F$-mode, 2D), $p=1.61$ ($F$-mode, 3D), and $p=2.34$ ($H_1$-mode, 1D). Note: The damping time for the highest resolution 2D $F$-mode run could not be determined due to insufficient evolution time (10 ms).

Figure 8

Comparison of the ${}^{2}f$- and ${}^{2}p_1$-mode eigenfunctions of model N in Cowling approximation, as obtained with our $\mathtt{P}\mathtt{i}\mathtt{z}\mathtt{z}\mathtt{a}$ code and the $\mathtt{C}\mathtt{o}\mathtt{C}\mathtt{o}\mathtt{N}\mathtt{u}\mathtt{T}$ code described in 4. The eigenfunctions are plotted in arbitrary units along the equatorial plane. The velocities are the physical radial velocities $v^r\sqrt{g_{rr}}$, and $r$ is the circumferential radius. The resolution is ${189}^2$ points for the $\mathtt{P}\mathtt{i}\mathtt{z}\mathtt{z}\mathtt{a}$ results, and $140\times 60$ points in $(r,\vartheta )$ for the $\mathtt{C}\mathtt{o}\mathtt{C}\mathtt{o}\mathtt{N}\mathtt{u}\mathtt{T}$ results.

Figure 9

Like Fig. 8, but for the radial $F$- and $H_1$-mode eigenfunctions. The $\mathtt{P}\mathtt{i}\mathtt{z}\mathtt{z}\mathtt{a}$ simulations were performed with a resolution of 1133 points in 1D. The $\mathtt{C}\mathtt{o}\mathtt{C}\mathtt{o}\mathtt{N}\mathtt{u}\mathtt{T}$ results were obtained in 2D with a resolution of $140\times 60$ points in $r$ and $\vartheta$.