Derivation Notes on viscous EoM in SPheRIO

Derivation Notes on viscous equation of motion in SPheRIO

Last updated 25 May, 2013

Summary
Here we derive the equations of motion (EoM) with shear and bulk viscosity in hyperbolic coordinate in terms smoothed particle hydrodynamics (SPH) degree of freedom. This note is based on previous notes of Gabriel, Philipe and Jaque, however, some misprints were found in the original documents, corrections were made accordingly in the present version.

Hyperbolic coordinates
Some definitions and identities on covariant derivative and metric.


 * $$\begin{align}

&\Gamma^i_{k\ell}=\frac{1}{2}g^{im} \left(\frac{\partial g_{mk}}{\partial x^\ell} + \frac{\partial g_{m\ell}}{\partial x^k} - \frac{\partial g_{k\ell}}{\partial x^m} \right) = {1 \over 2} g^{im} (g_{mk,\ell} + g_{m\ell,k} - g_{k\ell,m}) \\ &\Gamma^i_{ki}=\frac{1}{2}g^{im}\frac{\partial g_{im}}{\partial x^k}={1 \over 2g}\frac{\partial g}{\partial x^k}=\frac{\partial \log \sqrt{|g|}}{\partial x^k}\\ &D_{\mu}\phi\equiv(\phi)_{;\mu}=\partial_{\mu}\phi\\ &D_{\mu}V_{\nu}\equiv(V_{\nu})_{;\mu}=\partial_{\mu}V_{\nu}-\Gamma^{\alpha}_{\mu\nu}V_{\alpha}\\ &D_{\mu}V^{\nu}\equiv(V^{\nu})_{;\mu}=\partial_{\mu}V^{\nu}+\Gamma^{\nu}_{\mu\alpha}V^{\alpha}\\ &D_{\mu}V_{\alpha\beta}\equiv(V_{\alpha\beta})_{;\mu}=\partial_{\mu}V_{\alpha\beta}-\Gamma^{\lambda}_{\alpha\mu}V_{\lambda\beta}-\Gamma^{\lambda}_{\beta\mu}V_{\lambda\alpha}\\ &D_{\mu}V^{\mu}\equiv(V^{\mu})_{;\mu}=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}V^{\mu})\\ &D_{\mu}V^{\mu\nu}\equiv(V^{\mu\nu})_{;\mu}=\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}V^{\mu\nu})+\Gamma^{\nu}_{\alpha\mu}V^{\alpha\mu} \end{align}$$

The hyperbolic coordinates are defined as follows


 * $$\begin{align}

&\tau=\sqrt{t^2-z^2}\\ &\eta=\frac{1}{2}\ln\left(\frac{t+z}{t-z}\right) \end{align}$$

Here we briefly list the metric of the coordinates sistem and some useful expressions


 * $$\begin{align}

g_{\mu\nu}=\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & -1 & 0 & 0\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -{\tau^2} \end{pmatrix} \end{align}$$


 * $$\begin{align}

\sqrt{-g}=\tau \end{align}$$

The inverse of the metric reads,


 * $$\begin{align}

g^{\mu\nu}=\begin{pmatrix} 1 & 0 & 0 & 0\\ 0 & -1 & 0 & 0\\ 0 & 0 & -1 & 0\\ 0 & 0 & 0 & -\frac{1}{\tau^2} \end{pmatrix} \end{align}$$

One obtains the Christoff symbols


 * $$\begin{align}

\Gamma^i_{k\ell}=\frac{1}{2}g^{im} \left(\frac{\partial g_{mk}}{\partial x^\ell} + \frac{\partial g_{m\ell}}{\partial x^k} - \frac{\partial g_{k\ell}}{\partial x^m} \right) \end{align}$$

It is straightforward to show that only three terms are non-vanishing, they are


 * $$\begin{align}

\Gamma^0_{33}=\tau, \Gamma^3_{30}=\Gamma^3_{03}=\frac{1}{\tau} \end{align}$$

In the follows, we will consider the second order viscous EoM in hyperbolic coordinates.

SPH degree of freedom
In practice, one need to express everything in SPH degree of freedom. By the SPH algorithm, the density $$A$$ at point $$j$$ in the laboratory frame is obtained by the average among SPH particles at points $$i $$ in terms of its density ratio to that of a conserved charge $$\nu$$ in the laboratory frame.


 * $$\begin{align}

A_j \equiv\sum_i \left(\frac{A}{\sigma^*}\right)_i\nu_iW_{ji} =\sum_i \left(\frac{A}{\sigma\tau\gamma}\right)_i\nu_iW_{ji} \end{align}$$

Note that if $$A$$ itself is the density of a conserved charge $$B$$, namely, $$b^*$$, one has


 * $$\begin{align}

&A=b^* \\ &\left(\frac{A}{\sigma^*}\right)_i=\left(\frac{b^*}{\sigma^*}\right)_i=\left(\frac{b}{\sigma}\right)_i \\ &A_j=\sum_i \left(\frac{A}{\sigma^*}\right)_i\nu_iW_{ji}=\sum_i b_iW_{ji} \end{align}$$

For other quantities, such as pressure, one general has to compare its value to the density of a conserved charge by using Equation of State (EoS). In practice, the EoM may involve the spatial derivative of a quantity, usually it is


 * $$\begin{align}

(\partial A)_j\equiv\sum_i \left(\frac{A}{\sigma^*}\right)_i\nu_i\nabla_jW_{ji} \end{align}$$ However, in case of pressure $$\partial_i {\tilde P}$$, one has to consider a symmetric form (due to the third law of Newton). The following case will be encountered in the EoM in practice. This in part shows the complication of SPH method itself.


 * $$\begin{align}

&(\partial \tilde P)_j= \partial \left(\frac{\tilde P}{\phi} \phi\right)_j = \phi_j\partial \left(\frac{\tilde P}{\phi} \right)_j+\left(\frac{\tilde P}{\phi}\right)_j\partial \left( \phi\right)_j \\ &= \phi_j \sum_i \left(\frac{\tilde P}{\phi\sigma^*}\right)_i\nu_i\nabla_jW_{ji}+\left(\frac{\tilde P}{\phi}\right)_j \sum_i \left(\frac{\phi}{\sigma^*}\right)_i\nu_i\nabla_jW_{ji}\\ &= \sigma^*_j \sum_i \left(\frac{\tilde P}{\sigma^*\sigma^*}\right)_i\nu_i\nabla_jW_{ji}+\left(\frac{\tilde P}{\sigma^*}\right)_j \sum_i \left(\frac{\sigma^*}{\sigma^*}\right)_i\nu_i\nabla_jW_{ji}\\ &=\sum_i \sigma^*_j\frac{\tilde P_i}{{\sigma^*_i}^2}\nu_i\nabla_jW_{ji}+\sum_i \frac{\tilde P_j}{\sigma^*_j}\nu_i\nabla_jW_{ji}\\ &= \sum_i \nu_i\sigma^*_j\left(\frac{\tilde P_i}{{\sigma^*_i}^2}+ \frac{\tilde P_j}{{\sigma^*_j}^2}\right)\nabla_jW_{ji} \end{align}$$

We first introduce a density $$ \phi $$, then take it as $$\phi=\sigma^*$$ at the end, in short


 * $$\begin{align}

(\partial {\tilde P})_i = \sum_j \nu_j \sigma_i^* \left(\frac{{\tilde P}_i}{\sigma_i^{*2}}+\frac{{\tilde P}_j}{\sigma_j^{*2}}\right) \bigtriangledown_i W_{ij} \end{align}$$

There is another term involving spatial derivative, it is $$\partial_i v^i $$. We make use of the conserved current of $$\sigma u$$,


 * $$\begin{align}

\left(\sigma u^{\mu}\right)_{;\mu} =0 \end{align}$$

The l.h.s. can be written as


 * $$\begin{align}

&\partial_{\mu} (\sigma^* \frac{u^{\mu}}{u^0}) =\partial_{\tau}\sigma^* + \partial_i (\sigma^*   v^i) = 0   \\ &\partial_{\tau}\sigma^* + \sigma^* \partial_i v^i = 0 \end{align}$$

The second term on the r.h.s. contains what we need, we note


 * $$\begin{align}

\sigma^*_j =\sum_i {\nu}_i W_{ji} \end{align}$$

Now considering the l.h.s., for a given point $$ j$$, one has


 * $$\begin{align}

\partial_{\tau}\sigma^*_j = \sum_i {\nu}_i \partial_{\tau}W_{ji}(|r_j-r_i|)  =\sum_i{\nu}_i \left[\dot{r_j}\cdot \nabla_jW_{ji}(|r_j-r_i|)+\dot{r_i}\cdot\nabla_iW_{ji}(|r_j-r_i|)\right]   =\sum_i  {\nu}_i (\dot{r_j}-\dot{r_i}) \cdot  \nabla_jW_{ji}(|r_j-r_i|)) \end{align}$$

where


 * $$\begin{align}

\nabla_iW_{ji}(|r_j-r_i|)=-\nabla_jW_{ji}(|r_j-r_i|) \end{align}$$

We did not put a arrow on top of vectors, to be specific


 * $$\begin{align}

&W_{ji}\equiv W_{ji}(|r_j-r_i|)\equiv  W_{ji}(|\vec r_j-\vec r_i|) \\ &\nabla_jW_{ji}\equiv \nabla_jW_{ji}(|r_j-r_i|)\equiv \frac{\partial W_{ji}(|\vec r_j-\vec r_i|)}{\partial \vec r_j}   \\ &\nabla_iW_{ji}\equiv \nabla_iW_{ji}(|r_j-r_i|)\equiv \frac{\partial W_{ji}(|\vec r_j-\vec r_i|)}{\partial \vec r_i} \end{align}$$

Finally one obtains


 * $$\begin{align}

\left(\partial_i v^i\right)_j = - \frac{1}{\sigma^*}\sum_k \nu_k (\dot r_j - \dot r_k) \cdot  \nabla_j W_{jk} \end{align}$$

Symbols and basic relations
We will now define some symbols for furture use. But first, some comments on notations. The latter $$D$$ is reserved for covariant derivative, other symbols are understood as "normal" derivatives. The Greek letter $$\tau$$ indicates (for most of the time) the time-like coordinate in hyperbolic coordinates, rather than the proper time. An important exception is that $$D\tau$$ indicates the derivation with respect to proper time, rather than the time-like coordinate, the latter is denoted by $$d\tau$$. The symbol $$\partial$$ when by itself, usually implies "covariant", since this comes from a straightforward generalization of the formulae in flat space. In this context, we define the total derivative of proper time as follows


 * $$\begin{align}

\frac{D}{D\tau}\equiv u^{\mu}D_{\mu} \end{align}$$

This implies the following substitution


 * $$\begin{align}

&D\tau\longleftrightarrow \frac{1}{\gamma}d\tau\\ &\frac{D}{d\tau}\longleftrightarrow \frac{1}{\gamma}u^{\mu}D_\mu \end{align}$$

In other cases, $$\tau$$ means time-like coordinate, for instance,


 * $$\begin{align}

d_\tau f \equiv \frac{df}{d\tau} =\frac{1}{\gamma}u^\mu\partial_\mu f \end{align}$$

Make use of the projection operator of 4-velocity $$\Delta^{\mu\nu}\equiv g^{\mu\nu}-u^\mu u^\nu $$, we can write anything (in this case, the covariant derivative) as a summation of parallel and perpendicular parts:


 * $$\begin{align}

D^\mu=g^{\mu\nu}D_{\nu}=\Delta^{\mu\nu}D_\nu+u^\mu\frac{D}{D\tau}\equiv D^\mu_\perp+u^\mu\frac{D}{D\tau} \end{align}$$

Make use of $$\Delta_{\nu\beta}u^\beta=0,u_\beta u^\beta=1$$ andintegral in parts, as well as $$D_\mu g_{\alpha\beta}=0$$, one finds


 * $$\begin{align}

\Delta_{\nu\beta}D^\alpha u^\beta=-u^\beta D^\alpha\Delta_{\nu\beta}=D^\alpha u_\nu \end{align}$$

therefore


 * $$\begin{align}

\Delta_{\mu\alpha}\Delta_{\nu\beta}D^{\alpha} u^{\beta}=\Delta_{\mu\alpha}D^{\alpha} u_{\nu}=D_{\mu\perp} u_{\nu}=D_{\mu} u_{\nu}-u_{\mu}\frac{Du_{\nu}}{D\tau} \end{align}$$

Similarly, since the shear is perpendicular to the 4-velocity $$\pi^{\alpha\beta}u_alpha=\pi^{\alpha\beta}u_beta=0$$, one fines


 * $$\begin{align}

&\pi^{\alpha\beta}\frac{D}{D\tau}\Delta_{\alpha\beta}=0 \end{align}$$

The following expression is useful, we note that it is only correct in hyperbolic coordinates


 * $$\begin{align}

\Gamma^\lambda_{\mu\nu}u_\lambda=\Gamma^0_{\mu\nu}\gamma+\Gamma^3_{\mu\nu}u_3 =\tau\delta_{\mu 3}\delta_{\nu 3}\gamma+\frac{u_3}{\tau}(\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}) \end{align}$$

One can futher show


 * $$\begin{align}

u^{\mu}\Gamma^{\lambda}_{\mu\nu}u_{\lambda} =-\frac{1}{\tau^3}(u_3)^2\delta_{\nu 0} = -{\tau}(u^3)^2\delta_{\nu 0} \end{align}$$

We now introduce some notations, which will be used below


 * $$\begin{align}

&\partial \cdot u \equiv \left(u^{\mu}\right)_{;\mu} \\ &= \frac{1}{\tau} \partial_{\mu} \left(\tau u^{\mu}\right)  \\ &= \frac{1}{\tau} \left(\partial_{\tau} \left(\tau \gamma\right)+\partial_i\left(\tau\gamma v^i \right)\right)   \\ &= \frac{1}{\tau} \left(\partial_{\tau} \left(\tau \gamma\right)+v^i \partial_i\left(\tau\gamma  \right)\right)+\frac{1}{\tau}(\tau\gamma)\partial_i v^i   \\ &= \frac{1}{\tau} d_{\tau} (\tau\gamma) + \gamma \partial_i v^i  \\ &= \frac{\gamma}{\tau} + d_{\tau}\gamma+\gamma \partial_i v^{i} \end{align}$$

For instance, the time derivative of conserved charge can be written as


 * $$\begin{align}

&0= D_\mu(\sigma u^{\mu})=\sigma(\partial\cdot u)+\gamma d_{\tau}\sigma \\ &d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \end{align}$$

We can further express the r.h.s. of the above equation in terms of the final variables


 * $$\begin{align}

&d_{\tau}\gamma = \frac{d}{d\tau}\left(1-g_{ij}u^iu^j\right)^{1/2}  \\ &= -\frac{g_{ij}u^i}{\gamma}\frac{du^j}{d\tau}-\frac{u^iu^j}{2\gamma}\frac{dg_{ij}}{d\tau}  \\ &= -\frac{u_j}{\gamma}\frac{du^j}{d\tau}+\frac{\left(u^3\right)^2\tau}{\gamma} \end{align}$$

we note


 * $$\begin{align}

&u_\mu\partial^\mu u_0=u_\mu\partial^\mu\gamma=\gamma d_\tau\gamma\\ &\frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i} \end{align}$$

For the spatial derivative of velocity, one has


 * $$\begin{align}

&\partial_k(u^\mu u_\mu)=0\\ &0=u^\mu\partial_ku_\mu=\gamma\partial_k\gamma+u^i\partial_ku_i\\ &\partial_k\gamma=-v^i\partial_ku_i\\ &0=\partial_k\gamma+v^i\partial_ku_i=\partial_k\gamma+v^iv_i\partial_k\gamma+u^i\partial_kv_i=\partial_k\gamma+\frac{(1-\gamma^2)}{\gamma^2}\partial_k\gamma+u^i\partial_kv_i\\ &\partial_k\gamma=-{\gamma^2}u^i\partial_kv_i\\ \end{align}$$

The second relation above may introduce some numerical error, so in practice we will use the first one, which involves the spacial derivative of 4-velocity.

The proper time derivative of the spatial part of 4-velocity can be evaluated as follows


 * $$\begin{align}

&u\cdot \partial u_{i} \equiv \gamma \frac{du_i}{d\tau} - u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu}  \\ &= \gamma \frac{du_i}{d\tau}   \\ &= - \gamma\frac{du^1}{d\tau}\delta_{i1}- \gamma\frac{du^2}{d\tau}\delta_{i2}- \left(\gamma\tau^2\frac{du^3}{d\tau}+2\tau\gamma u^3 \right)\delta_{i3} \end{align}$$


 * $$\begin{align}

\gamma\frac{du_i}{d\tau}=\gamma\frac{du_1}{d\tau}\delta_{i1}+\gamma\frac{du_2}{d\tau}\delta_{i2}+\gamma\frac{du_3}{d\tau}\delta_{i3} \end{align}$$

where in the second equality one makes use of


 * $$\begin{align}

u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} = 0 \end{align}$$

One finds the following relation


 * $$\begin{align}

u\cdot \partial u_i \equiv u_{\mu}{(u_i)_;}^{\mu} = u_{\mu}{(u_i)_,}{^{\mu}} \end{align}$$

Now, in the hyperbolic coordinates, one can make use of the specific form of the metric


 * $$\begin{align}

\frac{\partial}{\partial x^{\beta}}\equiv\partial_\beta =g_{\mu\alpha}\frac{\partial}{\partial x_\alpha}=g_{\mu\alpha}\partial^\alpha \end{align}$$

and


 * $$\begin{align}

\partial^\beta = g^{\mu\alpha}\partial_\alpha \end{align}$$

Note that in a general coordinate system, covariant derivative is different from "normal" derivative, so one cannot change the order of subscript and superscript, since


 * $$\begin{align}

&\partial^\mu \equiv \frac{\partial}{\partial x_\mu}=\frac{\partial}{\partial(g_{\mu\alpha}x^\alpha)}=\frac{\partial x^\beta}{\partial(g_{\mu\alpha}x^\alpha)}\frac{\partial}{\partial x^\beta} \\ &\partial_\beta \equiv \frac{\partial}{\partial x^\beta}=\frac{\partial(g_{\mu\alpha}x^\alpha)}{\partial x^\beta}\frac{\partial}{\partial x_\mu} =(g_{\mu\beta}+g_{\mu\alpha,\beta}x^\alpha)\frac{\partial}{\partial x_\mu} \end{align}$$

therefore


 * $$\begin{align}

\partial_\beta \equiv \frac{\partial}{\partial x^\beta}\ne g_{\mu\beta}\frac{\partial}{\partial x_\mu}=g_{\mu\beta}\partial^\alpha \end{align}$$

However, we can show that for the time component


 * $$\begin{align}

\alpha=0=\tau \end{align}$$

due to the fact that $$x^0=x_0$$, one has


 * $$\begin{align}

&\partial_0=\frac{\partial}{\partial x^0}=\frac{\partial}{\partial x_0}=g_{0\mu}\partial^\mu \\ &\partial^0=g^{0\mu}\partial_\mu \\ \end{align}$$

where the last step we use the explicit form of the metric.

For the spatial component, since $$g_{\mu\alpha,i}=0$$, one has


 * $$\begin{align}

&\partial_i \equiv \frac{\partial}{\partial x^i}=\left(\frac{\partial(g_{\mu\alpha} x^\alpha)}{\partial x^i}\right)\frac{\partial}{\partial x_\mu} =\left(g_{\mu i}+g_{\mu\alpha,i} x^\alpha \right)\frac{\partial}{\partial x_\mu}=g_{\mu i}\frac{\partial}{\partial x_\mu}=g_{i\mu}\partial^\mu\\ &\partial^i=g^{i\mu}\partial_\mu \end{align}$$

Putting the two pieces together, one obtains a rather compact result which is only valid in the case of hyperbolic coordinates


 * $$\begin{align}

&\partial_\beta\equiv \frac{\partial}{\partial x^\beta}=g_{\mu\beta}\frac{\partial}{\partial x_\mu}=g_{\mu\beta}\partial^\alpha \\ &\partial^\beta= g^{\mu\alpha}\partial_\alpha \end{align}$$

Similarly, one can show that


 * $$\begin{align}

&T^{\mu\nu}\partial^\alpha(g_{\beta\nu}g_{\alpha\mu})\\ &=T^{\mu\nu}\partial^0(g_{\beta\nu}g_{0\mu})\\ &=T^{\mu\nu}\partial_0(g_{\beta\nu}g_{0\mu})\\ &=T^{\mu\nu}\delta_{\mu 0}\partial_0(g_{\beta\nu})\\ &=T^{\mu\nu}\delta_{\mu 0}\delta_{\beta 3}\delta_{\nu 3}\partial_0(-\tau^2)\\ &=T^{03}\delta_{\beta 3}\partial_0(-\tau^2)\\ &=T^{03}\delta_{\beta 3}(-2\tau)\\ &=\frac{2}{\tau}T^{03}g_{\beta 3} \end{align}$$

Here we summarize all the relations we derived above


 * $$\begin{align}

&D^\mu=g^{\mu\nu}D_{\nu}=\Delta^{\mu\nu}D_\nu+u^\mu\frac{D}{D\tau}\equiv D^\mu_\perp+u^\mu\frac{D}{D\tau}\\ &\Delta_{\nu\beta}D^\alpha u^\beta=D^\alpha u_\nu\\ &\Delta_{\nu\beta}D^\alpha u^\beta=-u^\beta D^\alpha\Delta_{\nu\beta}=D^\alpha u_\nu\\ &\Gamma^\lambda_{\mu\nu}u_\lambda =\tau\delta_{\mu 3}\delta_{\nu 3}\gamma+\frac{u_3}{\tau}(\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0})\\ &u^{\mu}\Gamma^{\lambda}_{\mu\nu}u_{\lambda} = -{\tau}(u^3)^2\delta_{\nu 0}\\ &\partial_\beta\equiv \frac{\partial}{\partial x^\beta}=g_{\mu\beta}\partial^\alpha \\ &\partial^\beta= g^{\mu\alpha}\partial_\alpha \\ &u\cdot \partial u_i = u_{\mu}{(u_i)_;}^{\mu} = u_{\mu}{(u_i)_,}{^{\mu}} \end{align}$$

and the symbols we introduced in the section


 * $$\begin{align}

&\theta\equiv Du\equiv\partial \cdot u = \left(u^{\mu}\right)_{;\mu} = \frac{1}{\tau} \partial_{\mu} \left(\tau u^{\mu}\right) =\frac{\gamma}{\tau}+\partial_{\mu}u^{\mu} =\frac{\gamma}{\tau} + d_{\tau}\gamma+\gamma \partial_i v^{i} \\ &d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &\frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &d_{\tau}\gamma = \frac{d}{d\tau}\left(1-g_{ij}u^iu^j\right)^{1/2} = -\frac{g_{ij}u^i}{\gamma}\frac{du^j}{d\tau}-\frac{u^iu^j}{2\gamma}\frac{dg_{ij}}{d\tau} = -\frac{u_j}{\gamma}\frac{du^j}{d\tau}+\frac{\left(u^3\right)^2\tau}{\gamma}\\ &\partial_k\gamma=-v^i\partial_ku_i\\ &u\cdot \partial u_{i} = \gamma \frac{du_i}{d\tau} - u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} = \gamma \frac{du_i}{d\tau} = - \gamma\frac{du^1}{d\tau}\delta_{i1}- \gamma\frac{du^2}{d\tau}\delta_{i2}- \left(\gamma\tau^2\frac{du^3}{d\tau}+2\tau\gamma u^3 \right)\delta_{i3} \\ &\gamma\frac{du_i}{d\tau}=\gamma\frac{du_1}{d\tau}\delta_{i1}+\gamma\frac{du_2}{d\tau}\delta_{i2}+\gamma\frac{du_3}{d\tau}\delta_{i3}\\ \end{align}$$

The number of degrees of freedom and equations
The total number of variables is 10, namely, 1 entropy density, 3 velocity, 1 bulk viscosity and 5 components of shear viscosity. Correspondingly, one has 10 EoS, namely, 1 for the variation of entropy density, 3 for the 3-momentum, 1 for bulk viscosity and 5 for the shear viscosity. Though one has one equation for each conserved charges, in terms of SPH algorithm, conserved charges are attached to the SPH particles, therefore their conservation is naturally implied. We will start our discussion from the equations for viscosity, then the conservation of energy momentum tensor. All the equations and variables will be listed out explicitly at the end of this notes. As pointed out before, all the spatial derivatives involved will be evaluated in terms of SPH degree of freedom, while the derivative with respect to the time-like coordinate $$\tau$$ gives rise to the time evolution.

Shear viscosity only contains 5 independent variables. Before writing down their equations, we explicitly pick out five independent components as final variables and express others in terms of them. The five shear coefficients chosen in this note are $$\pi^{ij}$$ with $$i,j=1,2,3$$ except $$\pi^{33}$$. By making use of the properties that $$\pi^{ij}$$ are symmetric and perpendicular to the 4-velocity, one has


 * $$\begin{align}

&0=u_\mu\pi^{\mu\nu}=\gamma\pi^{0\nu}+u_i\pi^{i\nu}\\ &\pi^{0\nu}=-\frac{u_j}{\gamma}\pi^{j\nu}\\ &\pi^{0i}=-\frac{u_j}{\gamma}\pi^{ji}\\ &\pi^{00}=-\frac{u_j}{\gamma}\pi^{j0}=-\frac{u_j}{\gamma}\pi^{0j}=\frac{u_iu_j}{\gamma^2}\pi^{ij} \end{align}$$

Therefore any $$\pi^{ij}$$ with at least one time-like index can be expressed in terms of those only with spatial indexes. Noting that $$\pi^{33}$$ can be written in terms of others by the traceless condition, hence


 * $$\begin{align}

&0=g_{\mu\nu}\pi^{\mu\nu}=\pi^{00}-\pi^{11}-\pi^{22}-\tau^2\pi^{33}\\ &\pi^{00}=\pi^{11}+\pi^{22}+\tau^2\pi^{33}=\frac{u_iu_j}{\gamma^2}\pi^{ij} \end{align}$$

Therefore


 * $$\begin{align}

&\left(1-\frac{u_1^2}{\gamma^2}\right)\pi^{11}+\left(1-\frac{u_2^2}{\gamma^2}\right)\pi^{22}+\left(1-\frac{u_3^2}{\gamma^2}\right)\tau^2\pi^{33}=\sum_{i\ne j}\frac{u_iu_j}{\gamma^2}\pi^{ij}\\ &\pi^{33}=\left(1-\frac{u_3^2}{\gamma^2}\right)^{-1}\frac{1}{\tau^2}\left[\sum_{i\ne j}\frac{u_iu_j}{\gamma^2}\pi^{ij}-\left(1-\frac{u_1^2}{\gamma^2}\right)\pi^{11}-\left(1-\frac{u_2^2}{\gamma^2}\right)\pi^{22} \right] \end{align}$$

The equation of bulk viscosity
In our calculation, the variable for bulk viscosity is chosen to be $$\frac{\Pi}{\sigma}$$. This is convenient, since it involves the ratio of two densities, volume does not enter the definition explicitly. By making use of the fact that $$\sigma$$ is the density of a conserved charge, one can easily obtain the total bulk viscosity. According to the model of "memory effect" proposed by Koide et. al., the bulk viscosity satisfies the following equation


 * $$\begin{align}

\widetilde{\Pi}=-\int_{\tau_0}^{\tau} d\tau'G(\tau,\tau')\frac{\zeta'\theta'}{\sigma'}+e^{-(\tau-\tau_0)/\tau_R}\widetilde{\Pi}_0 \end{align}$$

The Green function satisfies


 * $$\begin{align}

&G(\tau,\tau;\tau_R) \equiv \frac{1}{\tau_R}\\ &\frac{DG(\tau,\tau';\tau_R)}{D\tau}=-\frac{G(\tau,\tau';\tau_R)}{\tau_R} \end{align}$$

As a result,


 * $$\begin{align}

G(\tau,\tau';\tau_R) \equiv \frac{1}{\tau_R}e^{-(\tau-\tau')/\tau_R} \end{align}$$

In flat space, one finally obtains


 * $$\begin{align}

\dot{\widetilde{\Pi}}=-\frac{\widetilde{\Pi}}{\tau_R}-\frac{\zeta\theta}{\sigma\tau_R} \end{align}$$

In hyperbolic coordinates, one has to generate the above equation, the generation in fact contains some freedom. According to Gabriel's choice, the time derivative is taken as the covariant time derivative (for instance, in the place of derivative of time-like coordinate which gives rise to an additional factor of $$\gamma$$)


 * $$\begin{align}

\dot{\widetilde{\Pi}}\equiv\frac{D\widetilde{\Pi}}{D\tau} \end{align}$$

one also substitute the expansion rate by the following covariant divergence


 * $$\begin{align}

\theta\equiv\partial\cdot u \end{align}$$

On on hand, one has


 * $$\begin{align}

\frac{D\widetilde{\Pi}}{D\tau}=u^\mu D_\mu\widetilde{\Pi}=u^\mu \partial_\mu\widetilde{\Pi}=\gamma d_\tau\left(\frac{\Pi}{\sigma} \right) \end{align}$$

On the other hand,


 * $$\begin{align}

d_{\tau}\left(\frac{\Pi}{\sigma}\right) = \frac{d_{\tau}\Pi}{\sigma}-\frac{\Pi}{\sigma^2}d_{\tau}\sigma = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma^2}\left(\frac{\sigma}{\gamma}\partial \cdot u \right) = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma}\left(\frac{1}{\gamma}\partial \cdot u \right) \end{align}$$

where


 * $$\begin{align}

0 = \left(\sigma u^{\mu}\right)_{;\mu} = u \cdot \partial \sigma + \sigma \partial \cdot u \end{align}$$

Therefore


 * $$\begin{align}

d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \end{align}$$

Combining the two expressions, one has


 * $$\begin{align}

d_{\tau}\left(\frac{\Pi}{\sigma}\right) = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma^2}\left(\frac{\sigma}{\gamma}\partial \cdot u\right) =\frac{1}{\gamma}\left(-\frac{\widetilde{\Pi}}{\tau_R}-\frac{\zeta\theta}{\sigma\tau_R} \right) =-\frac{1}{\gamma\tau_R\sigma}\left(\Pi+{\zeta\partial\cdot u}\right) \end{align}$$

We note an useful relation


 * $$\begin{align}

{d_{\tau}\Pi}=-\frac{1}{\gamma\tau_R}\left(\Pi+{\left(\zeta+{\tau_R\Pi} \right)\partial\cdot u}\right) \end{align}$$

It is worth mention that we will not go any further once we write down some expression in terms of those symbols and basic expressions discussed in the last section. This is because all those symbols and expressions are already shown to be functions of the final variables, and this helps a lot to reduce the size of the formulae and to avoid unnecessary typos.

The equations of shear viscosity
The equations of shear viscosity have similar forms, which read


 * $$\begin{align}

\widetilde{\pi}^{\mu\nu}=\int^\tau_{\tau_0}d\tau'G(\tau,\tau')\frac{\eta\sigma^{\mu\nu}(\tau')}{\sigma}+e^{-(\tau-\tau_0)/\tau_{R_\pi}}\widetilde{\pi}^{\mu\nu}_0 \end{align}$$

where one makes use of the projection operator.


 * $$\begin{align}

&\sigma^{\mu\nu}\equiv \Delta^{\mu\nu\alpha\beta}D_\alpha u_\beta\\ &\Delta^{\mu\nu\alpha\beta} \equiv \frac{1}{2}\left[\Delta^{\mu\alpha}\Delta^{\nu\beta}+\Delta^{\mu\beta}\Delta^{\nu\alpha}-\frac{2}{3}\Delta^{\mu\nu}\Delta^{\alpha\beta}\right] \end{align}$$

In fact, both sides of the equation are affected by the projection operator, which is based on the Curie's principle. One obtains


 * $$\begin{align}

&\frac{D\widetilde{\pi}^{\mu\nu}}{D\tau}=-\frac{\widetilde{\pi}^{\mu\nu}}{\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\sigma\tau_R}\\ &\frac{D\widetilde{\pi}^{\mu\nu}}{D\tau}=-\frac{{\pi}^{\mu\nu}}{\sigma\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\sigma\tau_R} \end{align}$$

The l.h.s. of the equation can be written as


 * $$\begin{align}

&\frac{D\widetilde{\pi}^{\mu\nu}}{D\tau} =u_\alpha D^\alpha\widetilde{\pi}^{\mu\nu} =u_\alpha D^\alpha\left(\frac{\pi^{\mu\nu}}{\sigma}\right)\\ &=\frac{1}{\sigma}u_\alpha D^\alpha{\pi^{\mu\nu}}-\frac{\pi^{\mu\nu}}{\sigma^2}u_\alpha\partial^\alpha{\sigma}\\ &=\frac{1}{\sigma}u_\alpha D^\alpha{\pi^{\mu\nu}}+\frac{\pi^{\mu\nu}}{\sigma^2}\left({\sigma}\partial \cdot u \right)\\ &=\frac{1}{\sigma}u_\alpha D^\alpha{\pi^{\mu\nu}}+\frac{\pi^{\mu\nu}}{\sigma}\left(\partial \cdot u \right) \end{align}$$

therefore


 * $$\begin{align}

u_\alpha D^\alpha{\pi^{\mu\nu}}=-\frac{{\pi}^{\mu\nu}}{\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\tau_R}-{\pi^{\mu\nu}}\left(\partial \cdot u \right) \end{align}$$

The r.h.s. of the above equation is traceless and symmetric, one may use the projection operator on the l.h.s. of it in order to increase the precision when carrying out numerical calculations


 * $$\begin{align}

\Delta_{\mu\nu\alpha\beta}u_\lambda D^\lambda{\pi^{\alpha\beta}}=-\frac{{\pi}_{\mu\nu}}{\tau_R}+\frac{\eta\sigma_{\mu\nu}}{\tau_R}-{\pi_{\mu\nu}}\left(\partial \cdot u \right) \end{align}$$

In the following, we will evaluate each term individually, until everything is in terms of known quantities


 * $$\begin{align}

&\eta\sigma_{\mu\nu}=\eta\frac{1}{2}\left[\Delta_{\mu\alpha}\Delta_{\nu\beta}+\Delta_{\mu\beta}\Delta_{\nu\alpha}-\frac{2}{3}\Delta_{\mu\nu}\Delta_{\alpha\beta}\right]D^\alpha u^\beta\\ &=\frac{\eta}{2}\Delta_{\mu\alpha}\Delta_{\nu\beta}\left[D^\alpha u^\beta+D^\beta u^\alpha\right]-\frac{\eta}{3}\Delta_{\mu\nu}\Delta_{\alpha\beta}D^\alpha u^\beta\\ &=\frac{\eta}{2}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3}\Delta_{\mu\nu}(\partial\cdot u)-\eta\gamma\tau\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0}) \end{align}$$

where


 * $$\begin{align}

u_{\mu}u^{\alpha}\Gamma^{\lambda}_{\alpha\nu}u_{\lambda} =-u_\mu\frac{1}{\tau^3}(u_3)^2\delta_{\nu 0} = -u_\mu{\tau}(u^3)^2\delta_{\nu 0} \end{align}$$

The above expression derivative of 4-velocity, but since there is no summation, it is not too complicated, one only need to evaluate the spatial derivation of the 4-velocity.

On the other hand, starting from the l.h.s. of the equation of shear viscosity, one has


 * $$\begin{align}

\tau_R\Delta_{\mu\nu\alpha\beta}\frac{D\pi^{\alpha\beta}}{D\tau} =\tau_R\frac{D\pi_{\mu\nu}}{D\tau}-\tau_R\pi^{\alpha\beta}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau} \end{align}$$

where


 * $$\begin{align}

&-\tau_R\pi^{\alpha\beta}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}\\ &=-\frac{1}{2}\tau_R\pi^{\alpha\beta}\frac{D}{D\tau}\left[\Delta^{\mu\alpha}\Delta^{\nu\beta}+\Delta^{\mu\beta}\Delta^{\nu\alpha}-\frac{2}{3}\Delta^{\mu\nu}\Delta^{\alpha\beta}\right]\\ &=\tau_R(u_\mu\pi^\alpha_\nu+u_\nu\pi^\alpha_\mu)\frac{Du_\alpha}{D\tau}\\ &=\tau_R\gamma(u_\mu\pi^\alpha_\nu+u_\nu\pi^\alpha_\mu)\frac{du_\alpha}{d\tau}+\tau_R(u_\mu\pi_{\nu 0}+u_\nu\pi_{\mu 0})\frac{(u_3)^2}{\tau^3}\\ &=\tau_R\gamma(u_\mu\pi^0_\nu+u_\nu\pi^0_\mu)\frac{d\gamma}{d\tau}+\tau_R\gamma(u_\mu\pi^j_\nu+u_\nu\pi^j_\mu)\frac{du_j}{d\tau}+\tau_R(u_\mu\pi_{\nu 0}+u_\nu\pi_{\mu 0})\frac{(u_3)^2}{\tau^3}\\ &=\tau_R(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} \end{align}$$

The last step has made use of the explicit form of $$\gamma$$ factor, which is not necessary according to our principle, but doing this may further simplify the expression


 * $$\begin{align}

\tau_R\frac{D\pi_{\mu\nu}}{D\tau} =\gamma\tau_R\frac{d\pi_{\mu\nu}}{d\tau}+\frac{\tau_R}{\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})+\frac{\tau_R}{\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3}) \end{align}$$

The sum of the two terms reads


 * $$\begin{align}

&\tau_R\Delta_{\mu\nu\alpha\beta}\frac{D\pi^{\alpha\beta}}{D\tau}\\ &=\tau_R\gamma(u_\mu\pi^0_\nu+u_\nu\pi^0_\mu)\frac{d\gamma}{d\tau}+\tau_R\gamma(u_\mu\pi^j_\nu+u_\nu\pi^j_\mu)\frac{du_j}{d\tau}+\tau_R(u_\mu\pi_{\nu 0}+u_\nu\pi_{\mu 0})\frac{(u_3)^2}{\tau^3}\\ &+\gamma\tau_R\frac{d\pi_{\mu\nu}}{d\tau}+\frac{\tau_R}{\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3}) +\frac{\tau_R}{\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3}) \end{align}$$

Putting all pieces together, one obtains the equations for shear viscosity


 * $$\begin{align}

&\tau_R(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} +\gamma\tau_R\frac{d\pi_{\mu\nu}}{d\tau}+\frac{\tau_R}{\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})\\ &+\frac{\tau_R}{\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3})+\pi_{\mu\nu}\\ &=+\frac{\eta}{2}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3}\Delta_{\mu\nu}(\partial\cdot u)\\ &-\eta\gamma\tau\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0})-\tau_R\pi_{\mu\nu}(\partial\cdot u) \end{align}$$

If we instead choose $$\begin{align} \frac{\pi_{\mu\nu}}{\sigma} \end{align}$$ as the variable for shear viscosity, one notices that l.h.s. of the equation with projection operator is


 * $$\begin{align}

\tau_R\Delta_{\mu\nu\alpha\beta}\frac{D\widetilde{\pi}^{\alpha\beta}}{D\tau} =\tau_R\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau}-\frac{\tau_R\pi^{\alpha\beta}}{\sigma}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}\\ \end{align}$$

therefore


 * $$\begin{align}

&\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau}=\frac{\pi^{\alpha\beta}}{\sigma}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}-\frac{{\pi}^{\mu\nu}}{\sigma\tau_R}+\frac{\eta\sigma^{\mu\nu}}{\sigma\tau_R}\\ &\tau_R\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau}=\frac{\tau_R\pi^{\alpha\beta}}{\sigma}\frac{D\Delta_{\mu\nu\alpha\beta}}{D\tau}-\frac{{\pi}^{\mu\nu}}{\sigma}+\frac{\eta\sigma^{\mu\nu}}{\sigma} \end{align}$$

Comparing with previous equation, there is an additional term on the r.h.s. By making use of the previous results, the l.h.s. of the equation can be written in terms of "normal" derivatives, one has


 * $$\begin{align}

\tau_R\frac{D\widetilde{\pi}_{\mu\nu}}{D\tau} =\gamma\tau_R\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}+\frac{\tau_R}{\sigma\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})+\frac{\tau_R}{\sigma\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\sigma\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3}) \end{align}$$

The corresponding equations read


 * $$\begin{align}

&\frac{\tau_R}{\sigma}(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} +\gamma\tau_R\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}+\frac{\tau_R}{\sigma\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})\\ &+\frac{\tau_R}{\sigma\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\sigma\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3})+\frac{\pi_{\mu\nu}}{\sigma}\\ &=\frac{\eta}{2\sigma}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2\sigma}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3\sigma}\Delta_{\mu\nu}(\partial\cdot u)\\ &-\frac{\eta\gamma\tau}{\sigma}\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\sigma\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\sigma\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0}) \end{align}$$

For future use, we also write down the following relations


 * $$\begin{align}

&\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}=\frac{d\left(\frac{\pi_{\mu\nu}}{\sigma} \right)}{d\tau}=\frac{1}{\sigma}\frac{d{\pi}_{\mu\nu}}{d\tau}-\frac{{\pi}_{\mu\nu}}{\sigma^2}d_\tau\sigma=\frac{1}{\sigma}\frac{d{\pi}_{\mu\nu}}{d\tau}+\frac{{\pi}_{\mu\nu}}{\sigma\gamma}(\partial\cdot u)\\ &\frac{d{\pi}_{\mu\nu}}{d\tau}=\sigma\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}-\frac{{\pi}_{\mu\nu}}{\gamma}(\partial\cdot u) \end{align}$$

Conservation of energy-momentum tensor, spatial part
Since these three equations contain the time derivative of velocity, which is acceleration, the equations can be rearranged to read off the information on how acceleration is determined by the "force", namely, the second law of Newton. The energy-momentum conservation can be put into the following form


 * $$\begin{align}

\frac{1}{\tau}\partial^\alpha(\tau T_{\alpha\beta})+g_{\beta 0}\tau T^{33}=0 \end{align}$$

The zeroth component of the energy momentum conservation can be seen as a combination of parts parallel and perpendicular to the four velocity, thus can be used to check the numerical precision.

For the time being, we will only focus on the spatial part of the equation. We first write down the energy-momentum tensor in the following form


 * $$\begin{align}

T^{\mu\nu}=(\varepsilon+P+\Pi)u^\mu u^\nu-(P+\Pi)g^{\mu\nu}+\pi^{\mu\nu} \end{align}$$

From the covariant form of conservation, the spatial part reads


 * $$\begin{align}

\sigma\gamma\frac{d}{d\tau}\left(\frac{\varepsilon+P+\Pi}{\sigma}u_i\right)+\frac{1}{\tau}\partial^\mu(\tau\pi_{\mu i})=\partial_i(P+\Pi) \end{align}$$

The first term on the l.h.s. can be written as


 * $$\begin{align}

&\sigma\gamma\frac{d}{d\tau}\left(\frac{\varepsilon+P+\Pi}{\sigma}u_i\right)\\ &=\gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}-u_i(\varepsilon+P)\frac{\gamma}{\sigma}\frac{d\sigma}{d\tau}+u_i{\gamma}\frac{d(\varepsilon+P)}{d\tau}+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}\\ &=\gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}+u_i(\varepsilon+P)(\partial\cdot u)+u_i{\gamma}\frac{d(\varepsilon+P)}{d\tau}+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}\\ \end{align}$$

On the r.h.s. of the above equation, only the third term nees some further attention, the part contain the derivative is


 * $$\begin{align}

\frac{d(\varepsilon+P)}{d\tau} =\frac{d(\varepsilon+P)}{d\varepsilon}\frac{d\varepsilon}{d\tau} =\left(1+\frac{dP}{d\varepsilon}\right)\frac{d\varepsilon}{d\tau} =\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}u^\mu D_\mu\varepsilon =\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}\left({T\sigma\gamma}d_\tau\left(\frac{s}{\sigma}\right)-(\varepsilon+P)(\partial\cdot u)\right) \end{align}$$

The last step made use of the thermodynamic relation


 * $$\begin{align}

&u \cdot \partial {\varepsilon} = u \cdot \partial {\varepsilon}+\omega \partial \cdot u - \omega \partial \cdot u \\ &= T u \cdot \partial s + T s \partial \cdot u + \sum_p \mu_p u \cdot \partial n_p + \sum_p \mu_p n_p \partial \cdot u - \omega \partial\cdot u \\ &= T \partial \cdot (su) + \sum_p \mu_p \partial \cdot (n_pu) - \omega \partial \cdot u \\ &= T \partial \cdot (su) - \omega \partial \cdot u  \\ &={T \sigma\gamma}d_{\tau}\left(\frac{s}{\sigma}\right) - \omega \partial \cdot u \end{align}$$

where one has made use of the properties of conserved current


 * $$\begin{align}

\sum_p \mu_p \partial \cdot (n_p u)=0 \end{align}$$

and the first law of thermodynamics


 * $$\begin{align}

&\partial \varepsilon=\partial(\frac{E}{V})=\frac{V\partial E-E\partial V}{V^2} \\ &\partial E =T\partial S -P \partial V +\sum_p \mu \partial N_p\\ &\partial S = \partial (sV)=V\partial s+s\partial V\\ &\sum_p \mu \partial N_p=\sum_p \mu V\partial n_p+\sum_p \mu n_p \partial V\\ &\partial \varepsilon =T\partial s+\sum_p\mu_p \partial n_p \end{align}$$

where the terms involves the pressure cancel each other naturally. Finally, make use of the properties of extensive quantity


 * $$\begin{align}

\omega=\varepsilon+p=Ts-p+\sum_p \mu \partial n_p+p=Ts+\sum_p \mu \partial n_p \end{align}$$

The time derivative in the second term on the l.h.s. can be written in terms of time-like coordinate, the apacial derivative may remain unchanged since they can be evaluated in terms of SPH degree of freedom.


 * $$\begin{align}

&\frac{1}{\tau}\partial^\mu(\tau\pi_{\mu i})=\partial^0\pi_{0i}+\partial^j\pi_{ji}+\frac{1}{\tau}\pi_{\mu i}\partial^\mu\tau=\frac{d\pi_{0i}}{d\tau}-v^j\partial_j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}\\ &=\frac{d\pi_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau} =\sigma\frac{d\widetilde{\pi}_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}-\frac{{\pi}_{0i}}{\gamma}(\partial\cdot u)\\ \end{align}$$

Putting all pieces together, one has


 * $$\begin{align}

&\partial_i(P+\Pi)= \gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}+u_i(\varepsilon+P)(\partial\cdot u)+u_i{\gamma}\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}\left({T\sigma\gamma}d_\tau\left(\frac{s}{\sigma}\right)-(\varepsilon+P)(\partial\cdot u)\right)\\ &+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}+\sigma\frac{d\widetilde{\pi}_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}-\frac{{\pi}_{0i}}{\gamma}(\partial\cdot u) \end{align}$$

Conservation of the energy-momentum tensor, entropy flow
Instead of taking the zeroth component of energy-momentum conservation, we project the equation onto the direction of 4-velocity. It is well known that, the resulting equation describes the time variation of entropy current, which is not a conserved flow in the presence of viscosity. The projection gives


 * $$\begin{align}

&0=u_\mu D_\nu T^{\mu\nu}\\ &=u_\mu D_\nu((\varepsilon+P+\Pi)u^\mu u^\nu-(P+\Pi)g^{\mu\nu}+\pi^{\mu\nu})\\ &=\frac{D\varepsilon}{D\tau}+(\varepsilon+P+\Pi)\theta+u_\mu D_\nu\pi^{\mu\nu}\\ &=\frac{D\varepsilon}{D\tau}+(\varepsilon+P+\Pi)\theta-\pi^{\mu\nu}D_\nu u_\mu\\ &=\frac{D\varepsilon}{D\tau}+(\varepsilon+P+\Pi)\theta-\pi^{\mu\nu}\sigma_{\mu\nu}\\ \end{align}$$

where one involves to separate any vector into the sum of parallel and perpendicular parts


 * $$\begin{align}

&D_{\mu}u_{\nu}=u_{\mu}\frac{Du_{\nu}}{D\tau}+D_{\mu\perp}u_{\nu} \\ &=u_{\mu}\frac{Du_{\nu}}{D\tau}+\frac{1}{2}(D_{\mu\perp}u_{\nu}+D_{\nu\perp}u_{\mu})+\frac{1}{2}(D_{\mu\perp}u_{\nu}-D_{\nu\perp}u_{\mu})\\ &=u_{\mu}\frac{Du_{\nu}}{D\tau}+\sigma_{\mu\nu}+\frac{1}{3}\Delta_{\mu\nu}\theta+\Omega_{\mu\nu}\\ &\sigma_{\mu\nu}\equiv\frac{1}{2}(\nabla_{\mu}u_{\nu}+\nabla_{\nu}u_{\mu}-\frac{2}{3}\Delta_{\mu\nu}\theta)\\ &\theta\equiv D_{\mu}u^{\mu}\\ &\Omega_{\mu\nu}\equiv\frac{1}{2}(D_{\mu\perp}u_{\nu}-D_{\nu\perp}u_{\mu})\\ &A^{<\mu\nu>}\equiv\Delta^{\mu\nu\alpha\beta}A_{\alpha\beta}\\ &\nabla^{<\mu}u^{\nu>}\equiv 2\sigma^{\mu\nu}=\Delta^{\mu\nu\alpha\beta}D_\alpha u_\beta=2\nabla^{(\mu}u^{\nu)}-\frac{2}{3}\Delta^{\mu\nu}\theta\\ &\nabla^{(\mu}u^{\nu)}\equiv D_{\mu\perp}u^{\nu}+D_{\nu\perp}u^{\mu} \end{align}$$

By making use of the relations of thermodynamics and conserved current


 * $$\begin{align}

&d\varepsilon=Tds+\mu_P dn_P\\ &\frac{D\varepsilon}{D\tau}=T\frac{Ds}{D\tau}+\mu_P\frac{Dn_P}{D\tau}\\ &\frac{D\varepsilon}{D\tau}=u^\mu D_\mu\epsilon=Tu^\mu D_\mu s+\mu_P u^\mu D_\mu n_P =TD_\mu(su^\mu)+\mu_PD_\mu(n_Pu^\mu)-(Ts+\mu_Pn_P)\theta\\ &\varepsilon=-P+Ts+\mu_Pn_P \end{align}$$

one obtains


 * $$\begin{align}

TD_\mu(su^\mu)=\Pi(-\theta)+\pi^{\mu\nu}\sigma_{\mu\nu}=\Pi(-(\partial\cdot u))+\pi^{\mu\nu}\sigma_{\mu\nu} \end{align}$$

Now, we want to express anything in terms of the time derivative of final variables. The complication comes from the last term in the summation, we have to write down explicitly the summation and treat the time derivative differently from the spatial derivative


 * $$\begin{align}

&\pi^{\mu\nu}\sigma_{\mu\nu}=\pi^{\mu\nu}D_\mu u_\nu=\pi^{\mu\nu}\partial_\mu u_\nu-\tau\gamma\pi^{33}-\frac{2u_3}{\tau}\pi^{03}\\ &=\left(\pi^{0j}-\pi^{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}-\frac{(u_3)^2}{\gamma\tau^3}\pi^{00}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\tau\gamma\pi^{33}-\frac{2u_3}{\tau}\pi^{03}\\ &=\left(\pi_0^j-\pi_{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\frac{\gamma}{\tau^3}\pi_{33}+\frac{2u_3}{\tau^3}\pi_{03}-\frac{(u_3)^2}{\gamma\tau^3}\pi_{00} \end{align}$$

By taking into consideration


 * $$\begin{align}

T D_\mu(su^\mu)=T\sigma\gamma d_\tau(\frac{s}{\sigma})=\Pi(-\theta)+\pi^{\mu\nu}\sigma_{\mu\nu}=\Pi(-(\partial\cdot u))+\pi^{\mu\nu}\sigma_{\mu\nu} =\Pi(-(\partial\cdot u))+\pi^{\mu\nu}\sigma_{\mu\nu} \end{align}$$

For future use, we note that


 * $$\begin{align}

\frac{ds}{d\tau}=\sigma{d_\tau(\frac{s}{\sigma})}-\frac{s}{\gamma}(\partial\cdot u) \end{align}$$

Finally we obtain the equation for entropy flow


 * $$\begin{align}

\\d_\tau(\frac{s}{\sigma}) =\frac{1}{T\sigma\gamma}\Pi(-(\partial\cdot u))+\frac{1}{T\sigma\gamma} \left[\left(\pi_0^j-\pi_{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\frac{\gamma}{\tau^3}\pi_{33}+\frac{2u_3}{\tau^3}\pi_{03}-\frac{(u_3)^2}{\gamma\tau^3}\pi_{00}\right] \end{align}$$

The list of equations and variables
There are in total 10 variables, namely,


 * $$\begin{align}

u^i,\frac{s}{\sigma},\frac{\Pi}{\sigma},\frac{\pi^{\mu\nu}}{\sigma} \end{align}$$

among them, 5 variables of shear viscosity are chosen as those with spatial indexes except $$\pi^{33}$$. All other shear coefficients can be written in terms of these 5 variables by the following relations


 * $$\begin{align}

&\pi^{33}=\left(1-\frac{u_3^2}{\gamma^2}\right)^{-1}\frac{1}{\tau^2}\left[\sum_{i\ne j}\frac{u_iu_j}{\gamma^2}\pi^{ij}-\left(1-\frac{u_1^2}{\gamma^2}\right)\pi^{11}-\left(1-\frac{u_2^2}{\gamma^2}\right)\pi^{22} \right]\\ &\pi^{0i}=-\frac{u_j}{\gamma}\pi^{ji}\\ &\pi^{00}=\pi^{11}+\pi^{22}+\tau^2\pi^{33}=\frac{u_iu_j}{\gamma^2}\pi^{ij} \end{align}$$

The equations of motion are the following


 * $$\begin{align}

d_{\tau}\left(\frac{\Pi}{\sigma}\right) = \frac{d_{\tau}\Pi}{\sigma} + \frac{\Pi}{\sigma^2}\left(\frac{\sigma}{\gamma}\partial \cdot u\right) =\frac{1}{\gamma}\left(-\frac{\widetilde{\Pi}}{\tau_R}-\frac{\zeta\theta}{\sigma\tau_R} \right) =-\frac{1}{\gamma\tau_R\sigma}\left(\Pi+{\zeta\partial\cdot u}\right) \end{align}$$
 * $$\begin{align}

&\frac{\tau_R}{\sigma}(\gamma u_\mu\pi^j_\nu+\gamma u_\nu\pi^j_\mu-u_\mu\pi^0_\nu u^j-u_\nu\pi^0_\mu u^j)\frac{du_j}{d\tau} +\gamma\tau_R\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}+\frac{\tau_R}{\sigma\tau^3}u_3(\delta_{\mu 0}\pi_{\nu 3}+\delta_{\nu 0}\pi_{\mu 3})\\ &+\frac{\tau_R}{\sigma\tau}u_3(\delta_{\mu 3}\pi_{\nu 0}+\delta_{\nu 3}\pi_{\mu 0}) -\frac{\gamma\tau_R}{\sigma\tau}(\delta_{\mu 3}\pi_{\nu 3}+\delta_{\nu 3}\pi_{\mu 3})+\frac{\pi_{\mu\nu}}{\sigma}\\ &=\frac{\eta}{2\sigma}[\partial_\mu u_\nu+\partial_\nu u_\mu]-\frac{\eta\gamma}{2\sigma}\left[u_\mu\frac{du_\nu}{d\tau}+u_\nu\frac{du_\mu}{d\tau}\right]-\frac{\eta}{3\sigma}\Delta_{\mu\nu}(\partial\cdot u)\\ &-\frac{\eta\gamma\tau}{\sigma}\delta_{\nu 3}\delta_{\mu 3}-\eta\frac{u_3}{\sigma\tau}[\delta_{\mu 0}\delta_{\nu 3}+\delta_{\mu 3}\delta_{\nu 0}]-\frac{\eta(u_3)^2}{2\sigma\tau^3}(u_\mu\delta_{\nu 0}+u_\nu\delta_{\mu 0}) \end{align}$$
 * $$\begin{align}

&\partial_i(P+\Pi)= \gamma(\varepsilon+P+\Pi)\frac{du_i}{d\tau}+u_i(\varepsilon+P)(\partial\cdot u)+u_i{\gamma}\left(1+\frac{dP}{d\varepsilon}\right)\frac{1}{\gamma}\left({T\sigma\gamma}d_\tau\left(\frac{s}{\sigma}\right)-(\varepsilon+P)(\partial\cdot u)\right)\\ &+u_i\gamma\frac{d\widetilde{\Pi}}{d\tau}+\sigma\frac{d\widetilde{\pi}_{0i}}{d\tau}-v_j\partial^j\pi_{0i}+\partial^j\pi_{ji}+\frac{\pi_{0i}}{\tau}-\frac{{\pi}_{0i}}{\gamma}(\partial\cdot u) \end{align}$$
 * $$\begin{align}

\\d_\tau(\frac{s}{\sigma}) =\frac{1}{T\sigma\gamma}\Pi(-(\partial\cdot u))+\frac{1}{T\sigma\gamma} \left[\left(\pi_0^j-\pi_{00}\frac{u^j}{\gamma}\right)\frac{du_j}{d\tau}+(\pi^{ij}+\pi^{00}v^iv^j-\pi^{i0}v^j-\pi^{0j}v^i)\partial_iu_j-\frac{\gamma}{\tau^3}\pi_{33}+\frac{2u_3}{\tau^3}\pi_{03}-\frac{(u_3)^2}{\gamma\tau^3}\pi_{00}\right] \end{align}$$

where one has used the following notations


 * $$\begin{align}

&\theta\equiv Du\equiv\partial \cdot u = \left(u^{\mu}\right)_{;\mu} = \frac{1}{\tau} \partial_{\mu} \left(\tau u^{\mu}\right) =\frac{\gamma}{\tau}+\partial_{\mu}u^{\mu} =\frac{\gamma}{\tau} + d_{\tau}\gamma+\gamma \partial_i v^{i} \\ &d_\tau\sigma=\frac{d\sigma}{d\tau}=-\frac{\sigma}{\gamma}(\partial\cdot u) \frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &\frac{du_\mu}{d\tau}=\frac{d\gamma}{d\tau}\delta_{\mu 0}+\frac{du_i}{d\tau}\delta_{\mu i}\\ &d_{\tau}\gamma = \frac{d}{d\tau}\left(1-g_{ij}u^iu^j\right)^{1/2} = -\frac{g_{ij}u^i}{\gamma}\frac{du^j}{d\tau}-\frac{u^iu^j}{2\gamma}\frac{dg_{ij}}{d\tau} = -\frac{u_j}{\gamma}\frac{du^j}{d\tau}+\frac{\left(u^3\right)^2\tau}{\gamma}\\ &\partial_k\gamma=-v^i\partial_ku_i\\ &u\cdot \partial u_{i} = \gamma \frac{du_i}{d\tau} - u^{\mu} u_{\lambda} \Gamma^{\lambda}_{i\mu} = \gamma \frac{du_i}{d\tau} = - \gamma\frac{du^1}{d\tau}\delta_{i1}- \gamma\frac{du^2}{d\tau}\delta_{i2}- \left(\gamma\tau^2\frac{du^3}{d\tau}+2\tau\gamma u^3 \right)\delta_{i3} \\ &\gamma\frac{du_i}{d\tau}=\gamma\frac{du_1}{d\tau}\delta_{i1}+\gamma\frac{du_2}{d\tau}\delta_{i2}+\gamma\frac{du_3}{d\tau}\delta_{i3}\\ \end{align}$$

In some cases, the following relations may be found useful,


 * $$\begin{align}

{d_{\tau}\Pi}=-\frac{1}{\gamma\tau_R}\left(\Pi+{\left(\zeta+{\tau_R\Pi} \right)\partial\cdot u}\right) \end{align}$$
 * $$\begin{align}

\frac{d{\pi}_{\mu\nu}}{d\tau}=\sigma\frac{d\widetilde{\pi}_{\mu\nu}}{d\tau}-\frac{{\pi}_{\mu\nu}}{\gamma}(\partial\cdot u) \end{align}$$
 * $$\begin{align}

\frac{ds}{d\tau}=\sigma{d_\tau(\frac{s}{\sigma})}-\frac{s}{\gamma}(\partial\cdot u) \end{align}$$