ADAPTIVELY REFINED LARGE EDDY SIMULATIONS OF A GALAXY CLUSTER: TURBULENCE MODELING AND THE PHYSICS OF THE INTRACLUSTER MEDIUM

The dynamics of a compressible, viscous, self-gravitating fluid with density ρ(r_i, t), momentum density ρv_i(r_i, t), and total energy density ρe(r_i, t) at spatial position (r₁, r₂, r₃) is given by the following set of equations:

$$\begin{eqnarray} \frac{\partial }{\partial t}\rho + \frac{\partial }{\partial r_j}(v_j \rho) &= 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}(\rho v_i) + \frac{\partial }{\partial r_j}(v_j \rho v_i) &= -\frac{\partial }{\partial r_i}p + \frac{\partial }{\partial r_j}\sigma ^{\prime }_{ij}+\rho g_i, \end{eqnarray}\vspace*{-1pc} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}(\rho e) + \frac{\partial }{\partial r_j}(v_j \rho e) &= -\frac{\partial }{\partial r_j}(v_j p) + \frac{\partial }{\partial r_j}(v_i \sigma ^{\prime }_{ij})+ v_i \rho g_i,\hspace{9pt} \end{eqnarray} \tag{ 3 }$$

where p is the pressure, g_i the gravitational acceleration, and σ'_ij the viscous stress tensor. Note that the Einstein sum convention applies to repeated indices.

As shown by Schmidt et al. (2006), these equations can be decomposed into large-scale (resolved) and small-scale (unresolved) parts using the filter formalism proposed by Germano (1992) in terms of density-weighted quantities.⁴ By means of filtering, any field quantity a can be split into a smoothed part 〈a〉 and a fluctuating part a', where 〈a〉 varies only at scales greater than the prescribed filter length. We define density-weighted filtered quantities according to Favre (1969) by

$$\begin{eqnarray} &&\langle \rho a \rangle = \langle \rho \rangle \hat{a} \Rightarrow \hat{a}=\frac{\langle \rho a \rangle }{\langle \rho \rangle }. \end{eqnarray} \tag{ 4 }$$

Following Schmidt et al. (2006), filtered equations for compressible fluid dynamics can be derived:

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \rho \rangle + \frac{\partial }{\partial r_j}\hat{v}_j\langle \rho \rangle =&\ 0, \end{eqnarray}\vspace*{-1pc} \tag{ 5 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \rho \rangle \hat{v}_i + \frac{\partial }{\partial r_j}\hat{v}_j\langle \rho \rangle \hat{v}_i\,=& {-}\frac{\partial }{\partial r_i}\langle p \rangle +\frac{\partial }{\partial r_j}\langle \sigma ^{\prime }_{ij} \rangle \nonumber\\ &+\langle \rho \rangle \hat{g}_i-\frac{\partial }{\partial r_j}\hat{\tau }(v_i,v_j), \end{eqnarray}\vspace*{-1pc} \tag{ 6 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \rho \rangle e_{\mathrm{res}}+\frac{\partial }{\partial r_j}\hat{v}_j\langle \rho \rangle e_{\mathrm{res}}\,{=}&{-}\frac{\partial }{\partial r_i}\hat{v}_i\langle p \rangle +\frac{\partial }{\partial r_j}\hat{v}_i\langle \sigma ^{\prime }_{ij} \rangle\hspace{30pt} \nonumber\\ &+\langle \rho \rangle (\lambda +\epsilon)-\hat{v}_i\frac{\partial }{\partial r_j}\hat{\tau }(v_i,v_j)\nonumber\\ &+\langle \rho \rangle \hat{v}_i\hat{g}_i-\frac{\partial }{\partial r_j}\hat{\tau }(v_j,e_{\mathrm{int}}), \end{eqnarray} \tag{ 7 }$$

where we introduced the total resolved energy $e_{\mathrm{res}}=\hat{e}_{\mathrm{int}}+\frac{1}{2}\hat{v_i}\hat{v_i}$, $\hat{e}_{\mathrm{int}}$ being the filtered internal energy, and the generalized moments that are generically defined by

$$\begin{equation} \hat{\tau }(a,b)=\langle \rho a b \rangle - \langle \rho \rangle \hat{a} \hat{b}, \end{equation}\vspace*{-2pc} \tag{ 8 }$$

$$\begin{eqnarray} \hat{\tau }(a,b,c)=&\ \langle \rho a b c \rangle - \langle \rho \rangle \hat{a} \hat{b} \hat{c} \nonumber\\ &-\hat{a}\hat{\tau }(b,c)-\hat{b}\hat{\tau }(a,c)-\hat{c}\hat{\tau }(a,b), \end{eqnarray}\vspace*{-1pc} \tag{ 9 }$$

$$\begin{equation} \hat{\tau }(a,b,c,d)=\ldots \end{equation} \tag{ 10 }$$

for Favre-filtered quantities a, b, c, etc. Germano interpreted the trace of $\hat{\tau }(v_i,v_j)\langle \rho \rangle$as the squared velocity fluctuation, $q^{2}:=\hat{\tau }(v_i,v_i)/\langle \rho \rangle$. The evolution of the corresponding turbulent energy, $e_{\mathrm{t}}=\frac{1}{2}q^{2}$, is given by

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\langle \rho \rangle e_{\mathrm{t}}+\frac{\partial }{\partial r_j}\hat{v}_j \langle \rho \rangle e_{\mathrm{t}}=\ \mathbb {D}+\Sigma +\Gamma -\langle \rho \rangle (\lambda +\epsilon), \end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray} \mathbb {D} &= -\frac{\partial }{\partial r_j}{ \left[ \frac{1}{2}\hat{\tau }(v_j,v_i,v_i) -\mu +\kappa \right] }, \end{eqnarray}\vspace*{-1pc} \tag{ 12 }$$

$$\begin{eqnarray} \Sigma &= -\hat{\tau }_{ij} \partial v_i / \partial r_j, \end{eqnarray}\vspace*{-1pc} \tag{ 13 }$$

$$\begin{eqnarray} \Gamma &=+\hat{\tau }(v_i,g_i), \end{eqnarray}\vspace*{-1pc} \tag{ 14 }$$

$$\begin{eqnarray} \langle \rho \rangle \lambda &={ \left[ \langle p \rangle \frac{\partial }{\partial r_i}\hat{v}_i-\langle p \frac{\partial }{\partial r_i} v_i \rangle \right] }, \end{eqnarray}\vspace*{-1pc} \tag{ 15 }$$

$$\begin{eqnarray} \langle \rho \rangle \epsilon &=-{ \left[ \langle \sigma ^{\prime }_{ij} \rangle \frac{\partial }{\partial r_j}\hat{v}_i-\langle \sigma ^{\prime }_{ij} \frac{\partial }{\partial r_j} v_i \rangle \right] }, \end{eqnarray}\vspace*{-1pc} \tag{ 16 }$$

$$\begin{eqnarray} {-}\mu &=\langle v_i p \rangle -\hat{v}_i\langle p \rangle, \end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray} {-}\kappa &=\langle v_i \sigma ^{\prime }_{ij} \rangle -\hat{v}_i\langle \sigma ^{\prime }_{ij} \rangle. \end{eqnarray} \tag{ 18 }$$

The explicit forms of the quantities $\mathbb {D},\lambda,\epsilon,\Gamma$, and $\hat{\tau }(v_i,v_j)$ are unknown and have to be modeled in terms of closure relations, i.e., functions of the filtered flow quantities (or their derivatives) and the turbulent energy, e_t. The closures for all these terms represent the SGS model.

In the following, we consider a simplified set of equations to model the influence of the turbulent small-scale (SGS) motions on the numerically resolved scales ℓ ⩾ ℓ_Δ, neglecting the influence of the viscous stress tensor 〈σ'_ij〉 (which is a very good approximation for high Reynolds numbers) and the turbulent transport of heat given by the divergence of $\hat{\tau }(v_j,e_{\mathrm{int}})$ in Equations (6) and (7). Moreover, gravitational effects on unresolved scales are neglected, i.e., we set Γ = 0 in Equation (11) for the turbulent energy. For the terms (12), (13), and (16), we adopt SGS closures that have been applied in LES of incompressible turbulence. The numerical study by Schmidt et al. (2006) demonstrated that these closures can be carried over to transonic turbulence, for which the unresolved turbulent velocity fluctuations are small compared to the speed of sound. Additionally, we utilize the pressure-dilatation model of Sarkar (1992) in order to account for moderate compressibility effects. Since we concentrate on the dynamics of the gas in the ICM, where the Mach numbers may locally approach unity, the SGS closures outlined subsequently serve as a reasonable approximation. In supersonic flow regions, on the other hand, the SGS model is deactivated in order to maintain stability (see Section 2.4).

The flux of kinetic energy from resolved scales toward SGSs, i.e., the rate of turbulent energy production, is given by the contraction of the turbulent stress tensor and the Jacobian of the resolved velocity field. Since $\hat{\tau }_{kk}=\langle \rho \rangle q^2$, we split the tensor in a symmetric trace-free part $\hat{\tau }^*_{ij}$ and a diagonal part:

$$\begin{eqnarray} &&\hat{\tau }_{ij}=\hat{\tau }^*_{ij}+\frac{1}{3}\delta _{ij}\langle \rho \rangle q^2. \end{eqnarray} \tag{ 19 }$$

The model for $\hat{\tau }^*_{ij}$ is based on the turbulent viscosity hypothesis (Boussinesq 1877), which means that $\hat{\tau }^*_{ij}$ is assumed to be of the same form as the stress tensor σ'_ij of a Newtonian fluid. Hence,

$$\begin{eqnarray} &&\hat{\tau }^*_{ij}=- 2 \eta _{\mathrm{t}} S^*_{ij} \end{eqnarray} \tag{ 20 }$$

with a turbulent dynamic viscosity η_t = 〈ρ〉ν_t = 〈ρ〉C_νl_Δq and

$$\begin{eqnarray} &&S^*_{ij}=\frac{1}{2}{ \left(\frac{\partial }{\partial r_j}\hat{v}_i+\frac{\partial }{\partial r_i}\hat{v}_j \right) } -\frac{1}{3}\delta _{ij}\frac{\partial }{\partial r_k}\hat{v}. \end{eqnarray} \tag{ 21 }$$

The turbulence production term is therefore modeled as

$$\begin{eqnarray} &&\hat{\tau }(v_i,v_j)=-2\langle \rho \rangle C_{\nu } l_{\Delta } q S^*_{ij}+\frac{1}{3}\delta _{ij}\langle \rho \rangle q^2. \end{eqnarray} \tag{ 22 }$$

We set C_ν = 0.05 (Sagaut 2006).

The SGS transport of turbulent energy (Equation (12)) is modeled by a gradient-diffusion hypothesis, stating that the non-linear term is proportional to the turbulent velocity q² gradient (Sagaut 2006)

$$\begin{eqnarray} &&\mathbb {D}=\frac{\partial }{\partial r_i}C_{\mathbb {D}}\langle \rho \rangle l_{\Delta } q^2 \frac{\partial }{\partial r_i} q. \end{eqnarray} \tag{ 23 }$$

The diffusion coefficient has been calibrated to $C_{\mathbb {D}} \approx 0.4$ by numerical experiments (Schmidt et al. 2006).

For sufficiently high Reynolds numbers, viscous energy dissipation (Equation (16)) becomes entirely an SGS effect. The simplest expression that can be built from the characteristic turbulent velocity and length scale for dissipation is

$$\begin{eqnarray} &&\epsilon =C_{\epsilon }\frac{q^3}{l_{\Delta }}. \end{eqnarray} \tag{ 24 }$$

For our simulations, we set C = 0.5 (Sagaut 2006).

The effect of unresolved pressure fluctuations in compressible turbulence is described by the λ term (Equation (15)). A simple closure for subsonic turbulent flow is (Deardorff 1973)

$$\begin{eqnarray} &&\lambda =C_{\lambda } q^2 \frac{\partial }{\partial r_i}\hat{v}_i, \end{eqnarray} \tag{ 25 }$$

where $C_{\lambda }=-\frac{1}{5}$ (Fureby et al. 1997). Sarkar (1992) performed simulations of simple compressible flows and investigated the influence of the mean Mach number of the flow on the turbulent dissipation, , and the pressure dilatation, λ. Based on this analysis, he suggested different models for these terms, which we will describe in the following sections. These modifications have been proven to yield good results for transonic turbulence (Shyy & Krishnamurty 1997).

As a major effect of compressibility from direct numerical simulation (DNS), Sarkar (1992) identified that the growth rate of kinetic energy decreases if the initial turbulent Mach number increases. This means that the dissipation of kinetic energy (and, therefore, of the turbulent energy) increases with the turbulent Mach number M_t = q/c_s, where c_s is the speed of sound. Sarkar (1992) suggested to account for this effect by using

$$\begin{eqnarray} &&\epsilon =C_{\epsilon }\frac{q^3}{l_{\Delta }}\big(1+\alpha _1 M_{\mathrm{t}}^2\big) \end{eqnarray} \tag{ 26 }$$

with α₁ = 0.5 as a model for the dissipation of turbulent energy.

Based on a decomposition of all variables of the equation for instantaneous pressure,

$$\begin{eqnarray} &&\frac{\partial ^2}{\partial r_i^2} p = \frac{\partial ^2}{\partial t^2}\rho - \frac{\partial ^2}{\partial r_i \partial r_j}(\rho v_i v_j-\sigma ^{\prime }_{ij}), \end{eqnarray} \tag{ 27 }$$

into a mean and a fluctuating part and comparisons with DNSs of simple compressible flows, Sarkar (1992) proposed a different model for the pressure dilatation

$$\begin{eqnarray} &&\lambda =\alpha _2 M_{\mathrm{t}} \hat{\tau }^*_{ij} \frac{\partial \hat{v}_i}{\partial r_j} - \alpha _3 M_{\mathrm{t}}^2 C_{\epsilon }\frac{q^3}{l_{\Delta }} - 8 \alpha _4 M_{\mathrm{t}}^2 p_{\mathrm{t}} \frac{\partial \hat{v}_k}{\partial r_k} \end{eqnarray} \tag{ 28 }$$

with α₂ = 0.15 and α₃ = 0.2 obtained from a curve fit of the model with DNS. Unfortunately, Sarkar (1992) does not specify a value for α₄, so there is some confusion in the literature about it. For example, Shyy & Krishnamurty (1997) set α₄ = 0 and still found the Sarkar model in good agreement with their DNS simulation. In this work, we adopt α₄ = α²₂/2. With this choice, the effective production of turbulent energy vanishes for a turbulent Mach number, M_t = 1/α₂. In the following, we will sometimes refer to the "Sarkar SGS" when the Equations (26) and (28) are used in the model.

Simulations with a comoving cosmological background require a formulation of the filtered fluid dynamical equations in comoving coordinates. Applying the Germano decomposition (Section 2.1) in a comoving coordinate system with spatially homogeneous scale factor a(t), we obtain⁵

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \tilde{\rho } \rangle + \frac{1}{a}\frac{\partial }{\partial x_j}\hat{u}_j \langle \tilde{\rho } \rangle =&\ 0, \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \tilde{\rho } \rangle \hat{u}_i +\frac{1}{a}\frac{\partial }{\partial x_j}\hat{u}_j \langle \tilde{\rho } \rangle \hat{u}_i& =& {-} \frac{1}{a}\frac{\partial }{\partial x_i}\langle \tilde{p} \rangle +\langle \tilde{\rho } \rangle \hat{g}^*_i\nonumber\\ &&{-} \frac{1}{a}\frac{\partial }{\partial x_j}\hat{\tau }(u_i,u_j)- \frac{\dot{a}}{a}\langle \tilde{\rho } \rangle \hat{u}_i,\nonumber\\ \end{eqnarray}\vspace*{-1.50pc} \tag{ 30 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \tilde{\rho } \rangle e_{\mathrm{res}}+\frac{1}{a}\frac{\partial }{\partial x_j}\hat{u}_j\langle \tilde{\rho } \rangle e_{\mathrm{res}}& =& {-}\frac{1}{a}\frac{\partial }{\partial x_i}\hat{u}_i\langle \tilde{p} \rangle -\frac{1}{a}\langle \tilde{\rho } \rangle \hat{u}_i\hat{g}^*_i\nonumber\\ &&{-}\frac{\dot{a}}{a}(\langle \tilde{\rho } \rangle e_{\mathrm{res}} +\frac{1}{3}\langle \tilde{\rho } \rangle \hat{u}_i\hat{u}_i+\langle \tilde{p} \rangle)\nonumber\\ &&+\langle \tilde{\rho } \rangle (\lambda +\epsilon) -\frac{1}{a}\hat{u}_i\frac{\partial }{\partial x_j}\hat{\tau }(u_i,u_j),\nonumber\\ \end{eqnarray}\vspace*{-1.50pc} \tag{ 31 }$$

$$\begin{eqnarray} \frac{\partial }{\partial t}\langle \tilde{\rho } \rangle e_{\mathrm{t}}+\frac{1}{a}\frac{\partial }{\partial x_j}\hat{u}_j\langle \tilde{\rho } \rangle e_{\mathrm{t}}&=& \mathbb {D} +\Gamma -\langle \tilde{\rho } \rangle (\lambda +\epsilon)\nonumber\\ &&{-}\frac{1}{a}\hat{\tau }(u_j,u_i)\frac{\partial }{\partial x_j}\hat{u}_i -2\frac{\dot{a}}{a}\langle \tilde{\rho } \rangle e_{\mathrm{t}}.\nonumber\\ \end{eqnarray} \tag{ 32 }$$

With respect to the non-comoving equations listed in Section 2.1, the only term to additionally implement is the last one on the right-hand side of Equation (32). Furthermore, SGS closures have to be expressed in terms of the Jacobian of the velocity in comoving coordinates,

$$\begin{eqnarray} J_{ij} = \frac{\partial }{\partial r_i}v_j&=\frac{1}{a}\frac{\partial }{\partial x_i}u_j + \frac{\dot{a}}{a}\delta _{ij}. \end{eqnarray} \tag{ 33 }$$

In particular, the trace-free rate of strain tensor in comoving coordinates is given by

$$\begin{eqnarray} &&S^*_{ij}=\frac{1}{2}{ (J_{ij}+J_{ji}) }-\frac{1}{3}\delta _{ij}J_{kk}. \end{eqnarray} \tag{ 34 }$$

Numerical difficulties result from constraints on the validity of SGS closures. In particular, the turbulent viscosity hypothesis expressed by Equation (22) was devised to account for the production of turbulence by shear in moderately compressible flow. This is typically encountered in the dense, central regions of galaxy clusters. However, the surrounding low-density gas can be accelerated very quickly in the gravitational field of the cluster. Moreover, high-velocity gradients are encountered in the vicinity of shocks, which are produced by gas accretion onto filaments, sheets, and halos as well as by the merging of substructures. In order to inhibit unphysical production of turbulent energy by these mechanisms, which are not accommodated in the present formulation of the SGS model, several numerical safeguarding mechanisms have been introduced.

First of all, we implemented a simple shock detector, which identifies strong negative divergence. A cell is marked if the velocity jump corresponding to the negative divergence becomes greater than the sound speed across the cell width l_Δ:

$$\begin{eqnarray} &&{-}\frac{\partial v_i}{\partial r_i} > \frac{c_{\mathrm{s}}}{l_{\Delta }}. \end{eqnarray} \tag{ 35 }$$

In cells satisfying the above criterion, the source and transport terms of the SGS model (Equations (22) and (23)) are disabled. The turbulent energy is only advected in these cells, and no coupling to the velocity and the resolved energy is applied.

Besides the previous check, an additional constraint is imposed on the magnitude of the turbulent energy, via the turbulent Mach number. Basically, the SGS model breaks down once M_t becomes large compared to unity, therefore a threshold for this quantity is set to $M_\mathrm{{t,max}}=\sqrt{2}$. This value for the maximal turbulent Mach number is motivated by the theory of isothermal turbulence, where the effective gas pressure can be expressed as

$$\begin{eqnarray} &&p_{\mathrm{eff}}=\rho c_{\mathrm{s}}^2 + \frac{1}{3}\rho q^2 = \gamma _{\rm {eff}} \rho c_{\mathrm{s}}^2, \end{eqnarray} \tag{ 36 }$$

and, consequently, p_eff is limited to the adiabatic value $\frac{5}{3}\rho c_{\mathrm{s}}^2$. We verified that this threshold does not harm our cluster simulations, because in the hotter gas phases (T > 10⁵ K) turbulence is largely subsonic, and the threshold is rarely reached.

A supplementary low-temperature cutoff ensures that the sound speed does not drop to excessively low values, which occur in cosmological simulations especially in the low-density voids. We set the lower limit of the temperature to T_min = 10 K. This threshold ensures numerical stability and does not affect the baryon physics appreciably, apart from possibly making the shocks on accreting structures weaker.

In LES, the filtered Equations (29)–(32) are solved using an SGS model as outlined in the previous section. However, the closure relations we use and, in fact, the very concept of SGS turbulence energy only applies if the velocity fluctuations on SGSs are nearly isotropic. This limits the LES methodology to flows where all anisotropies stemming from large-scale features, like boundary conditions or external forces, can be resolved. In the fearless method, the grid resolution ℓ_Δ is locally adjusted by AMR in order to ensure that the anisotropic, energy-containing scales are resolved everywhere. On the other hand, it is assumed that turbulence is asymptotically isotropic on length scales comparable to or less than ℓ_Δ. It is very difficult to justify the latter assumption a priori, because there are no refinement criteria that would guarantee asymptotic isotropy on the smallest resolved length scales. By careful analysis of simulation results, however, one can gain confidence that AMR resolves turbulent regions appropriately.

As an infrastructure for the implementation of fearless, we chose Enzo v. 1.0 (O'Shea et al. 2005), an AMR, grid-based hybrid (hydrodynamics plus N-body) code based on the PPM solver (Colella & Woodward 1984) and especially designed for simulations of cosmological structure formation. When a grid location is flagged for refinement in Enzo, a new finer grid is created, and the cell values on the finer grid are generated by interpolating them from the coarser grid using a conservative interpolation scheme. At each time step of the coarse grid, the values from the fine grid are averaged and the values computed on the coarse grid (in the region where fine and coarse grid overlap) are replaced. However, this approach does not account for the inherent scale dependence of the turbulent energy. Assuming Kolmogorov scaling (Kolmogorov 1941; Frisch 1995), the turbulent energies at two different levels of refinement with cell size l_Δ,1 and l_Δ,2, respectively, are statistically related by

$$\begin{eqnarray} &&\frac{e_{\mathrm{t},1}}{e_{\mathrm{t},2}}=\frac{q_1^2}{q_2^2}\sim { \left(\frac{l_{\Delta,1}}{l_{\Delta,2}} \right) }^{2/3}. \end{eqnarray} \tag{ 37 }$$

The assumption of Kolmogorov scaling has a caveat: since grid regions that are flagged for refinement cannot be preconditioned to ensure local isotropy, the above scaling relation is an estimate that applies to regions in which turbulence is fully developed. Near pronounced inhomogeneities or in regions undergoing rapid temporal evolution deviations from Kolmogorov scaling are inevitable. For this reason, there is still room for improvement of the method.

Using the above scaling relation, we implemented a simple algorithm to adjust the turbulent energy budget when grids are refined or derefined. The following procedure is used once a grid is refined.

• 1.  
  Interpolate the values from the coarse to the fine grid using the standard interpolation scheme from Enzo.
• 2.  
  On the finer grid, correct the values of the velocity components, $\hat{v}_i$, and the turbulent energy, $e_{\mathrm{t}}=\frac{1}{2}q^2$, as follows

  $$\begin{eqnarray} \hat{v}^{\prime }_i&=\hat{v}_i \sqrt{1+\frac{e_{\mathrm{t}}}{\hat{e}_{\mathrm{kin}}}{ \big(1-r_{\Delta }^{-2/3} \big) }}, \end{eqnarray}\vspace*{-1pc} \tag{ 38 }$$

  $$\begin{eqnarray} e^{\prime }_{\mathrm{t}}&=e_{\mathrm{t}} r_{\Delta }^{-2/3}, \end{eqnarray} \tag{ 39 }$$

  where $\hat{e}_{\mathrm{kin}}$ is the resolved kinetic energy, r_Δ is the refinement factor of the mesh, and the primed quantities are the final values on the fine grid. The resolved energy is adjusted such that the sum of resolved energy and turbulent energy remains conserved.

Apart from adjusting the energy budget, the resolved flow should feature velocity fluctuations on length scales smaller than the cutoff length of the parent grid if a refined grid is generated. To address this problem, we observe that the smallest pre-existing eddies that are inherited from the parent grid will produce new eddies of smaller size within a turnover time. Although this implies a small delay because of the higher time resolution of the refined grid, the flow will rapidly adjust itself to the new grid resolution.

For grid derefinement we reverse this procedure.

• 1.  
  Average the values from the fine grid and replace the corresponding values on the coarse grid.
• 2.  
  In the regions of the coarse grid covered by finer grids, correct the velocity components and the turbulent energy:

  $$\begin{eqnarray} e^{\prime }_{\mathrm{t}}&=e_{\mathrm{t}}r_{\Delta }^{2/3}, \end{eqnarray}\vspace*{-1pc} \tag{ 40 }$$

  $$\begin{eqnarray} \hat{v}^{\prime }_i&=\hat{v}_i \sqrt{1-\frac{e_{\mathrm{t}}}{\hat{e}_{\mathrm{kin}}}{ \big (r_{\Delta }^{2/3}-1 \big) }}. \end{eqnarray} \tag{ 41 }$$

  Here, primed quantities denote the final values on the coarse grid. The resolved energy is adjusted to maintain energy conservation and a positive kinetic energy.

For global energy conservation, it turned out to be important to compute the turbulent stress term in Equation (6) indirectly as

$$\begin{eqnarray} &&\hat{v}_i\frac{\partial }{\partial r_j}\hat{\tau }(v_i,v_j) = \frac{\partial }{\partial r_j}\hat{v}_i\hat{\tau }(v_i,v_j) - \hat{\tau }(v_j,v_i)\frac{\partial }{\partial r_j}\hat{v_i}. \end{eqnarray} \tag{ 42 }$$

The reason behind it is that only by using this rearrangement of the terms we can ensure that we do not introduce small numerical errors which would violate the local sum,

$$\begin{eqnarray} &&\frac{\partial }{\partial r_j}\hat{v}_i\hat{\tau }(v_i,v_j) = \hat{v}_i\frac{\partial }{\partial r_j}\hat{\tau }(v_i,v_j)+\hat{\tau }(v_j,v_i)\frac{\partial }{\partial r_j}\hat{v_i}, \end{eqnarray} \tag{ 43 }$$

leading to a big error in global energy conservation.

As a testing case, we run a LES of driven turbulence as outlined above on a static grid of 256³ grid points and periodic boundary conditions. The adiabatic index $\gamma =\frac{5}{3}$ and M_f = 0.68.

Figure 1 shows the typical time development of the mass-weighted mean energies in our simulation including the energy e_f injected into the system by random forcing. It is evident from the curve of the turbulent energy that after one integral timescale, our simulation reaches an equilibrium between production and dissipation of turbulent energy.

Zoom In Zoom Out Reset image size

In Figure 2, we plot the time development of the relative error $\frac{\Delta e(t)}{e(0)}$ of the mean total energy, defined as the sum of the mass-weighted means of internal energy, kinetic energy, turbulent energy minus the injected energy by the forcing,

$$\begin{eqnarray} &&\hat{e}_{\mathrm{tot}}= \hat{e}_{\mathrm{int}}+\hat{e}_{\mathrm{kin}}+e_{\mathrm{t}}-\hat{e}_{\mathrm{f}}, \end{eqnarray} \tag{ 44 }$$

where $\hat{e} = \frac{\langle \rho e \rangle }{\langle \rho \rangle }$. It demonstrates that with our basic model, the relative error in energy is comparable to the error without SGS model, and is around 1%. The energy conservation of the model using the Sarkar modifications is equally good.

Zoom In Zoom Out Reset image size

It is also instructive to plot the difference between the energy contributions (internal and kinetic energy) in the simulations with and without the SGS model. These differences are shown in Figure 3. One can conclude from this figure that, at the beginning of the turbulent driving, the turbulent energy produced in our simulation with SGS model is found in the kinetic energy of the simulation without SGS model. In contrast, from t = 1.2 t_int on, most of the turbulent energy can be found in the internal energy of the simulation without SGS model. Turbulent energy can therefore be interpreted as a kind of buffer which prevents the kinetic energy in our simulation to be converted instantly into thermal energy.

Zoom In Zoom Out Reset image size

A necessary condition for the validity of the turbulent energy transfer algorithms explained in Section 2.5 is that an AMR simulation should approximately reproduce the results of static grid simulations corresponding to the different levels of refinement. To test this, we compared an AMR simulation of driven turbulence with a 32³ root grid resolution and three additional levels (with a refinement factor of 2 between each level) to three static grid simulations with resolutions of 32³, 64³, 128³, and 256³. In order to allow for the comparison of averaged quantities, refinement of the entire domain was enforced at all levels of the AMR simulation. In order to reach a statistically stationary root-mean-square (rms) Mach number, we ran the simulations for nearly isothermal gas with γ = 1.01 (for adiabatic turbulence, after an initial rise, the rms Mach number gradually decreases with time because of the dissipative heating of the gas). We used a supersonic forcing Mach number, M_f = 2.7, to check whether a consistent turbulent energy budget could be achieved for highly compressible turbulence, albeit the scaling relations (Equation (37)) of incompressible turbulence were utilized.

The results of this consistency check can be seen in Figure 4. We observe that the time development of the mean turbulent energy is very similar on the different levels of the AMR simulation compared to the static grid simulations, except for some deviations at the root level. On the other hand, comparing these results to a simulation without correcting turbulent energy at grid refinement/derefinement (Figure 5), it is evident that the scaling of the turbulent energy in the latter case is inconsistent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We performed simulations of cluster formation with Enzo, following Iapichino & Niemeyer (2008). We will compare two runs: one of them was done with the public version of Enzo (without SGS model), and the other with fearless implemented, in its version including the Sarkar correction (Equations (26) and (28)). The simulations were done using a flat ΛCDM background cosmology with a dark energy density Ω_Λ = 0.7, a total (including baryonic and dark matter) matter density Ω_m = 0.3, a baryonic matter density Ω_b = 0.04, the Hubble parameter set to h = 0.7, the mass fluctuation amplitude σ₈ = 0.9, and the scalar spectral index n = 1. Both simulations were started with the same initial conditions at redshift z_ini = 60, using the Eisenstein & Hu (1999) transfer function, and evolved to z = 0. The simulations are adiabatic with a heat capacity ratio γ = 5/3 assuming a fully ionized gas with a mean molecular weight m_μ = 0.6 u. Cooling physics, magnetic fields, feedback, and transport processes are neglected.

The simulation box has a comoving size of 128 Mpc h⁻¹. It is resolved with a root grid (level l = 0) of 128³ cells and 128³ N-body particles. A static child grid (l = 1) is nested inside the root grid with a size of 64 Mpc h⁻¹, 128³ cells and 128³ N-body particles. The mass of each particle in this grid is 9 × 10⁹ M_☉ h⁻¹. Inside this grid, in a volume of 38.4 Mpc h⁻¹, adaptive grid refinement from level l = 2 to l = 7 is enabled using an overdensity refinement criterion as described in Iapichino & Niemeyer (2008) with an overdensity factor f = 4.0. The refinement factor between two levels was set to r_Δ = 2, allowing for an effective resolution of 7.8 kpc h⁻¹.

The static and dynamically refined grids were nested around the place of formation of a galaxy cluster, identified using the HOP algorithm (Eisenstein & Hut 1998). Since the realization of the initial conditions was chosen identical to Iapichino & Niemeyer (2008), this study is based on the same cluster analyzed in that work. The cluster has a virial mass of M_vir = 5.95 × 10¹⁴ M_☉ h⁻¹ and a virial radius of R_vir = 1.37 Mpc h⁻¹.

In static grid simulations, one often chooses to use the grid resolution l_Δ as characteristic length scale to compute a characteristic velocity or eddy turnover time for this scale. However, in an AMR code it is not trivial to compute the turbulent velocity q_l associated with a characteristic length scale l = l_Δ, since l_Δ varies in time and space. To circumvent this difficulty, we assume that below the grid resolution, turbulent velocity locally scales according to Kolmogorov

$$\begin{eqnarray} &&q(l) \sim l^{1/3}. \end{eqnarray} \tag{ 45 }$$

We thereby assume that locally a Kolmogorov-like energy cascade sets in, at a length scale given by the resolution of the grid at that position. This local hypothesis holds here only for the analysis of our simulations, and is similar to the assumption done in Section 2.5 for managing the grid refinement and derefinement.

As a characteristic scale of our analysis, we choose the length scale of our highest resolved regions, which is l_min = 7.8 kpc h⁻¹. The turbulent velocity in the most finely resolved regions can be computed directly from the values of the turbulent energy $q(l) = \sqrt{2 e_{\mathrm{t}}}$ on the grid; the turbulent velocity in less finely refined regions is scaled down according to our local Kolmogorov hypothesis as

$$\begin{eqnarray} &&q(l_{\mathrm{min}}) = q(l_{\Delta }) { \left(\frac{l_{\mathrm{min}}}{l_ {\Delta }} \right) }^{2/3}. \end{eqnarray} \tag{ 46 }$$

In Section 3.2, we studied the temporal evolution of the turbulent energy at different resolutions in a simulation of driven turbulence. In this section, we repeat this analysis for our fearless cluster simulation. Figure 6 shows the evolution of the mass-weighted mean turbulent energy for every level of our AMR simulation. We see from the plot that the turbulent energy on the higher AMR levels l (meaning at smaller scales) is higher at early times (2 Gyr < t < 6 Gyr). Later this picture changes, but not completely. For example, the turbulent energy at l = 4 stays above the turbulent energy at l = 3 throughout the simulation. The magnitude of e_t along the AMR levels suggests to locate the turbulence injection length scale between 125 and 250 kpc h⁻¹, corresponding to the effective resolutions of levels 3 and 4. This is only a rough qualitative estimate, but nevertheless in agreement with theoretical expectations (e.g., Subramanian et al. 2006).

Zoom In Zoom Out Reset image size

Particularly noteworthy are the turbulent energy fluctuations on smaller scales at the time 2 Gyr < t < 6 Gyr, corresponding to a redshift z = 3–1. We can interpret these large fluctuations as evidence for violent major mergers that happen at that time, producing turbulent energy which is then dissipated into internal energy, heating up the cluster gas. However, at t > 12 Gyr, the simulation seems to globally reach some kind of stable state, comparable to what has been found in the driven turbulence simulations.

Before performing a quantitative analysis of the cluster properties in the fearless run, it is useful to visually inspect the generation and the spatial distribution of the turbulent energy in the ICM and around the cluster, in order to compare the simulations with the theoretical expectations of cluster mergers. Figures 7 and 8 present a time series of density and turbulent velocity slices, where several merger events in the cluster outskirts can be identified.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the density slice at redshift z = 0.15 (Figure 7(a)), we can see a filament extending from the lower left to the upper right corner of the figure; material is falling onto the cluster along this structure. On both sides, the inflow of relatively cold gas from the filament onto the ICM produces a moderate increase of turbulent energy (cf. Nagai & Kravtsov 2003). From the upper right side there is not only a smooth inflow of matter, but two small clumps are approaching the cluster. During the simulation, these two clumps merge with the main cluster (Figures 7(c)–8(c)) and one of them (on the left) is assimilated completely at redshift z = 0. A substructure approaches the cluster from the lower left corner along the filament (Figure 8(a)) and another one is visible only at z = 0.05 just to the right of the cluster core, when it crosses the slicing plane.

The merging process can be followed much more easily in terms of production of turbulent energy, visualized by the turbulent velocity $q = \sqrt{2 e_{\mathrm{t}}}$. In Figure 7(b) at z = 0.15, a marked peak of turbulent energy in the center of our cluster, resulting from a former massive merger, can be seen. The turbulent energy produced by this merger declines (Figures 7(d)–8(d)), and at z = 0 it is dissipated into internal energy nearly completely, as confirmed by our further analysis in Section 5.4.

The two approaching clumps described above continue to drive turbulence in the cluster. Thereby, the left clump can be identified in the turbulent velocity slice at z = 0.15 (Figure 7(b)) as a ring-like structure, showing that turbulence is not produced at the center of the infalling clump but at the front (behind a bow shock) and in the wake of the infalling material. The right clump only shows some turbulence production in its wake, which might be due to its smaller size and smaller velocity. On their way toward the main cluster and through its ICM, both clumps show a relevant production of turbulent energy (Figures 7(d)–8(b)). The considerable amount of turbulent energy can even be identified after the two clumps have merged with the main cluster (Figure 8(d)) and the left one is not easily visible in the corresponding density slice (Figure 8(c)). From this point of view, the distribution of turbulent energy traces the local merging history of a galaxy cluster until it is dissipated into internal energy completely.

The morphological evolution of the cluster gives a clear sense of the markedly local behavior of the production and dissipation of turbulence, which is confirmed to be an intermittent process in the ICM.

It is extremely difficult to apply an energy analysis similar to that performed in Section 3.1 to a galaxy cluster. Different than a periodic box, a galaxy cluster is an open system, with a growth over time of negative gravitational potential energy. Nevertheless, a comparison of the energy contributions of the two simulations at z = 0 is useful to understand the role of the SGS model.

The results are summarized in Table 1, where the energies in a sphere centered at the cluster center and with r = R_vir are reported. Different than elsewhere in this work, in Table 1 the energies are not specific, i.e., E_tot = ρe_tot, etc. This choice allows a better evaluation of the energy budget but it does not differ appreciably from the analysis of specific energies, since the baryon masses in the two runs agree within 1%.

The total energy remains basically unaltered in the two simulations, whereas the most important change is the decrease of E_kin in the fearless run. The missing kinetic energy is transferred mostly to E_t, which acts as a buffer between the resolved kinetic energy and E_int. The SGS model transfers energy from E_t to E_int either adiabatically (via the pressure dilatation term, Equation (15)) or irreversibly (via the dissipative term, Equation (16)). Turbulent dissipation is thus added to the dominant numerical dissipation, resulting in a moderate increase of E_int, though it is less relevant than the variation of E_kin.

In both runs, the kinetic energy contribution in the cluster E_kin is smaller than E_int. The mass-weighted average of the Mach number in the cluster is about 0.6, in agreement with the known fact that the flow in the ICM is, on average, mildly subsonic. In this regime, it is not surprising that the subgrid energy, E_t, is about 2 orders of magnitude smaller than E_int on the global level. The energy contribution from unresolved scales is globally negligible, though locally turbulence can play a more significant role, as will be discussed in Section 5.4.

Finally, we note that the ratio of the turbulent production term Σ (Equation (13)) to the turbulent dissipation term (Equation (16)) in the cluster core is

$$\begin{eqnarray} &&\frac{\Sigma }{\epsilon } = 0.93, \end{eqnarray} \tag{ 47 }$$

suggesting that the turbulent flow in the ICM is globally in a regime of near equilibrium of production and dissipation of turbulent energy.

The results of Section 5.3, referring to global features of the galaxy cluster, will be complemented by a local comparison in terms of radial profiles of selected physical quantities and an analysis of the cluster core in this section.

As a first interesting result of the comparison between radial profiles in the fearless simulation and in the standard adiabatic run, one can notice (Figure 9) that the temperature profile of the fearless run deviates slightly from that of the adiabatic run. This is especially apparent at the center, where T is larger in the cluster core for r ≲ 0.07 R_vir, with respect to the standard run. Consequently (Figure 10), the core in the fearless run is less dense, so that the ICM remains in hydrostatic equilibrium. The local energy budget in the cluster core is therefore modified by the SGS model.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Table 2, we explore this feature in detail, reporting the mass-weighted averages of selected variables in a sphere within 0.07 R_vir from the cluster center. Because of the adjustment of the cluster hydrostatic equilibrium, the mass enclosed in this sphere is significantly different in the two runs (it decreases by 10% in the fearless run), thus it is more convenient to present specific energies in Table 2.

Table 2. Mass-weighted Averages in a Sphere of r = 0.07 R_vir, Centered at the Cluster Core, at z = 0

Quantity	Adiabatic Run	fearless Run
Σ/	⋅⋅⋅	0.59
etot(1016 cm2 s-2)	1.3781	1.4189
eint(1016 cm2 s-2)	1.2529	1.3030
ekin(1015 cm2 s-2)	1.252	1.138
et(1013 cm2 s-2)	⋅⋅⋅	2.078
pres/(pres + ptherm)(%)	1.46	1.52
vrms(km s-1)	196	204

Download table as:  ASCIITypeset image

First, the low value of the Σ/ ratio indicates that, at z = 0, in the cluster core region the dissipation of turbulence is dominant with respect to its production. This confirms the result of the morphological analysis in Section 5.2 that turbulence is not produced locally in the core by mergers at z < 0.15, but that it decays in this region. The impact of this turbulent dissipation on the local energy budget of the cluster core can seen from the comparison of the energy contributions in Table 2. Similar to the global analysis in Table 1, there is a clear decrease of e_kin, transferred both to e_t and e_int. Both e_tot and e_int are higher in the fearless run, pointing to the existence of an energy flux from the resolved scales to the thermal reservoir through the turbulent buffer, leading to the increase of the internal energy. We interpret this additional energy contribution as caused by the turbulent dissipation introduced by the SGS model.

In the cluster core, the energy content at the SGSs is marginal. Apparently the relative contribution to the total energy is even smaller than in Table 1, but one should notice that, for consistency, in that table both E_kin and E_t are reported according to the original scale separation introduced by the AMR resolution, and without rescaling E_t as described in Section 5.1. In the cluster core, the refinement level is maximum, therefore the unresolved part of the turbulent cascade is relatively smaller than elsewhere, and so is e_t. In Table 2, we use the scaled definition of e_t; but in the core, it differs from the unscaled one only marginally, because in this region the resolution is l_min almost everywhere.

To further quantitatively appreciate the contribution of e_t to the energy budget, Figure 11 reports the profile of the turbulent pressure support p_t/(p_t + p_therm) in the cluster, where the turbulent pressure is defined as p_t = 1/3 ρq², and p_therm is the usual thermodynamical pressure. This ratio is also equal to the ratio of the corresponding energies (e_t/(e_t + e_int)). At the length scale of the effective spatial resolution of the simulations l_min = 7.8 kpc h⁻¹, the contribution of the turbulent pressure (or energy) is well below 1%, although it increases at larger central distances.

Zoom In Zoom Out Reset image size

In analogy with the turbulent pressure, we define a "resolved pressure"

$$\begin{eqnarray} &&p_{\mathrm{res}}=\frac{1}{3} \rho v_{\mathrm{rms}}^2, \end{eqnarray} \tag{ 48 }$$

where the rms velocity (Iapichino & Niemeyer 2008) is defined as

$$\begin{eqnarray} &&v_{\mathrm{rms}} = \frac{\sum _i m_i { \left(v_i - \langle v \rangle \right) }^2}{\sum _i m_i}. \end{eqnarray} \tag{ 49 }$$

Here, 〈v〉 is the mass-weighted average of the velocity in the analysis volume. This quantity essentially probes the contribution of turbulent motions at length scales of the order of 0.07 R_vir ∼ 90 kpc h⁻¹. As shown in Table 2, the pressure contribution at these length scales is at the percent level, and is slightly higher in the fearless simulation. Interestingly, the rms velocity in the fearless run is also somewhat larger than in the adiabatic case.

The changes in the temperature and density profiles are also reflected on the entropy which is defined, as is customary in astrophysics, as

$$\begin{eqnarray} &&K=\frac{T}{\rho ^{\gamma -1}} \end{eqnarray} \tag{ 50 }$$

with γ = 5/3. The entropy in the cluster core is higher in the fearless run as compared to the standard run (Figure 12). This result is consistent both with the locally increased dissipation of turbulent to internal energy provided by the SGS model and with the higher degree of mixing induced in the cluster core, shown by v_rms in Table 2.

Zoom In Zoom Out Reset image size

The radial distribution of turbulent energy is displayed in Figure 13 in terms of the turbulent velocity scaled to l_min. The turbulent velocity at this length scale is below 100 km s^-1. There is a pronounced peak at r = 0.6 R_vir, which is correlated with analogous trends in the turbulent pressure (Figure 11) and in the radial velocity profile (Figure 14). This structure is clearly linked with the most prominent merging clump shown in Figure 8(d), and analyzed below in Table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Mass-weighted Averages in a Volume of (512 × 768 × 1280) kpc h⁻¹, Containing a Subclump and its Wake, at z = 0.05

Quantity	Adiabatic Run	fearless Run
Σ/	⋅⋅⋅	1.13
etot(1016 cm2 s-2)	1.7447	1.5746
eint(1015 cm2 s-2)	5.4281	5.9607
ekin(1016 cm2 s-2)	1.2032	0.9786
et(1014 cm2 s-2)	⋅⋅⋅	2.290
pt/(pt + ptherm)(%)	⋅⋅⋅	3.70

Download table as:  ASCIITypeset image

There is an appealing similarity between the intervals of radii where q and the pressure ratio are larger (r < 0.1 R_vir and 0.4 R_vir < r < 0.8 R_vir), and the corresponding intervals where the temperature and entropy of the fearless run are slightly larger than those computed for the adiabatic run. The opposite trend occurs in the interval in between, where v_turb is comparatively smaller. The effects are very small, but suggest that the SGS model plays the same role in the ICM that was shown above for the cluster core, and in Section 5.3 for the global quantities. In case of radial profiles, the spherical averaging combined with the intermittent behavior of turbulence tends to mask the turbulent effects. This can be better understood with a comparison of the values of q in Figure 13 and in the right-hand panels of Figures 7 and 8: the peak values in the slices are much larger than the spherical averages in the profile.

The idea that, locally, the turbulence and its modeling can play a sizeable role is further corroborated by the data in Table 3, reporting the analysis at z = 0.05 of a small volume (512 × 768 × 1280 kpc h⁻¹) that contain one of the clumps presented in Section 5.2 and its wake (cf. Figure 8(a)). The morphology of this accreting subcluster in the fearless run is not substantially different from the adiabatic one. The energy content, however, is rather different from that in the cluster core: in the region under consideration e_kin is dominant with respect to e_int. The importance of the turbulence, injected by the hydrodynamic instabilities in the wake of the moving clump, is testified by the large ratio Σ/ and by the turbulent pressure support, which is at the level of some percent, about 1 order of magnitude larger than the spherical averages in Figure 11. Despite of the slightly smaller average of e_tot in the fearless run with respect to the adiabatic one, one can see an increase of e_int, mostly at the expenses of e_kin, resulting from the turbulent dissipation. The decrease of e_kin is rather large (19%) but can be partially ascribed to the difficulty of comparing energy budgets in such open volumes.

It is important to stress the deep difference between the turbulent velocity profile in Figure 13 and the profile of v_rms, defined by Equation (49) (see also Norman & Bryan 1999; Iapichino & Niemeyer 2008). In the former case, the mass-weighted average of a local quantity (i.e., defined in every cell) is computed for each spherical shell, whereas in the latter case v_rms is interpreted shell-wise as the standard deviation with respect to the average 〈v〉. Clearly, the latter definition does not retain any information related to a length scale, and can be interpreted as turbulent velocity only in a loose sense. From this point of view, the turbulent velocity provided by the SGS model is a more powerful probe of the features of a turbulent flow. On the other hand, spherically averaged velocity dispersions (and the derived turbulent pressure) are meaningful in comparison with observations, for example, in the procedure for estimating the cluster mass (cf. Rasia et al. 2006; Piffaretti & Valdarnini 2008). According to this different definition, the spherically averaged turbulent pressure of the simulated cluster (in a run similar to the adiabatic one presented here) is reported in Iapichino & Niemeyer (2008). It reaches values around 10%, in agreement with the values found recently in simulations by Lau et al. (2009) for relaxed clusters. The turbulent pressure is somewhat smaller in the fearless run because of its slightly reduced content in kinetic energy, but the difference is small, and is not expected to significantly affect the estimates of the cluster mass.