The notation is complicatedly mixed from Refs. [1-3]. The relations are written as much as possible.

Fluid quantities and dynamical eqs in 3+1 form

We consider the perfect fluid whose energy-momentum tensor is given by

$$ T_{\mu\nu}=(\rho + P)u_\mu u_\nu +P g_{\mu\nu}, $$

where $u^\mu$ is the fluid four-velocity and the Lorentz factor $\Gamma$ is given by [^1]

$$ \Gamma=-u^\mu n_\mu $$

with $n_\mu$ being the normal one-form for the time slice. We also introduce the velocity relative to the Eulerian observer $U^\mu$ as [^2]

$$ u^\mu=\Gamma(n^\mu+U^\mu) $$

with $n^\mu U_\mu=0$ and

$$ \begin{equation} \Gamma=\left(1-U^iU_i\right)^{-1/2}. \quad (1) \end{equation}$$

The fluid energy density measured by the Eulerian observer is given by

$$ \begin{equation} E=T^{\mu\nu}n_\mu n_\nu=\Gamma^2(\rho+P)-P. \quad (2) \end{equation} $$

Let us write the proper rest mass energy density as $\rho_0$. Then, the relativistic specific enthalpy is defined by [^3]

$$ h=\frac{\rho+P}{\rho_0}=\frac{\rho_0(1+\varepsilon)+P}{\rho_0}, $$

where $\rho_0\varepsilon$ is the internal energy. Let $D$ denote the baryon rest mass density measured by the Eulerian observer as [^4]

$$ D=\rho_0\Gamma. $$

For later convenience, we also introduce the fluid momentum density measured by the Eulerian observer $p_\mu$ as

$$ p_\mu=(E+P)U_\mu. $$

Then, the fluid equations are written as [1]

$$\begin{align} \left(\partial_t-\mathcal L_\beta\right)D+D_i(\alpha D U^i)-\alpha K D=0, \ \left(\partial_t-\mathcal L_\beta\right)E+\alpha\left[D_ip^i-(E+P)(K+K_{ij}U^iU^j)\right] +p^iD_i\alpha=0, \ \left(\partial_t-\mathcal L_\beta\right)p_i+ \alpha D_j(P\delta^j_i+p_iU^j)+\left[P\delta^j_i+p_iU^j\right] D_j\alpha -\alpha K p_i+ED_i\alpha=0. \end{align} $$

These equations can be rewritten as

$$\begin{align} \partial_t(\sqrt{\gamma}D)+\partial_i\left[ \alpha\sqrt{\gamma}D\left( U^i-\frac{\beta^i}{\alpha} \right) \right]=0, \ \partial_t(\sqrt{\gamma}E)+\partial_i\left[ \alpha\sqrt{\gamma}\left( p^i-\frac{\beta^i}{\alpha}E \right) \right]+\sqrt{\gamma}\left( p^i\partial_i\alpha-\alpha S_{ij}K^{ij} \right)=0, \ \partial_t(\sqrt{\gamma}p_i)+\partial_j \left[ \alpha\sqrt{\gamma}\left( p_i\left(U^j-\frac{\beta^j}{\alpha}\right)+\delta^j_iP \right) \right]\cr +\sqrt{\gamma}\left( E\partial_i\alpha-p_j\partial_i\beta^j+\frac{1}{2}\alpha S_{jk}\partial_i\gamma^{jk} \right)=0, \end{align}$$

where

$$ S_{ij}=(E+P)U_iU_j+P\gamma_{ij}. $$

Let us define the following variables:

$$\begin{align} \rho_*&=&\sqrt{\gamma}D,\ S_0&=&\sqrt{\gamma}E,\ S_i&=&\sqrt{\gamma}p_i. \end{align}$$

Then, the equations can be rewritten as

$$ \begin{align} \partial_t \rho_{*}&+\partial_i\left[\rho_{*}V^i\right]=0, \quad (3) \ \partial_t S_0&+\partial_i\left[S_0V^i+P\sqrt{\gamma}(V^i+\beta^i)\right] +S^iD_i\alpha-\alpha \sqrt{\gamma}S_{ij}K^{ij}=0, \ \partial_tS_i&+\partial_j \left[S_iV^j+\alpha\sqrt{\gamma}\delta^j_iP\right]+S_0D_i\alpha-S_j\partial_i\beta^j+\frac{1}{2}\alpha \sqrt{\gamma}S_{jk}\partial_i\gamma^{jk}=0, \end{align} $$

where we have introduced $V^i$ as

$$ V^i=\alpha U^i-\beta^i, $$

and used the following relations

$$\begin{equation} \alpha p^i-E\beta^i=EV^i+P(V^i+\beta^i). \end{equation}$$

Note that $V^i=u^i/u^0$. [^5] We solve the above equations for the dynamical variables $\rho_{*}$, $S_0$ and $S_i$. The so-called primitive variables are $\rho$, $V^i$ and $\varepsilon$. The fluxes are given by

$$\begin{align} f_{\rho_{*}}^i&=\rho_{*}V^i, \ f_{S_0}^i&=S_0V^i+\sqrt{\gamma}P(V^i+\beta^i), \ f_{S_j}^i&=S_jV^i+\alpha\sqrt{\gamma}\delta^i_jP. \end{align} $$

We note that the variable $\gamma$ in the fluxes can be evaluated by $\gamma={\rm e}^{12\psi}\tilde \gamma ={\rm e}^{12\psi}f$ with $f$ being the determinant of the reference flat metric.

From the Jacobi matrix of the fluxes, we obtain the following expression for the three characteristic speeds for three directions:

$$\begin{align} \lambda_0^i&=V^i,\ \lambda_\pm^i&= \frac{\alpha}{1-U^2c_{\rm s}^2}\left(U^i(1-c_{\rm s}^2) \pm c_{\rm s}\sqrt{(1-U^2)\left[\gamma^{ii}(1-U^2c_{\rm s}^2)-(1-c_{\rm s}^2)U^iU^i\right]}\right), \end{align}$$

where $c_{\rm s}$ is the sound velocity defined by

$$\begin{equation*} c_{\rm s}^2=\left(\frac{\partial P}{\partial\rho}\right)_s \end{equation*}$$

with fixed specific entropy $s$.

Primitive variables from dynamical variables

The equations between the primitive variables and the conserved variables are given by

$$\begin{align} P&=P(\rho,\varepsilon), \nonumber \ \rho_*&=\sqrt{\gamma}\Gamma\frac{\rho}{1+\varepsilon}, \nonumber \ S_0&=\sqrt{\gamma}\left[\Gamma^2(\rho+P)-P\right], \nonumber \ S_i&=\sqrt{\gamma}\left(E+P\right)U_i =\frac{1}{\alpha}\left(S_0+\sqrt{\gamma}P\right)\gamma_{ij}\left(V^j+\beta^j\right). \quad (4) % \label{eq:ptod4} \end{align}$$

We need to invert these equations to obtain the primitive variables from the dynamical variables.

From Eqs. (1) and (3), we obtain

$$\begin{equation} \Gamma^2p^2=(E+P)^2(\Gamma^2-1), \quad (5) % \label{eq:Gam2p2} \end{equation}$$

where $p^2=p^\mu p_\mu$. From Eq. (2), we find

$$\begin{equation} \Gamma^2=\frac{E+P}{\rho+P}. \quad (6) % \label{eq:Gam2} \end{equation}$$

Substituting Eq. (6), into Eq. (5), we obtain

$$ \begin{equation} p^2-E^2-(P-\rho)E+\rho P=0. \quad (7) \end{equation}$$

Barotropic EoM

Let us assume $P=P(\rho)$. In this case, we need not solve the continuity equation (3) unless we are interested in the value of the internal energy. Since the values of $E$ and $p_i$ can be calculated from the dynamical variables $S_0$ and $S_i$, Eq. (7) can be regarded as an equation to get the value of $\rho$ with a given equation of state $P=P(\rho)$. Once the value of $\rho$ is calculated, the value of $\Gamma$ is given by Eq. (6). From Eq. (4), the value of $V^i$ is given by

$$\begin{equation*} V^i=\alpha U^i-\beta^i=\alpha\frac{\gamma^{ij}S_j}{S_0+\sqrt{\gamma}P}-\beta^i. \end{equation*}$$

The value of $\varepsilon$ is given by

$$ \begin{equation} \varepsilon=\frac{\rho-\rho_0}{\rho_0}=\frac{\sqrt{\gamma}\Gamma\rho-\rho_{*}}{\rho_{*}}. \quad (8) % \label{eq:varep} \end{equation} $$

For the case $P=w\rho$, we obtain

$$ \begin{equation*} \rho=\frac{1}{2w}\left[-(1-w)E+\sqrt{E^2(1-w)^2+4w(E^2-p^2)}\right]. \end{equation*} $$

From Eq. (6), the value of $\Gamma$ is given by

$$\begin{equation*} \Gamma^2=\frac{E+w\rho}{(1+w)\rho}. \end{equation*}$$

From the conservation equation

$u^\mu\nabla_\nu T_\mu^\nu=0$ and the continuity equation $\nabla_\mu(\rho_0 u^\mu)=0$, we obtain

$$\begin{equation*} \mathrm{d}\ln(1+\varepsilon)=w\mathrm{d}\ln\rho_0. \end{equation*}$$

Then, the value of $\varepsilon$ is given by

$$\begin{equation} \varepsilon=C\rho_0^w-1, \end{equation}$$

or, equivalently,

$$\begin{equation} \varepsilon=\varepsilon_\rho:=C'\rho^{w/(1+w)}-1, % \label{eq:vareps} \quad(9) \end{equation}$$

where $C$ and $C'$ are integration constants fixed by the initial condition. If we solve Eq. (3), the value of $\varepsilon$ calculated by Eq. (8) should be consistent with Eq. (9).

Flux calculation scheme

For the flux calculation, we use a central scheme with MUSCL (Mono Upstream-centered Scheme for Conserva-tion Laws) [4,5].

For simplicity, let us consider the following 1-dimensional equation in a conserved form:

$$\begin{equation*} \partial_t u +\partial_x f(u)=0. \end{equation*}$$

We evaluate the flux in a finite difference formula through the values at mid-points between two grid points as follows:

$$\begin{equation*} \partial_t u=-\frac{1}{\Delta x}\left(f_{i+1/2}-f_{i-1/2}\right). \end{equation*}$$

In order to evaluate $f_{i\pm1/2}$, we first introduce the following variables for the value of $u$ at the point specified by $i+1/2$:

$$\begin{align} (u_L){i+1/2}&=u_i+\frac{1}{4}\left((1-\kappa)(\bar \Delta-)i+(1+\kappa)(\bar \Delta+)i\right),\ (u_R){i+1/2}&=u_{i+1}-\frac{1}{4}\left((1-\kappa)(\bar \Delta_+){i+1}+(1+\kappa)(\bar \Delta-)_{i+1}\right), \end{align}$$

where $\kappa$ is a parameter to specify the way to the interpolation, and $(\bar \Delta_\pm)_i$ is defined by

$$\begin{align} (\bar \Delta_+)i&=&{\rm minmod}((\Delta+)i,b(\Delta-)i),\ (\bar \Delta-)i&=&{\rm minmod}((\Delta-)i,b(\Delta+)_i) \end{align}$$

with

$$\begin{align} (\Delta_+)i&=u{i+1}-u_i,\ (\Delta_-)i&=u_i-u{i-1}. \end{align}$$

The function minmod($a$,$b$) returns 0 if $a$ and $b$ have different signs, and if not, the value of one of the arguments which has a smaller absolute value. That is,

$$\begin{equation*} {\rm minmod}(a,b)={\rm sign}(a){\rm max}(0,{\rm min}(|a|,{\rm sign}(a) b)). \end{equation*}$$

Then, the flux is evaluated as

$$\begin{equation} f=\frac{1}{2}\left(f(u_L)+f(u_R)-a^{*}(u_R-u_L)\right), \end{equation} $$

where $a^*$ is the value of the fastest characteristic speed given by

$$\begin{equation*} a^{*}={\rm max}{|\lambda_{0L}|,|\lambda_{+L}|,|\lambda_{-L}|,|\lambda_{0R}|,|\lambda_{+R}|,|\lambda_{-R}|}. \end{equation*}$$

::: thebibliography

[1]: E. Gourgoulhon, (Springer, Berlin New York, 2012).

[2]: J. A. Font, Living Rev. Rel. 11, 7 (2008), Numerical Hydrodynamics and Magnetohydrodynamics in General Relativity.

[3]: M. Shibata, (World Scientific Publishing Co. Pte. Ltd, Singapore, 2016).

[4]: A. Kurganov and E. Tadmor, Journal of Computational Physics 160, 241 (2000), New High-Resolution Central Schemes for Nonlinear Conservation Laws and Convection-Diffusion Equations.

[5]: M. Shibata and J. A. Font, Phys. Rev. D72, 047501 (2005), arXiv:gr-qc/0507099, Robustness of a high-resolution central scheme for hydrodynamic simulations in full general relativity. :::

[^1]: $\Gamma$ is $W$ in Font [2] and $w$ in Shibata [3].

[^2]: $U^i$ is $v^i$ in Font [2].

[^3]: $\rho_0$ is $\rho$ in Shibata [3].

[^4]: $D$ is $m_{\rm B}\mathcal N_{\rm B}$ in Gourgoulhon [1].

[^5]: $V^i$ is $v^i$ in Shibata [3].