Combining spectral and shock-capturing methods: A new numerical approach for 3D relativistic core collapse simulations

We present a new three-dimensional general relativistic hydrodynamics code which is intended for simulations of stellar core collapse to a neutron star, as well as pulsations and instabilities of rotating relativistic stars. Contrary to the common approach followed in most existing three-dimensional numerical relativity codes which are based in Cartesian coordinates, in this code both the metric and the hydrodynamics equations are formulated and solved numerically using spherical polar coordinates. A distinctive feature of this new code is the combination of two types of accurate numerical schemes specifically designed to solve each system of equations. More precisely, the code uses spectral methods for solving the gravitational field equations, which are formulated under the assumption of the conformal flatness condition (CFC) for the three-metric. Correspondingly, the hydrodynamics equations are solved by a class of finite difference methods called high-resolution shock-capturing schemes, based upon state-of-the-art Riemann solvers and third-order cell-reconstruction procedures. We demonstrate that the combination of a finite difference grid and a spectral grid, on which the hydrodynamics and metric equations are, respectively, solved, can be successfully accomplished. This approach, which we call Mariage des Maillages (French for grid wedding), results in high accuracy of the metric solver and, in practice, allows for fully three-dimensional applications using computationally affordable resources, along with ensuring long-term numerical stability of the evolution. We compare our new approach to two other, finite difference based, methods to solve the metric equations which we already employed in earlier axisymmetric simulations of core collapse. A variety of tests in two and three dimensions is presented, involving highly perturbed neutron star spacetimes and (axisymmetric) stellar core collapse, which demonstrate the ability of the code to handle spacetimes with and without symmetries in strong gravity. These tests are also employed to assess the gravitational waveform extraction capabilities of the code, which is based on the Newtonian quadrupole formula. The code presented here is not limited to approximations of the Einstein equations such as CFC, but it is also well suited, in principle, to recent constrained formulations of the metric equations where elliptic equations have a preeminence over hyperbolic equations.

Figure 1

Radial setup of the initial spectral grid (collocation points are marked by plus symbols) and the time-independent finite difference grid (cell centers are marked by filled circles, separated by cell interfaces symbolized by vertical dashes) for a typical core collapse simulation. The upper panel shows the innermost 500 km containing the nucleus (ending at $r_{\mathrm{d}}\approx 200\mathrm{km}$), the first shell, and a part of the second shell of the spectral grid. In the lower panel a part of the last regular shell (which is confined by the outer boundary of the finite difference grid at $r_{\mathrm{fd}}\approx 2200\mathrm{km}$) and the beginning of the compactified domain of the spectral grid are plotted. The domain boundaries are indicated by vertical dotted lines.

Figure 2

Comparison of the convergence behavior for the three-metric solvers in 2D. For solver 1 (filled circles), the maximum increment $\Delta {\widehat{u}}_{\max }^s$ per iteration $s$ decreases to the threshold $\Delta {\widehat{u}}_{\mathrm{thr}}^s={10}^{-15}$ (lower horizontal dotted line) within less than 10 iterations, while solver 3 (asterisks) needs more than 40 iterations to reach its (less restrictive) threshold (upper horizontal dotted line) of ${10}^{-6}$. The very low relaxation factor needed for solver 2 (filled squares) results in a remarkably slow convergence, requiring more than 700 iterations. The solid lines mark the approximate linear decrease of $\log \Delta {\widehat{u}}_{\max }^s$.

Figure 3

Equatorial profile of the shift vector component ${\beta }_{\phi \mathrm{e}}$ obtained by different metric solvers compared with the correct profile from the initial data solver (solid line) for the rotating neutron star model RNS. Because of its approximate boundary value, the profile from solver 1 (dashed line) shows large deviation from the correct solution, particularly for a grid boundary $r_{\mathrm{fd}}$ close to the stellar equatorial radius $r_{\mathrm{s}\mathrm{e}}$ (upper panel). As solver 2 (dashed-dotted line) needs no explicit boundary conditions, its solution matches well with the correct solution, with improving agreement as $r_{\mathrm{fd}}$ is at larger distance from $r_{\mathrm{s}\mathrm{e}}$ (lower panel). The compactified radial grid of solver 3 (dotted line) fully accounts for noncompact support terms, and thus agrees very well with the correct solution, independent of the location of $r_{\mathrm{fd}}$. The radii $r_{\mathrm{s}\mathrm{e}}$ and $r_{\mathrm{fd}}$ are indicated by vertical dotted lines.

Figure 4

Time evolution of the central conformal factor ${\varphi }_{\mathrm{c}}$ (upper panel) for the core collapse model SCC, using metric solver 1 (solid line), 2 (dashed line), and 3 (dashed-dotted line), respectively. All three solvers yield similar results. The small relative differences of less than ${10}^{-3}$ in ${\varphi }_{\mathrm{c}}$ (lower panel) obtained with solvers 1 and 3 (solid line) and solvers 2 and 3 (dashed line) prove that numerical variations of the metric from each solver are of the order of the small overall discretization error of the entire evolution code. The time of bounce $t_{\mathrm{b}}$ is indicated by the vertical dotted line.

Figure 5

Two different methods for determining the domain radii of the spectral grid boundary. The upper panel shows the time evolution of the domain radius parameter $r_{\mathrm{d}}$ for the core collapse model SCC, where $r_{\mathrm{d}}$ is either set by the sonic point method (solid line; sonic point first detected at $t\sim 23\mathrm{ms}$) or by the rest-mass fraction method (dashed line). The boundary of the finite difference grid $r_{\mathrm{fd}}$, the stellar equatorial radius $r_{\mathrm{s}\mathrm{e}}$, the minimal domain radius $r_{\mathrm{d}\min }$ (set to $r_{40}$), and the approximate location of shock formation $r_{\mathrm{sh}}$ are indicated by horizontal dotted lines. The relative difference between the values of ${\varphi }_{\mathrm{c}}$ from simulations using the two methods (lower panel) is less than ${10}^{-4}$ throughout the evolution. The time of bounce $t_{\mathrm{b}}$ is indicated by the vertical dotted line.

Figure 6

Importance of the correct spectral domain setup for highly dynamic simulations, shown for the core collapse model SCC. If the domain radius parameter is not reasonably adjusted (upper panel), e.g. $r_{\mathrm{d}}$ is held fixed at 10% of the initial stellar equatorial radius (dashed line), or if the minimal domain radius is too large, $r_{\mathrm{d}\min }=r_{100}$ (dashed-dotted line), the central conformal factor ${\varphi }_{\mathrm{c}}$ deviates strongly from the correct value (solid line; cf. Fig. 4). If the number of domains is too small (lower panel), e.g. $n_{\mathrm{d}}=3$ (dashed line) instead of $n_{\mathrm{d}}=6$ (solid line), the metric inside the star (here the equatorial conformal factor ${\varphi }_{r_{100}\mathrm{e}}$ at the 100th radial grid point) shows strong oscillations after core bounce. The time of bounce $t_{\mathrm{b}}$ is indicated by the vertical dotted line.

Figure 7

Influence of the density used in the wave extraction equations (upper panel) and of small differences in the initial model (lower panel) on the gravitational waveforms from rotational core collapse. If ${\rho }^{*}=\rho h{W}^2$ is used in the quadrupole formula (solid line) instead of $\rho$ (dashed line), the wave amplitude $A_{+}^{\mathrm{e}}$ increases by about 20% at core bounce (upper panel). A change from model SCC (solid line) to model $\mathrm{SC}{\mathrm{C}}_{\mathrm{SS}}$ (dashed line), which corresponds solely to a difference in the initial configuration, results in a qualitatively different waveform, in particular, during the ring-down phase (lower panel). The times of bounce $t_{\mathrm{b}}$ are indicated by the vertical dotted lines.

Figure 8

Influence of differences in the initial model on the evolution of the central density ${\rho }_{\mathrm{c}}$ for rotational core collapse. Changing from the collapse model SCC (solid line) to $\mathrm{SC}{\mathrm{C}}_{\mathrm{SS}}$ (dashed line) only slightly shifts the time of bounce $t_{\mathrm{b}}$ (indicated by the vertical dotted line), but leads to much stronger post-bounce ringdown oscillations. Nuclear matter density ${\rho }_{\mathrm{nuc}}$ is indicated by the horizontal dotted line.

Figure 9

Nonaxisymmetric density and velocity perturbation of the rapidly rotating neutron star equilibrium model RNS. By applying the perturbations described in the text, the original profiles (dashed lines) of the density $\rho$ along the azimuthal direction $\phi$ (upper panel), the radial velocity $v_r$ along the meridional direction $\theta$ (center panel), and the rotation velocity $v_{\phi }$ along the radial direction $r$ (lower panel) become strongly distorted (solid lines). The $\phi$-dependence of $\rho$ in the upper panel shows the nonaxisymmetric character of the perturbation.

Figure 10

Left (solid line) and right (dashed line) hand sides (computed on the finite difference grid) of the equation for the metric components $\varphi$ along the azimuthal direction $\phi$ (upper panel), ${\beta }^1$ along the meridional direction $\theta$ (center panel), and ${\beta }^3$ along the radial direction (lower panel). Even for strong nonaxisymmetric perturbations of the rotating neutron star model RNS, the metric solver 3 yields a highly accurate matching, such that the lines almost lie on top of one another. The insets show the relative difference ${\Delta }_{\mathrm{rel}u}$ between the left and right-hand sides of the equation for the same metric components. The relative differences are $\lesssim {10}^{-2}$, except where they exhibit a pole.

Figure 11

Time evolution of a symmetry violating perturbation. The upper two panels correspond to the spherically symmetric model SNS, and the lower two panels to the axisymmetric model RNS. The relative variation in density $\Delta \rho$ (solid line), radial velocity $\Delta v_r$ (dashed line), rotational velocity $\Delta v_{\phi }$ (dotted line), and conformal factor $\Delta \varphi$ (dashed-dotted line) show a remarkable constancy in time (note that $\Delta v_{\phi }$ is nonzero only for the rotating model RNS). The symmetry violating variation of the different fields scale with the initial perturbation amplitude (horizontal dotted lines; left panels: $a={10}^{-3}$; right panels: $a={10}^{-6}$).

Figure 12

Radial profile of the equatorial rotation velocity $v_{\phi \mathrm{e}}$ for the unperturbed axisymmetric rapidly rotating neutron star model RNS evolved in 3D. The profile of $v_{\phi \mathrm{e}}$ at $t=10\mathrm{ms}$ (dashed line) closely reproduces the initial profile (solid line). The stellar equatorial radius $r_{\mathrm{s}\mathrm{e}}$ and the boundary of the finite difference grid $r_{\mathrm{fd}}$ are indicated by vertical dotted lines.

Figure 13

Time evolution of the radial velocity at half the stellar equatorial radius $v_{r\mathrm{he}}$ for the perturbed rapidly rotating neutron star model RNS. The radial velocity shows regular oscillations with neither a noticeable drift nor damping when the 3D code is used in low resolution (solid line) as well as for the 2D code with high-resolution (dashed line). For comparison, the dashed-dotted line shows $v_{r\mathrm{he}}$ when no explicit perturbation is added. In this case the oscillations are triggered by truncation errors and (mostly) by the error resulting from using the CFC approximation in the evolution code.

Figure 14

Evolution of a strongly distorted nonaxisymmetric rotating neutron star model. The color coded distribution of $\log \rho$ on the equatorial plane shows how the initial perturbation (left panel) is partly accreted by the neutron star, and partly stretched into spiral arms (center panel). After about one rotation period of the neutron star, the trailing spiral arms have grown considerably in size (right panel). The rotation velocity $v_{\phi }$ is indicated by white arrows. Note that the atmosphere (color coded in black) has a density of much less than ${10}^7{\mathrm{g}\mathrm{cm}}^{-3}$, and that only the innermost 60 km of the computational domain are shown.

Figure 15

Time evolution of the central density ${\rho }_{\mathrm{c}}$ for distorted nonaxisymmetric rotating neutron star models. If the distortion is strong ($a=0.1$, solid line), matter accretion from the rotating bars results in a steep initial increase of ${\rho }_{\mathrm{c}}$, which slowly settles down to a new equilibrium state (indicated by the horizontal dotted line). For a small perturbation ($a=0.01$, dashed line), the evolution of ${\rho }_{\mathrm{c}}$ follows very closely that of an unperturbed model (dashed-dotted line).

Figure 16

Gravitational wave signal for distorted nonaxisymmetric rotating neutron star model. If the distortion is strong ($a=0.1$, solid lines), the nonzero gravitational wave amplitudes $A_{+}^{\mathrm{e}}$ (upper panel), $A_{+}^{\mathrm{p}}$ (center panel), and $A_{\times }^{\mathrm{p}}$ (lower panel) reach peak values of up to $\sim 15000\mathrm{cm}$. The amplitudes reduce significantly for $a=0.01$ (dashed lines). If an axisymmetric perturbation with $a=0.1$ is applied (dashed-dotted line), only the $A_{+}^{\mathrm{e}}$ gravitational wave mode is present.