VALIDITY OF HYDROSTATIC EQUILIBRIUM IN GALAXY CLUSTERS FROM COSMOLOGICAL HYDRODYNAMICAL SIMULATIONS

Our analysis method to examine the validity of HSE using gas in a simulated cluster is based on the Euler equations:

$$\begin{equation} \frac{\partial \boldsymbol{v}}{\partial t}+(\boldsymbol{v}\cdot \nabla)\boldsymbol{v}=-\frac{1}{\rho _{\mathrm{gas}}}\nabla {p}-\nabla {\phi }, \end{equation} \tag{ 1 }$$

where ϕ is the gravitational potential, and ρ_gas, $\boldsymbol{v}$, and p are the density, velocity, and pressure of gas. While the Jeans equations do not describe the gas dynamics in their original form, Rasia et al. (2004) and Lau et al. (2009) added the gas pressure gradient term to the Jeans equations and adopted the resulting equations when analyzing simulated clusters. We argue in the next subsection that this is not justified, and present the result based on the Jeans equations but using collisionless dark matter particles in Appendix B, for reference.

We define the total mass M_tot of a cluster inside a volume V as

$$\begin{equation} M_{\rm tot}=\int _V d^3x\ \rho _{\mathrm{tot}}, \end{equation} \tag{ 2 }$$

where ρ_tot is the total density of the cluster. For the simulated cluster considered throughout this paper, ρ_tot consists of densities of gas, dark matter, and stars, i.e., ρ_tot = ρ_gas + ρ_dm + ρ_star. The total mass can be rewritten in terms of p and $\boldsymbol{v}$ using Poisson's equation and Gauss's theorem:

$$\begin{eqnarray} M_{\rm tot} &=& \frac{1}{4\pi G}\int _{\partial V} d\boldsymbol{S}\cdot \nabla \phi \nonumber\\ &=&\frac{1}{4\pi G}\int _{\partial V} d\boldsymbol{S} \cdot \left[-\frac{1}{\rho _{\mathrm{gas}}}\nabla p -(\boldsymbol{v}\cdot \nabla)\boldsymbol{v} -\frac{\partial \boldsymbol{v}}{\partial t}\right],\qquad \end{eqnarray} \tag{ 3 }$$

where ∂V is the surface surrounding the volume V and we have used Equation (1) in the second equality. Now the total mass is evaluated by the gas quantities alone, without any knowledge of dark matter and stars. This is why the present method is applicable, in principle, to the X-ray data of galaxy clusters.

If we adopt a spherical surface as ∂V, then the total mass can be decomposed into the following four effective mass terms:

$$\begin{equation} M_{\rm tot}=M_{\rm therm}+M_{\rm rot}+M_{\rm stream}+M_{\rm accel}, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} M_{\rm therm}=-\frac{1}{4\pi G}\int _{\partial V} dS\ \frac{1}{\rho _{\mathrm{gas}}}\frac{\partial p}{\partial r}, \end{equation} \tag{ 5 }$$

$$\begin{equation} M_{\rm rot}=\frac{1}{4\pi G}\int _{\partial V} dS\ \frac{v_\theta ^2+v_\varphi ^2}{r}, \end{equation} \tag{ 6 }$$

$$\begin{equation} M_{\rm stream}=-\frac{1}{4\pi G}\int _{\partial V} dS\ \left[v_r\frac{\partial v_r}{\partial r}+\frac{v_\theta }{r}\frac{\partial v_r}{\partial \theta }+\frac{v_\varphi }{r\sin \theta }\frac{\partial v_r}{\partial \varphi }\right], \end{equation} \tag{ 7 }$$

and

$$\begin{equation} M_{\rm accel}=-\frac{1}{4\pi G}\int _{\partial V} d S\ \frac{\partial v_r}{\partial t}. \end{equation} \tag{ 8 }$$

We emphasize here that the above set of equations does not assume spherical symmetry of the system; we just take a spherical surface as the integral surface and write down Equation (3) in spherical coordinates.

The first term, M_therm, originates from the thermal pressure gradient of gas. If the gas motion is negligible, then M_therm should be equal to the total mass M_tot, and thus regarded as a cluster mass estimated under the HSE assumption. In other words, the difference between M_tot and M_therm is a quantitative measure of the departure from the HSE assumption.

The inertial term $(\boldsymbol{v}\cdot \nabla)\boldsymbol{v}$ in the Euler equations reduces to M_rot and M_stream. If there exists a coherent rotational motion around the center of the cluster, then M_rot can be interpreted as the centrifugal force term. Without such a motion, however, the local tangential velocity of different directions at different locations on the sphere could make M_rot very large. During the course of cluster evolution, gas generally falls toward the center of the cluster with larger streaming speed and increasing radius from the center. In this case, M_stream becomes negative and cannot be neglected. On the other hand, it becomes positive and/or negligible, especially in the innermost region where the gas velocity is more randomized than that in the outer regions.

Finally, the acceleration term, M_accel, corresponds to the −∂v_r/∂t term in the Euler equations and becomes positive/negative when gas is decelerating/accelerating.

All the mass terms are invariant with respect to the choice of the axis of the spherical coordinates, but are not necessarily positive. Note also that M_rot and M_stream, corresponding to the inertial term, are not invariant with respect to the Galilean transformation. Thus, we evaluate those in the center-of-mass frame of the entire simulated cluster.

The set of basic equations that we adopt in this paper is essentially identical to that of Fang et al. (2009), except for the fact that they interpreted the difference between M_tot and M_therm + M_rot + M_stream as "turbulent gas motion" while we call it the acceleration term M_accel, as directly implied from Equation (1), and evaluate it from the residual, M_accel = M_tot − M_therm − M_rot − M_stream.

We are not sure why they ascribed the term to the turbulent motion. It is true that part of the gas acceleration would be due to the turbulent gas motion, but not entirely. Furthermore, the numerical simulation does not include any physical processes directly related to the turbulent motion. Even if the turbulent motion might be important for real clusters, it should come from some physics below the subgrid scales that cannot be properly resolved in the numerical simulation. The effects of gas random motion above the resolved scales should be included in M_rot and M_stream. We note, however, that the different interpretation of M_accel does not at all affect the conclusion of Fang et al. (2009) that gas rotation is the most important term to describe the origin of departure from HSE in their simulated clusters. As we will show below, this is not consistent with our result.

In previous literature, the Jeans equations are sometimes used to analyze the gas motion in simulated clusters. For that purpose, a thermal pressure gradient term is added by hand to the basic equations (e.g., Rasia et al. 2004; Lau et al. 2009). Assuming the steady state, i.e., $\partial \boldsymbol{v}/\partial t=0$, the Jeans equations in the r-direction (Equation (A20)) are now replaced by the following equation:

$$\begin{eqnarray} && \left[v_r\frac{\partial }{\partial r}+\frac{v_\theta }{r}\frac{\partial }{\partial \theta }+\frac{v_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]v_r \nonumber\\ && \quad +\frac{1}{\rho _{\mathrm{gas}}} \left[\frac{1}{r}\frac{\partial \left(\rho _{\mathrm{gas}}\sigma _{r\theta }^2 \right)}{\partial \theta } +\frac{1}{r\sin \theta }\frac{\partial \left(\rho _{\mathrm{gas}}\sigma _{r\varphi }^2 \right)}{\partial \varphi }\right] +\frac{\sigma _{r\theta }^2\cot \theta }{r} \nonumber\\ && \quad -\frac{v_\theta ^2+v_\varphi ^2}{r}=-\frac{1}{\rho _{\mathrm{gas}}}\frac{\partial \left(\rho _{\mathrm{gas}}\sigma _{rr}^2 \right)}{\partial r} -\frac{2\sigma _{rr}^2-\sigma _{\theta \theta }^2-\sigma _{\varphi \varphi }^2}{r} \nonumber\\ && \quad -\frac{1}{\rho _{\mathrm{gas}}}\frac{\partial p}{\partial r}-\frac{\partial \phi }{\partial r}, \end{eqnarray} \tag{ 9 }$$

where $\sigma _{ij}^2$ denotes the ij-component of the velocity dispersion tensor.

For example, Lau et al. (2009) converted each term of Equation (9) to the corresponding mass term; see their Equations (6)–(11). Then, they interpreted that the first two terms on the right-hand side of the above equation originate from "random gas motion." We suspect that the implicit assumption underlying Equation (9) is that the ICM consists of two distinct gas components: the thermal gas and the unthermalized gas, such as the cold gas accreting along with galaxies into clusters. While it may be reasonable to take into account both components separately when dealing with real galaxy clusters, it remains to be determined if the extra pressure gradient term can be added by hand to the conventional Jeans equations assuming that the two components share the same density ρ_gas. (In Appendix A, we show that diagonal components of the velocity dispersion tensor in the Jeans equations correspond to thermal pressure in the Euler equations.) It should also be noted that the current numerical simulations are performed entirely on the basis of the Euler equations. Therefore, for the sake of consistency, it is appropriate to use the Euler equations, instead of Equation (9), to evaluate the mass terms of the simulated galaxy clusters. The proper treatment of the unthermalized gas in numerical simulations is certainly an important issue, but is beyond the scope of the present paper.

We briefly describe the main features of the current simulation and refer readers to Cen (2012) for more detail. The simulation was run first with a low-resolution mode in a periodic box of 120 h⁻¹ Mpc on a side. Then, a region centered on a cluster with a mass of ∼2 × 10¹⁴ h⁻¹ M_☉ was resimulated with a higher resolution in an adaptively refined manner. The size of the refined region is 21 × 24 × 20 (h⁻¹ Mpc)³. The mean interparticle separation and the dark matter particle mass in the refined region are 117 h⁻¹ kpc comoving and 1.07 × 10⁸ h⁻¹ M_☉, respectively.

Star particles are created according to the prescription of Cen & Ostriker (1992). Their typical mass is ∼10⁶ M_☉. The simulation includes a metagalactic UV background (Haardt & Madau 1996), shielding of UV radiation by neutral hydrogen (Cen et al. 2005), and metallicity-dependent radiative cooling (Cen et al. 1995). While supernova feedback is modeled following Cen et al. (2005), active galactic nucleus feedback is not included in this simulation. The cosmological parameters used in this simulation are $(\Omega _b, \Omega _m, \Omega _\Lambda , h, n_s, \sigma _8) = (0.046, 0.28, 0.72, 0.70, 0.96, 0.82)$, following the WMAP7-normalized ΛCDM model (Komatsu et al. 2011).

Then, a cluster is identified and the cubic box of a side of 3.8 h⁻¹ Mpc surrounding the entire cluster is extracted from the simulation data. The dark matter and stars are represented by particles, and the temperature and density of gas are given on the 520³ grids (the grid length is 7.324 h⁻¹ kpc).

The radius r₅₀₀ of the cluster is ∼640 h⁻¹ kpc (r₅₀₀ is defined so that the mean density inside r₅₀₀ is 500 times the critical density of the universe). The center-of-mass velocity of the cluster within r₅₀₀ is set to vanish. The total mass M₅₀₀ within r₅₀₀ is ∼2 × 10¹⁴ M_☉. The average ICM temperature at r₅₀₀ is ∼2 keV, and the circular speed there is $v_{500}=\sqrt{GM_{500}/r_{500}}\sim 1000$ km s⁻¹.

Projected surface densities of gas, dark matter, and stars on the x–z plane are plotted in the left panels of Figure 1. The right panels of Figure 1 show the three-dimensional view of the three components. The left (right) plots are color-coded according to the surface (space) densities normalized by the fraction of each component averaged over the box, $\tilde{\Omega }_k$ (k =  gas, dark matter, and stars). Note that the fraction $\tilde{\Omega }_k$ is different from the density parameter Ω_k because the box is selected preferentially around the cluster.

Zoom In Zoom Out Reset image size

As is clear from Figure 1, the gas distribution is smoother but traces the underlying dark matter distribution very well. In contrast, stars are more significantly concentrated in high density regions, and exhibit numerous small clumps, most of which are not identified/resolved in the gas distribution.

Figure 2 plots the radial density and mass profiles of the cluster. The stellar fraction in the inner region (r < 200 h⁻¹ kpc) is significantly higher than the typical observed value. This is a well-known common problem among current high-resolution cosmological simulations, and implies that some important baryon physics including high-energy phenomena and star formation is still missing in the simulation. We perform the analysis of cluster gas, assuming that the excess star densities in the inner region do not affect our conclusions at the outer radius.

Zoom In Zoom Out Reset image size

Figure 3 represents velocity fields in the x–y, y–z, and z–x planes passing through the center of the cluster. The red/blue arrows have negative/positive radial velocity, showing that the gas in the outer regions (r ≳ 1 h⁻¹ Mpc) falls toward the center while its direction is randomized in the inner region.

Zoom In Zoom Out Reset image size

We evaluate the effective mass terms defined in Section 2 for the simulated cluster, which is plotted in Figure 5; the left panel shows the mass profiles, while the right panel indicates their fractional contribution to the total mass within the radius.

The total mass M_tot(r) is computed by directly summing all the dark matter and star particles and gas of grids within the sphere of r. The other terms, M_therm, M_rot, and M_stream, require the pressure and velocity fields evaluated at r. For that purpose, we use the density, velocity, and temperature of gas defined at 520³ original grid points, and first bin the cluster into 50 logarithmically equal radial intervals between 90 and 1900 h⁻¹ kpc, 10 × 10 linearly equal angular intervals.

Integrands of Equations (5)–(7) are calculated in each bin. Derivatives of the physical variables such as ∂p/∂r are calculated as follows; first, ∂/∂x, ∂/∂y, and ∂/∂z are computed from the difference of the adjacent and the original Cartesian grid points. Then, ∂/∂r, ∂/∂θ, and ∂/∂φ in our spherical coordinates are calculated by applying the chain rule.

In this way, M_therm(r), M_rot(r), and M_stream(r) are computed by integrating the corresponding integrands evaluated above. Finally, we estimate M_accel simply from the residual of M_accel = M_tot − M_therm − M_rot − M_stream, since we have the cluster data at z = 0 alone. This estimate for M_accel may be different from the original definition, i.e., Equation (8). Indeed, when we attempted to compute M_accel directly from the box mentioned in Section 3, it turned out to be too small to obtain a correct gravitational potential for the entire cluster. Thus, we go back to a larger simulation box with a side of 22.5 h⁻¹ Mpc and a grid length of 29.34 h⁻¹ kpc that encloses our cluster. Then, we compute the gravitational potential using a fast Fourier transform to obtain the gas acceleration at each grid point. This enables us to directly calculate M_accel. Figure 4 is a comparison of M_accels calculated by two methods. The directly calculated M_accel (magenta line) is in good agreement with M_tot − M_therm − M_rot − M_stream (black line), although there is a large difference between the two within 200 h⁻¹ kpc. Also, we make sure that the sum M_therm + M_rot + M_stream + M_accel reproduces M_tot within ∼2%, except for the innermost region (r < 200 h⁻¹ kpc), where it deviates from M_tot by up to ∼9%. Thus, the estimation of M_accel by M_tot − M_therm − M_rot − M_stream is sufficiently good given the quoted errors of our conclusion below. Although it seems better to use M_accel directly calculated from Equation (8), the grid size of the larger box is so coarse that we cannot take advantage of the high resolution of the simulations in this study. Therefore, we decided to use the smaller box explained in Section 3 and M_accel is calculated by M_tot − M_therm − M_rot − M_stream in the following analysis.

Zoom In Zoom Out Reset image size

The left panel of Figure 5 implies that M_therm agrees with M_tot reasonably well. Each of the other three terms contributes less than 10% of the total mass (the dotted curves correspond to the case in which each term becomes negative and its absolute value is plotted instead).

Zoom In Zoom Out Reset image size

In order to consider the validity of HSE more quantitatively, we plot the fractional contribution of each mass term in the right panel of Figure 5. The rotation term, M_rot, is always positive (by definition) and contributes approximately 10% almost independent of radius. In contrast, the streaming velocity term, M_stream, is mostly negative, and varies a lot at different radial bins. As a result, the difference of the total mass M_tot and the HSE mass M_therm is mostly explained by the acceleration term M_accel alone; compare the black and magenta curves in the right panel of Figure 5. At r = r₅₀₀ and r₂₀₀, the deviation from HSE in terms of the mass difference (M_therm − M_tot)/M_tot is about 10%. Nevertheless, the value significantly varies at different radii and it is safe to conclude that (M_therm − M_tot)/M_tot ranges approximately 10%–20% at r < r₂₀₀. Also, there is no systematic trend of the validity of the HSE assumption as a function of radius. Even though the reliability of the simulation is suspicious for r < 200 h⁻¹ kpc due to the excessive stellar concentration (Section 3.1), (M_therm − M_tot)/M_tot fluctuates between −10% and +25% for 300 h⁻¹ kpc < r < r₅₀₀. Thus, there is no guarantee that HSE becomes a better approximation toward the inner central region.

Physically speaking, M_rot + M_stream + M_accel(= M_tot − M_therm) corresponds to the term intergrating the Lagrangian derivative of the gas velocity over the sphere. Therefore, the fact that it is small compared with M_therm and M_tot is simply translated into the condition of HSE that the gas acceleration from a Lagrangian point of view is negligible compared with the pressure gradient and the total gravity.

Since we analyze a single simulated cluster, it is not clear to what extent our interpretation that M_accel determines the departure from HSE holds in general. Thus, we use a different set of smoothed particle hydrodynamic (SPH) simulation clusters (Dolag et al. 2009) kindly provided by Klaus Dolag. For these clusters, we assume spherical symmetry and do not divide the spherical surface at a given radius just for simplicity. Klaus Dolag also provides us with M_accel directly computed from the acceleration data. We find that in most cases, |M_accel| is larger than M_rot and |M_stream|, especially where M_therm deviates from M_tot. We also confirm that M_therm + M_rot + M_stream + M_accel reproduces M_tot within ∼5% for most regions (see Appendix C).

The above results basically support our conclusion for the cluster from the AMR simulations, but we do not have a statistical discussion including clusters from the SPH simulations because of the differences in resolutions and simulation methods. Instead, we divide the AMR cluster into two regions; upper and lower hemispheres with respect to the x–y plane. Then, we duplicate each hemisphere into one cluster. We call the synthetic cluster constructed from the z > 0 (z < 0) hemisphere "z +" ("z −"). Although these clusters are of course not independent of the original cluster and we cannot make a statistical argument on their properties, we can briefly look at the effects of substructures or the inhomogeneity of temperature and velocity field.

We repeat the same analysis on these two synthetic clusters, and the results are plotted in Figure 6. The amplitudes of the different terms vary significantly between the two clusters, and the degree of the validity of HSE is also very different. Nevertheless, the generic trend is clear; the relation of M_accel ≈ M_tot − M_therm holds almost independently of r.

Zoom In Zoom Out Reset image size

It is not clear, however, why the two hemispheres have such different values of (M_tot − M_therm)/M_tot; HSE holds very well for "z +," while this is not the case for "z −." The visual inspection of Figure 1 does not reveal any significant difference between the two. It may be because some local concentrations of dark matter enhance the acceleration/deceleration of gas and influence the overall non-sphericity of the gas density. Thus, an analysis that accounts for the ellipticity may provide a deeper insight into the validity of HSE, but it is beyond the scope of the present paper.

We define the distribution function f such that $f(\boldsymbol{x},\boldsymbol{v},t)d^3xd^3v$ is the probability that a randomly chosen particle in the system lies in the phase space volume d³xd³v at position $(\boldsymbol{x},\boldsymbol{v})$ and time t. The motion of such particles under the gravitational potential ϕ is described by the Boltzmann equation:

$$\begin{equation} \frac{\partial f}{\partial t}+v^i\frac{\partial f}{\partial x^i}-\frac{\partial \phi }{\partial x_i}\frac{\partial f}{\partial v^i} =\left(\frac{\delta f}{\delta t}\right)_{\rm {coll}}, \end{equation} \tag{ A1 }$$

where the collision term on the right-hand side takes into account the collisions between particles. Note that v_i (i = 1, 2, 3) represents a coordinate in the phase space and should not be confused with the velocity field at the spatial point xⁱ. For simplicity, we assume that all particles have the same mass m in the following.

First, we consider a collisionless case with (δf/δt)_coll = 0. Multiplying Equation (A1) by m and integrating it over the velocity space yields the continuity equation:

$$\begin{equation} \frac{\partial \rho }{\partial t}+\frac{\partial (\rho \bar{v}^i)}{\partial x^i}=0, \end{equation} \tag{ A2 }$$

where

$$\begin{equation} \rho (\boldsymbol{x},t)=\int d^3v\ mf(\boldsymbol{x},\boldsymbol{v},t) \end{equation} \tag{ A3 }$$

and we introduce the average over the velocity space:

$$\begin{equation} \bar{q}(\boldsymbol{x},t)=\frac{1}{\rho (\boldsymbol{x},t)}\int d^3v\ mq(\boldsymbol{x}, \boldsymbol{v}, t)f(\boldsymbol{x},\boldsymbol{v},t) \end{equation} \tag{ A4 }$$

for an arbitrary variable q such as v_i. Multiplying Equation (A1) by mv^j and integrating over the velocity space gives the momentum equations:

$$\begin{equation} \frac{\partial (\rho \bar{v}^i)}{\partial t}+\frac{\partial \tau ^{ij}_{\rm {J}}}{\partial x^j}=-\rho \frac{\partial \phi }{\partial x_i}, \end{equation} \tag{ A5 }$$

where

$$\begin{equation} \tau ^{ij}_{\rm {J}}=\rho \overline{v^iv^j}=\rho \sigma ^{2, ij}+\rho \bar{v}^i\bar{v}^j \end{equation} \tag{ A6 }$$

and σ^{2, ij} is the velocity dispersion tensor. Equations (A2) and (A5) reduce to the Jeans equations:

$$\begin{equation} \frac{\partial \bar{v}^i}{\partial t}+\bar{v^j}\frac{\partial \bar{v}^i}{\partial x^j}=-\frac{1}{\rho }\frac{\partial (\rho \sigma ^{2, ij})}{\partial x^j}-\frac{\partial \phi }{\partial x_i}. \end{equation} \tag{ A7 }$$

Next, we consider a collisional case. Rigorous handling of the collisional term is rather complicated and simplified models are often used. A conventional one is the Bhartnagar–Gross–Krock equation, which employs a linearized collisional term:

$$\begin{equation} \frac{\partial f}{\partial t}+v^i\frac{\partial f}{\partial x^i}-\frac{\partial \phi }{\partial x_i}\frac{\partial f}{\partial v^i}=-\frac{f-f_0}{\tau }, \end{equation} \tag{ A8 }$$

where τ is the relaxation time of the system considered, f₀ is the Maxwellian distribution function characterized by the local temperature $T(\boldsymbol{x}, t)$:

$$\begin{equation} f_0(\boldsymbol{x},\boldsymbol{v},t)=\left[\frac{m}{2\pi k_BT(\boldsymbol{x},t)}\right]^{3/2}\exp \left[-\frac{m(\boldsymbol{v}-\bar{\boldsymbol{v}}(\boldsymbol{x},t))^2}{2k_BT(\boldsymbol{x},t)}\right], \end{equation} \tag{ A9 }$$

and k_B is the Boltzmann constant. If we assume that mean values of conservatives such as mass and momentum are the same as those in local thermal equilibrium, then the collision term vanishes in the continuity and momentum equations. In local thermal equilibrium, pressure is defined from the diagonal components of $\sigma _{ij}^2$ by $p\delta _{ij} \equiv \rho \sigma _{ij}^2$ and the dispersion tensor can be written as

$$\begin{equation} \tau ^{ij}_{\rm {E}}= p\delta ^{ij} + \rho \bar{v}^i\bar{v}^j. \end{equation} \tag{ A10 }$$

If we replace τ_J in Equation (A5) with τ_E and combine them with Equation (A2), then we obtain the Euler equations:

$$\begin{equation} \frac{\partial \bar{v}^i}{\partial t}+\bar{v}^j\frac{\partial \bar{v}^i}{\partial x^j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}-\frac{\partial \phi }{\partial x_i}. \end{equation} \tag{ A11 }$$

The difference between the Euler and the Jeans equations resides only in the form of the dispersion tensor.

If we retain the off-diagonal components of the dispersion tensor, then they can be interpreted as viscosity, and the equations reduce to the Navier–Stokes equations (Choudhuri 1998; Chapman & Cowling 1970).

One can rewrite the continuity and momentum equations in the previous section in general coordinates

$$\begin{equation} \frac{\partial \rho }{\partial t}+\nabla _i(\rho \bar{v}^i)=0 \end{equation} \tag{ A12 }$$

and

$$\begin{equation} \frac{\partial (\rho \bar{v}^i)}{\partial t}+\nabla _j\tau ^{ij}=-\rho g^{ij}\nabla _j\phi , \end{equation} \tag{ A13 }$$

by using the covariant derivative operator ∇_i. For spherical coordinates (x¹ = r, x² = θ, x³ = φ), the non-zero components of the metric tensor g^ij are

$$\begin{equation} g_{11}=1,\quad g_{22}=r^2,\quad g_{33}=r^2\sin ^2\theta \end{equation} \tag{ A14 }$$

and the corresponding non-zero connection coefficients are

$$\begin{eqnarray} &\Gamma ^1_{22}=-r,\quad \Gamma ^1_{33}=-r\sin ^2\theta ,\quad \Gamma ^2_{12}=\Gamma ^2_{21}=\frac{1}{r}, \nonumber\\ &\Gamma ^2_{33}=-\sin \theta \cos \theta ,\quad \Gamma ^3_{13}=\Gamma ^3_{31}=\frac{1}{r},\quad \Gamma ^3_{23}=\Gamma ^3_{32}=\cot \theta. \nonumber\\ \end{eqnarray} \tag{ A15 }$$

The velocity vector is now given by $\bar{v}^i=(\bar{v}_r, \bar{v}_\theta /r, \bar{v}_\varphi /r\sin \theta)$.

In spherical coordinates, the continuity equation reads

$$\begin{equation} \frac{\partial \rho }{\partial t}+\frac{1}{r^2}\frac{\partial (r^2\rho \bar{v}_r)}{\partial r}+\frac{1}{r\sin \theta }\frac{\partial (\sin \theta \rho \bar{v}_\theta)}{\partial \theta }+\frac{1}{r\sin \theta }\frac{\partial (\rho \bar{v}_\varphi)}{\partial \varphi }=0. \end{equation} \tag{ A16 }$$

Setting $\tau ^{ij}=\tau ^{ij}_{\rm {E}}=\rho \bar{v}^i\bar{v}^j+pg^{ij}$ gives the Euler equations:

$$\begin{eqnarray} &&\left[\frac{\partial }{\partial t}+\bar{v}_r\frac{\partial }{\partial r}+\frac{\bar{v}_\theta }{r}\frac{\partial }{\partial \theta }+\frac{\bar{v}_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]\bar{v}_r-\frac{\bar{v}_\theta ^2+\bar{v}_\varphi ^2}{r}\nonumber\\ &&\quad =-\frac{1}{\rho }\frac{\partial p}{\partial r}-\frac{\partial \phi }{\partial r} \end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray} &&\left[\frac{\partial }{\partial t}+\bar{v}_r\frac{\partial }{\partial r}+\frac{\bar{v}_\theta }{r}\frac{\partial }{\partial \theta }+\frac{\bar{v}_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]\bar{v}_\theta +\frac{\bar{v}_r\bar{v}_\theta -\bar{v}_\varphi ^2\cot \theta }{r}\nonumber\\ &&\quad =-\frac{1}{\rho r}\frac{\partial p}{\partial \theta }-\frac{1}{r}\frac{\partial \phi }{\partial \theta } \end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray} &&\left[\frac{\partial }{\partial t}+\bar{v}_r\frac{\partial }{\partial r}+\frac{\bar{v}_\theta }{r}\frac{\partial }{\partial \theta }+\frac{\bar{v}_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]\bar{v}_\varphi +\frac{\bar{v}_r\bar{v}_\varphi +\bar{v}_\theta \bar{v}_\varphi \cot \theta }{r}\nonumber\\ &&\quad =-\frac{1}{\rho r\sin \theta }\frac{\partial p}{\partial \varphi }-\frac{1}{r\sin \theta }\frac{\partial \phi }{\partial \varphi }. \end{eqnarray} \tag{ A19 }$$

On the other hand, putting $\tau ^{ij}=\tau ^{ij}_{\rm {J}}=\rho \overline{v^iv^j}$ leads to the Jeans equations:

$$\begin{eqnarray} &&\left[\frac{\partial }{\partial t}+\bar{v}_r\frac{\partial }{\partial r}+\frac{\bar{v}_\theta }{r}\frac{\partial }{\partial \theta }+\frac{\bar{v}_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]\bar{v}_r\nonumber\\ &&\quad +\frac{1}{\rho }\left[\frac{\partial \left(\rho \sigma _{rr}^2\right)}{\partial r}+\frac{1}{r}\frac{\partial \left(\rho \sigma _{r\theta }^2\right)}{\partial \theta }+\frac{1}{r\sin \theta }\frac{\partial \left(\rho \sigma _{r\varphi }^2\right)}{\partial \varphi }\right] \nonumber\\ &&\quad +\frac{1}{r}\left(2\sigma _{rr}^2-\sigma _{\theta \theta }^2-\sigma _{\varphi \varphi }^2-\bar{v}_\theta ^2-\bar{v}_\varphi ^2+\sigma _{r\theta }^2\cot \theta \right) =-\frac{\partial \phi }{\partial r} \nonumber\\ \end{eqnarray} \tag{ A20 }$$

$$\begin{eqnarray} &&\left[\frac{\partial }{\partial t}+\bar{v}_r\frac{\partial }{\partial r}+\frac{\bar{v}_\theta }{r}\frac{\partial }{\partial \theta }+\frac{\bar{v}_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]\bar{v}_\theta\nonumber\\ &&\quad +\frac{1}{\rho }\left[\frac{\partial \left(\rho \sigma _{r\theta }^2\right)}{\partial r}+\frac{1}{r}\frac{\partial \left(\rho \sigma _{\theta \theta }^2\right)}{\partial \theta }+ \frac{1}{r\sin \theta }\frac{\partial \left(\rho \sigma _{\theta \varphi }^2\right)}{\partial \varphi }\right] \nonumber\\ &&\quad+\frac{1}{r}\left(3\sigma _{r\theta }^2-\sigma _{\varphi \varphi }^2\cot \theta +\bar{v}_r\bar{v}_\theta -\bar{v}_\varphi ^2\cot \theta +\sigma _{\theta \theta }^2\cot \theta \right)\nonumber\\ &&\quad =-\frac{1}{r}\frac{\partial \phi }{\partial \theta } \end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray} &&\left[\frac{\partial }{\partial t}+\bar{v}_r\frac{\partial }{\partial r}+\frac{\bar{v}_\theta }{r}\frac{\partial }{\partial \theta }+\frac{\bar{v}_\varphi }{r\sin \theta }\frac{\partial }{\partial \varphi }\right]\bar{v}_\varphi\nonumber\\ && \quad+\frac{1}{\rho }\left[\frac{\partial \left(\rho \sigma _{r\varphi }^2\right)}{\partial r}+\frac{1}{r}\frac{\partial \left(\rho \sigma _{\theta \varphi }^2\right)}{\partial \theta }+\frac{1}{r\sin \theta }\frac{\partial \left(\rho \sigma _{\varphi \varphi }^2\right)}{\partial \varphi }\right] \nonumber\\ &&\quad+\frac{1}{r}\left(3\sigma _{r\varphi }^2+\bar{v}_r\bar{v}_\varphi +\bar{v}_\theta \bar{v}_\varphi \cot \theta +2\sigma _{\theta \varphi }^2\cot \theta \right)\nonumber\\ &&\quad =-\frac{1}{r\sin \theta }\frac{\partial \phi }{\partial \varphi }. \end{eqnarray} \tag{ A22 }$$