Stable numerical simulation of Einstein equations
in gravitational collapse space–time

Numerical relativity, which solves the Einstein equations numerically, has been widely studied [1, 2, 3, 4]. In particular, for the direct observation of gravitational waves (e.g., [5]), numerical relativity makes important contributions.

Stable numerical simulations are important to clarify the details of phenomena. Since the Einstein equations are the nonlinear partial differential equations, numerical errors tend to accumulate. Although there are some research studies to reduce the numerical errors, they are mainly in the vacuum case. Thus, we suggest a modified system to reduce numerical errors by modifying the evolution equation in a perfect fluid.

The structure of this paper is as follows. We review the space–time decomposition of the Einstein equations in a perfect fluid in Sec. 2. In Sec. 3, we introduce the equations for the numerical simulations of the Einstein equations as the initial value problem in the perfect fluid. We perform some simulations in the dust case and propose a modified system to reduce the numerical errors in Sec. 4. We summarize this paper in Sec. 5. In this paper, indices such as $(\mu,\nu,\lambda,\cdots)$ and $(i,j,k,\cdots)$ run from 0 to 3 and 1 to 3, respectively. We use the Einstein convention of summation of repeated up–down indices.

The Einstein equations are as follows.

where $g_{\mu\nu}$ is the four-dimensional metric, ${}^{(4)}R_{\mu\nu}$ is the four-dimensional Ricci tensor, ${}^{(4)}R\equiv g^{\mu\nu}{}^{(4)}R_{\mu\nu}$, $g^{\mu\nu}$ is the inverse of $g_{\mu\nu}$, and $T_{\mu\nu}$ is the stress–energy tensor. Eq. (1) is not a dynamical form because the time and space components are mixed. Generally, we often carry out space–time decomposition of Eq. (1) for performing the simulations.

There are some dynamical forms of Eq. (1), the most basic formulation is the ADM formulation [6, 7]:

where $\gamma_{ij}\equiv g_{ij}$ is the induced metric, $\gamma^{ij}$ is the inverse of $\gamma_{ij}$, and $K_{ij}$ is the extrinsic curvature defined as Eq. (2). $\alpha\equiv 1/\sqrt{-g^{00}}$ and $\beta^{i}\equiv\alpha^{2}g^{0i}$ are the lapse function and the shift vector, respectively. $D_{i}$ is the covariant derivative associated with $\gamma_{ij}$, $R_{ij}$ is the Ricci tensor in three dimensions, $K\equiv\gamma^{ij}K_{ij}$, and $R\equiv\gamma^{ij}R_{ij}$. $\rho_{\mathrm{H}}\equiv\alpha^{2}g^{0\mu}g^{0\nu}T_{\mu\nu}$ is the mass density, $J_{i}\equiv\alpha g^{0\mu}T_{\mu i}$ is the momentum density, $S_{ij}\equiv T_{ij}$ is the stress tensor, and $S\equiv\gamma^{ij}S_{ij}$. $\mathcal{H}$ and $\mathcal{M}_{i}$ are the Hamiltonian constraint and the momentum constraint, respectively. The symbol $\approx$ means zero in the mathematical sense but nonzero in the numerical sense.

In the calculation in the nonvacuum case, we also solve the evolution equations of matter fields. These equations are given by the space–time decomposition of $g^{\mu\lambda}(\nabla_{\mu}T_{\lambda\nu})=0$, where $\nabla_{\mu}$ is the covariant derivative associated with $g_{\mu\nu}$. For the perfect fluid case, the stress–energy tensor is given as

where $\rho$ is the rest mass, $\epsilon$ is the inner energy, $u^{\mu}$ is the four velocity, and $p$ is the pressure. The evolution equations of $\epsilon$ and $\rho u^{i}$ are given by $g^{\mu\lambda}(\nabla_{\mu}T_{\lambda\nu})=0$. The evolution equation of $\rho$ is given by the continuous equation $\nabla_{\mu}(\rho u^{\mu})=0$. On the other hand, $p$ is usually given by other conditions such as the equation of state.

In the dust case, we solve the Einstein equations and the matter evolution equations. With reference to [14], we set the initial data as follows.

where $a$ is the gravitational mass $M_{g}=\int d^{3}x(\rho_{\text{H}}\psi^{5}+(1/16\pi)\psi^{-7}\hat{A}_{ij}\hat{A}_{mn}\hat{\gamma}^{im}\hat{\gamma}^{jn})$ as the unit at the initial time and $r=\sqrt{x^{2}+y^{2}+z^{2}}$. This time, we set $a=4.129\times 10^{-3}$, $\delta r^{2}=0.18M_{g}^{2}$, and $r_{0}=4.0M_{g}$. The numerical ranges are $-9\leq x,y,z\leq 9$. This case is static at the initial time, so we set $X_{i}=\hat{A}_{ij}=J_{i}=\beta^{i}=u^{i}=0$. In the dust case, we set the inner energy as $\epsilon=0$ and the pressure $p=0$. We assume $\hat{\gamma}_{ij}=\delta_{ij}$ and $K=0$. Then, the extrinsic curvature $K_{ij}=0$ and $\hat{R}=0$. We obtain the initial data $(\alpha,\gamma_{ij})$ with the above conditions by solving Eqs. (7) and (9). With reference to the Schwarzschild metric, the initial step values of $\psi$ and $\alpha$ are set as $1+1/(1+6r)$ and $1-1/(3/2+r)$, respectively.

Recently, the maximal slicing condition, Eq. (9), and the minimal distortion gauge condition are almost never used in the evolution because the numerical costs are high. The evolution equation of $\alpha$ often uses the $1+\log$ slicing condition [15],

Since this condition has the characteristic of singularity avoidance, it is widely used in the simulations for the gravitational collapse models. We choose $\beta^{i}=0$. For the rotating models of neutron stars and/or black holes, the Gamma-driver condition [16] is often used in the evolution equations of $\beta^{i}$.

This time, we select the grid as $\Delta x=\Delta y=\Delta z=1/60$, $\Delta t=1/240$, and the boundary condition as the approximate asymptotic flat boundary. We use the fourth-order Runge–Kutta scheme with mainly the second-order centered space difference. However, for only the advection term, which is the first term on the right-hand side in each of Eqs. (21)–(23), we use the second-order upwind scheme. The direction is defined by the signature of the velocity $v^{i}$.

Fig. 1 shows the weighted rest mass $\rho_{*}$ at $t=12$ and $t=13$ obtained by solving Eqs. (10)–(14), Eqs. (21)–(23), and Eq. (28). We see that $\rho_{*}<0$ at $t=13$ in the $-2<z<2$ range. We check $\rho_{*}\geq 0$ during $0\leq t\leq 12$ in all the simulation ranges. Since $\rho_{*}\geq 0$ is a necessary condition for successful simulations, this simulation fails after $t=13$.

For more stable simulations, we modify Eq. (23) as

where $\kappa$ is a damping parameter. This modification is based on the following ideas. By using this modification, we obtain the following: (i) the positive rest mass condition $\rho_{*}\geq 0$ is satisfied, (ii) the constraints, Eqs. (16)–(20), are conserved, and (iii) the total weighted rest mass

is conserved. The negative sign of $\kappa$ makes the simulations stable because the dynamical equations of $\tilde{\mathcal{M}}_{i}$ become

because of the adjusted terms of Eq. (29). The adjusted term of Eq. (31) has the dissipation effect if $\kappa<0$ because $\rho_{*}/(e^{6\varphi}w)\geq 0$. The set of Eqs. (10)–(14), Eqs. (21)–(22), Eq. (28), and Eq. (29) is called as the modified system hereafter.

We show the numerical results of $\rho_{*}$ obtained using the modified system in Fig. 2. The numerical settings are consistent with those in Fig. 1 without the evolution equation of $\rho_{*}\hat{u}_{i}$ and we set $\kappa=-0.1$. For this simulation, we check $\rho_{*}\geq 0$ during $0\leq t\leq 17$ in all simulation ranges. In Fig. 3, we show the L2 norm of the following values:

We see that the norm of the constraints $\mathcal{C}^{2}$ with $\kappa=-0.1$ is less than those in other cases until $0\leq t\leq 17$. On the other hand, the norm with $\kappa=0.1$ is larger than those in other cases. Thus, these results are consistent with the analytical results using Eq. (31). Fig. 4 shows the relative errors against the initial values of the total weighted mass $I$ in Eq. (30). They are not markedly different between the cases of $\kappa=0.0$ and $\kappa=-0.1$.

We reviewed the equations for construction of the initial values, the dynamical equations of the Einstein equations called the BSSN formulation, and the dynamical equations of the matter fields in the perfect fluid. With these equations, we performed the simulations in the dust case. We modified the system by modifying the evolution equations of the matter field, investigated the stability analytically, and performed the simulations with the modified system to confirm the consistency of the analytical results. In addition, the lifetime of the simulations was extended from $t=12$ to $t=17$ with the modification.

T.T. was partially supported by JSPS KAKENHI Grant Number 21K03354 and a Grant for Basic Science Research Projects from The Sumitomo Foundation. G.Y. and T.T. were partially supported by JSPS KAKENHI Grant Number 20K03740. G.Y. was partially supported by a Waseda University Grant for Special Research Projects 2021C-138.