Adaptive mesh refinement simulations of a galaxy cluster merger - I. Resolving and modelling the turbulent flow in the cluster outskirts

TL;DR

This paper presents adaptive mesh refinement simulations of a galaxy cluster merger, focusing on resolving turbulent flows in the outskirts using a novel AMR strategy based on vorticity, including a subgrid turbulence model.

Contribution

The study introduces an AMR refinement criterion based on vorticity variability, improving turbulence resolution in cluster outskirts over previous methods.

Findings

The turbulence volume-filling factor in outskirts is about 60% at low redshift.

The AMR strategy effectively resolves shocks and turbulence in the cluster outskirts.

Implications for radio relic modeling and cosmic-ray interactions are discussed.

Abstract

The outskirts of galaxy clusters are characterised by the interplay of gas accretion and dynamical evolution involving turbulence, shocks, magnetic fields and diffuse radio emission. The density and velocity structure of the gas in the outskirts provide an effective pressure support and affect all processes listed above. Therefore it is important to resolve and properly model the turbulent flow in these mildly overdense and relatively large cluster regions; this is a challenging task for hydrodynamical codes. In this work, grid-based simulations of a galaxy cluster are presented. The simulations are performed using adaptive mesh refinement (AMR) based on the regional variability of vorticity, and they include a subgrid scale model (SGS) for unresolved turbulence. The implemented AMR strategy is more effective in resolving the turbulent flow in the cluster outskirts than any previously used criterion based on overdensity. We study a cluster undergoing a major merger, which drives turbulence in the medium. The merger dominates the cluster energy budget out to a few virial radii from the centre. In these regions the shocked intra-cluster medium is resolved and the SGS turbulence is modelled, and compared with diagnostics on larger length scale. The volume-filling factor of the flow with large vorticity is about 60% at low redshift in the cluster outskirts, and thus smaller than in the cluster core. In the framework of modelling radio relics, this point suggests that upstream flow inhomogeneities might affect pre-existing cosmic-ray population and magnetic fields, and the resulting radio emission.

Tables

Table 1. Table 1: Summary of the four simulations analysed in the present work. All runs follow the evolution of the same realisation of a cosmological volume and differ only for the AMR criterion from z = 2 𝑧 2 z=2 . The AMR overdensity factor f i subscript 𝑓 i f_{\text{i}} , both for DM and baryons (third column), used in the runs is indicated for clarity also by the number in the simulation name. The turbulence SGS model is used in all simulations.

Simulation	$f_{\text{i}}$	Cells [$\times {10}^6$]	$f_{\text{V}}$ [$\times {10}^{-7}$]	AMR
		with $l\ge 3$	($l=8$)	criteria
$OD4$+	4	13.35	148.9	OD and
				vorticity
$OD2$	2	8.34	28.2	OD only
$OD4$	4	3.69	6.23	OD only
$OD8$	8	1.76	0.76	OD only

Table 2. Table 2: Velocities computed from the mass-weighted energy data in Figure 6 , for the run O ​ D ​ 4 𝑂 𝐷 4 OD4 +. The reported velocities are v 𝑣 v (computed from e kin subscript 𝑒 kin e_{\text{kin}} ) and q 𝑞 q (from e t subscript 𝑒 t e_{\text{t}} ). The velocities are reported at z = 0 𝑧 0 z=0 and for the redshift of maximum e t subscript 𝑒 t e_{\text{t}} in the core/outskirts.

	$v$	$q$
	$[$km s${}_{}{}^{-1}]$	$[$km s${}_{}{}^{-1}]$
core
$z=0$	658	39.7
$z=0.350$	1151	73.1
outskirts
$z=0$	990	36.3
$z=0.275$	1283	82.2

Table 3. Table 3: Velocities computed from the mass-weighted energy data in Figure 6 , for the run O ​ D ​ 4 𝑂 𝐷 4 OD4 +, distinguishing between "turbulent” gas (second and third column) or ”non-turbulent” gas (fourth and fifth column). The reported velocities are v 𝑣 v (computed from e kin subscript 𝑒 kin e_{\text{kin}} ) and q 𝑞 q (from e t subscript 𝑒 t e_{\text{t}} ). The velocities are reported in the outskirts at z = 0 𝑧 0 z=0 and 0.275 0.275 0.275 .

	$v$ (turb)	$q$ (turb)	$v$ (no turb)	$q$ (no turb)
	$[$km s${}_{}{}^{-1}]$	$[$km s${}_{}{}^{-1}]$	$[$km s${}_{}{}^{-1}]$	$[$km s${}_{}{}^{-1}]$
outskirts
$z=0$	1030	40.0	845	19.7
$z=0.275$	1307	88.0	1142	39.4

Equations

$${\rho }_{\text{i}}>f_{\text{i}}{\rho }_0{\Omega }_{\text{i}}{N}^l,$$

$$Q(\text{mn@boldsymbol}x,t)\ge ⟨Q{⟩}_i(t)+\alpha {\lambda }_i(t),$$

$$\text{mn@boldsymbol}\omega =\nabla \times \text{mn@boldsymbol}v.$$

$$⟨\rho f⟩=⟨\rho ⟩\hat{f}\Rightarrow \hat{f}=\frac{⟨\rho f⟩}{⟨\rho ⟩}.$$

$${\hat{e}}_{\text{kin}}=\frac{1}{2}\widehat{v_i}\widehat{v_i}+\frac{1}{2}\hat{\tau }(v_i,v_j)/⟨\rho ⟩,$$

$$e_{\text{t}}=\frac{1}{2}{q}^2=\frac{1}{2}\hat{\tau }(v_i,v_i)/⟨\rho ⟩$$

$$p_{\text{t}}=\frac{2}{3}⟨\rho ⟩e_{\text{t}}.$$

$$\frac{\partial }{\partial t}⟨\rho ⟩e_{\text{t}}+\frac{\partial }{\partial r_j}{\hat{v}}_j⟨\rho ⟩e_{\text{t}}=\mathcal{D}+\Sigma +\Gamma -⟨\rho ⟩(\lambda +ϵ),$$

$$\frac{M(<R_{\Delta })}{4/3\pi R_{\Delta }^3}=\Delta {\rho }_c(z),$$

$$v/q\sim R_{\text{l}}^n$$

$$\mathcal{E}=\frac{1}{2}{\omega }^2.$$

$${\omega }_i>N/t_{\text{age}}(z),$$

$$v_{\text{rms}}(r)=\sqrt{\frac{{\sum }_i{m}_i{\sum }_j(v_{i,j}-⟨v_j(r)⟩{)}^2}{{\sum }_i{m}_i}},$$

$$v_1=v_2{\left(\frac{l_{\Delta ,1}}{l_{\Delta ,2}}\right)}^{1/3},$$

$$\frac{e_{\text{t,1}}}{e_{\text{t,2}}}=\frac{v_1^2}{v_2^2}\sim {\left(\frac{l_{\Delta ,1}}{l_{\Delta ,2}}\right)}^{2/3}.$$

$$\begin{array}{r}{e}_{\text{t,scaled}}\equiv e_{\text{t},l_{\text{max}}}=e_{\text{t}}{\left(\frac{l_{\Delta ,\text{max}}}{l_{\Delta }}\right)}^{2/3}=\\ =e_{\text{t}}{\left(\frac{2^l}{2^{l_{\text{max}}}}\right)}^{2/3}=e_{\text{t}}\times 2^{\frac{2}{3}(l-l_{\text{max}})}\end{array}$$

$${\mathcal{M}}_{\text{turb}}=\frac{v}{c_{\text{s}}}=\frac{v}{\sqrt{\gamma (\gamma -1)e_{\text{int}}}},$$

Full text

††pubyear: 2017††pagerange: Adaptive mesh refinement simulations of a galaxy cluster merger – I. Resolving and modelling the turbulent flow in the cluster outskirts–References

1 Introduction

The evolution of the large-scale structure of the universe proceeds through a sequence of merger events, by which collapsed objects contribute to form larger entities (e.g. Ostriker 1993; White et al. 1993). At the top of this hierarchy, clusters of galaxies reach a mass up to a few times $10^{15}$ M*☉*. The baryonic gas falls in the corresponding potential wells and is shock-heated to temperatures up to about $10^{8}$ K in the intra-cluster medium (ICM). A fraction of the gas tracks the filaments interconnecting and surrounding the collapsed halos, and can be found in a phase of lower density and temperature than the ICM, called the warm-hot intergalactic medium (WHIM).

Since the hot ICM gas emits in the X-ray wavebands, mainly by thermal bremsstrahlung emission, and the X-ray emissivity scales with $n_{\text{e}}^{-2}$, where $n_{\text{e}}$ is the electron number density, the X-ray emitting gas in the outskirts is very hard to detect. It is therefore not surprising that observations of the cluster periphery have received attention only in recent years (see the review by Reiprich et al. 2013).

The observations of the cluster outskirts have raised important questions about the physical conditions of the gas at those locations. In particular, the most striking discrepancy with analytical predictions (Tozzi & Norman, 2001) and simulations (e.g., Voit et al. 2005) is in the entropy profiles of the outskirts, flattening to values smaller than the theoretical expectations. In order to explain this feature, gas clumping has been invoked (Simionescu et al., 2011; Morandi et al., 2013; Ichinohe et al., 2015; Walker et al., 2013; Vazza et al., 2013; Roncarelli et al., 2013; Zhuravleva et al., 2013), as well as electron-ion equilibration timescales (see recently Avestruz et al. 2015, and references therein). An alternative possibility is that part of the pressure support in the cluster periphery comes from turbulent gas flows. This can be easily understood, since the outskirts are subject to accretion flows delivering material from the surrounding regions as part of cosmic structure formation, and had no time to settle in hydrostatic equilibrium, unlike the cluster core. This scenario is supported by many numerical simulations (Lau et al., 2009; Vazza et al., 2009, 2011b; Parrish et al., 2012; Shi et al., 2015). Accretion has been proposed to drive turbulence also on scales of galaxies or molecular clouds within these galaxies (see, e.g., the discussion by Klessen & Hennebelle 2010).

Probing the physics of the cluster peripheries with realistic numerical simulations has become mandatory for using galaxy clusters as tools in cosmology. Among the points touched above, the level of turbulence is of special importance because the presence of non-thermal pressure will result in underestimating the cluster mass, if hydrostatic equilibrium in the ICM is assumed (Rasia et al., 2004). The first direct detection of turbulence-induced X-ray line broadening in the ICM, performed on the Perseus cluster by the Hitomi satellite (Hitomi Collaboration et al., 2016), indicates that the correction due to turbulence pressure is not of leading order in the core regions of that relaxed object. The unfortunate loss of the satellite leaves the problem of direct detection of turbulence in other clusters with different dynamical histories open, most likely until one of the next upcoming missions like Athena111Website: http://www.the-athena-x-ray-observatory.eu/ is launched.

The important role of the turbulent kinetic energy component in the cluster outer regions has intriguing consequences for the study of non-thermal processes. Both MHD simulations (e.g., Dolag et al. 2002; Dubois & Teyssier 2008; Ryu et al. 2008; Dolag & Stasyszyn 2009; Xu et al. 2010; Bonafede et al. 2011; Xu et al. 2011; Miniati & Beresnyak 2015; Beresnyak & Miniati 2016) and analytical estimates222Also noteworthy in this framework is the development of laboratory experiments for the study of turbulence and turbulent magnetic amplification in astrophysical plasmas (Meinecke et al., 2014, 2015). (Iapichino & Brüggen 2012; see also Schober et al. 2013; Schleicher et al. 2013) indicate that magnetic fields, amplified by the turbulent dynamo, can be in the range between 0.1 and 1 $\mu$G. Indeed, observational analyses from different techniques (Clarke et al., 2001; Govoni et al., 2006; Bonafede et al., 2010) support this idea. Also, connected with shocks and turbulence is the acceleration of cosmic ray (CR) electrons (Brunetti & Jones, 2014). The presence of this particle component and of magnetic fields is indicated by the observations of diffuse emission in the radio wavelengths in a growing number of cluster peripheries; these objects are known as radio relics (Ferrari et al. 2008, for a review). For some relics, X-ray observations have detected the presence of a shock corresponding to the location of the radio emission (Giacintucci et al., 2008; Finoguenov et al., 2010; Macario et al., 2011; Akamatsu et al., 2012; Akamatsu & Kawahara, 2013; Botteon et al., 2016; Basu et al., 2016). Propagating shocks in the ICM have also an impact on the star formation rate of cluster galaxies (Stroe et al., 2014), which is in agreement with the dependence of star formation on Mach number in molecular clouds, studied by Federrath & Klessen (2012) and Renaud et al. (2012); Renaud et al. (2014).

A leading model for relic emission implies the acceleration of CR electrons at the shock (diffusive shock acceleration, henceforth DSA), and the subsequent synchrotron emission in the downstream magnetic field. The details of the acceleration and of the dynamics of CRs in relics go beyond the scope of the present study and are extensively referred elsewhere (Brüggen et al., 2012; Kang et al., 2012; Brunetti & Jones, 2014). Here it is sufficient to recall that there has been tension between the DSA mechanism as formation model for the relics and observational data. For example, the acceleration efficiency of cluster merger shocks (with a Mach number of a few) has been questioned (e.g., Kang & Ryu 2011). An increasing interest (Pinzke et al., 2013) is attracted by acceleration models invoking a pre-existing CR population, either produced by earlier stirring events or ejected by active galaxies (Shimwell et al., 2015). Turbulent reacceleration in relics has been proposed by Fujita et al. (2015). In principle, the propagation of a shock can thus expose the upstream conditions of the ICM, as far as CR population is concerned. Furthermore, we notice that inhomogeneities of the upstream magnetic field can mimic the same effect.

Since shocks and turbulence in the ICM are related both to the acceleration of the CRs and to the amplification of the magnetic fields, we take a first step in the study of the problem and investigate the properties of the flow in the cluster outskirts purely from the viewpoint of hydrodynamics, putting aside the MHD treatment and the acceleration model. This requires the choice of a numerical scheme able to properly model and resolve the ICM flow, and in particular in the outer cluster regions.

Dealing with the cluster outskirts and the surrounding filaments is extremely demanding for numerical simulations. Smoothed Particle Hydrodynamics is self-adaptive in simulations of collapsed objects like galaxy clusters, which is convenient to spatially resolve the dense central regions, but can be problematic in the mildly overdense surroundings (e.g. Vazza et al. 2011a). On the other hand, Adaptive Mesh Refinement (AMR) simulations can in principle allow to refine on any arbitrary variable, resulting in a matchless advantage in problems involving turbulent flows. Resolving turbulent flows has been motivating the AMR simulations performed by using refinement criteria based on gas velocity jumps (Vazza et al., 2009) or on regional variability of vorticity modulus and compression rate (Iapichino & Niemeyer 2008; Schmidt et al. 2009, 2015; see Section 2 for a detailed description). In recent years, the increased availability of computational power has revived the static grid approach (with some technical variants), in simulations able to provide an inertial range (the interval of length scales between the turbulence driving scale and the dissipation scale) of at least one decade in length scale, suitable to resolve the turbulent flow in the ICM and around the cluster (Vazza et al., 2014a; Vazza et al., 2014b; Miniati, 2014, 2015; Vazza et al., 2017).

Previous numerical studies of turbulence in galaxy clusters (Iapichino et al., 2011) point to an important difference between the ICM and the WHIM. In those simulations, the time evolution of the turbulence energy in the ICM shows a peak at redshift between $z\simeq 1$ and $0.6$, while in the WHIM the turbulence energy grows steadily until the current epoch. The explanation for these different evolutions is provided by the turbulence injection mechanisms in the two gas phases: turbulent flow is induced mainly by cluster major mergers in the ICM, and by accretion through structure formation shocks in the WHIM. For the latter, and in particular for the cluster outer regions, a similar result has been developed analytically by Cavaliere et al. (2011).

The simulation discussed by Iapichino et al. (2011) makes use of a subgrid scale (henceforth SGS) model for the computation of the energy content on unresolved length scales (Schmidt et al., 2006a; Maier et al., 2009). This tool, to be discussed in Section 2, is very suitable in simulations of turbulent flows in contracting objects with a large dynamical range. On the other hand, a limit of that work is a modest effective spatial resolution ($48.8$ kpc $h^{-1}$), necessary to simulate a relatively large cosmological volume with a feasible dynamical range. Zoomed simulations of single clusters, with a combination of static grids and AMR allowing a much finer spatial resolution, are thus needed for a more complete and complementary understanding of the cluster energy budget and its time evolution.

The current work and its companion study (henceforth Paper II) aim to focus on the cluster outskirts and further elaborate on the properties of the turbulence flow in those regions, as opposed to the cluster core, during a major merger. We perform a number of zoomed, grid-based cosmological simulations of the evolution of a massive cluster, including turbulence SGS model (Schmidt et al. 2006a, with the addition of the improvements by Schmidt & Federrath 2011) and AMR criteria suitable for refining turbulent flows (Section 2). By using these tools and analysing the flow properties, we will study the following relevant questions for the ICM physics and for the turbulence driving mechanisms during cosmological structure formation:

Is the turbulence in the cluster outer regions controlled by accretion of pristine gas, as Iapichino et al. (2011) stated for the WHIM, or is it rather affected by merger events occurring in the ICM? 2. 2.

What are, during the cluster evolution, the dominant stirring modes (solenoidal or compressional)? 3. 3.

How intense and volume-filling is the turbulent flow in the cluster outskirts?

This first paper focuses mainly on the numerical methods, namely on the effectiveness of the AMR strategy based on regional variability of vorticity in refining the cluster outer region, when compared with other methods based on overdensity (Section 2.1). Furthermore, we will show that the employed turbulence SGS model provides a quantity (the subgrid turbulent kinetic energy, defined in Section 2.2) which is useful in the characterisation of the turbulent flow during the cluster evolution and, being computed locally (cell-by-cell), it is more easily accessible than global (volume-averaged) diagnostics. It will be discussed how these observables of turbulence on different length scales compare to each other during a merger event. We will limit the study of the general properties of turbulence in the ICM to the minimum which is needed to interpret the results of the SGS model, well aware that results on large scales are available in literature and cited above. In Paper II, more devoted to the physical interpretation of simulation results, we further analyse the velocity and density field both of the innermost and the outer cluster regions, to understand how the turbulent flow evolves during a major merger, and which stirring modes are dominant.

The structure of the current paper is the following: in Section 2 we describe our numerical tools, in particular the AMR criteria and the turbulence SGS model. Our results are presented in Section 3, where we focus in particular on comparisons of different refinement strategies and diagnostics of the turbulent flow, both on large and subgrid scales. The results are discussed in Section 4, and our conclusions are summarised in Section 5.

2 Numerical methods

The numerical simulations presented here have been performed using the enzo code (O’Shea et al., 2005a; Bryan et al., 2014), in a customised version based on the public release v. 2.1. enzo is a hybrid, $N$-body plus hydrodynamical grid-based code featuring block-structured AMR (Berger & Colella, 1989). The $N$-body section of the code is based on a particle-mesh scheme (Hockney & Eastwood, 1988), and the hydrodynamics on the Piecewise Parabolic Method (PPM; Colella & Woodward 1984).

The simulated box has a size of $256$ Mpc $h^{-1}$ on a side, and is evolved from redshift $z=60$ to $z=0$. The initial conditions have been produced according to the Eisenstein & Hu (1999) transfer function. We adopted the cosmological parameters coming from the Planck 2013 results (Planck Collaboration et al., 2014), and consequently set $\Omega_{\Lambda}=0.693$, $\Omega_{\text{m}}=0.307$, $\Omega_{\text{b}}=0.0482$, $h=0.678$, $\sigma_{8}=0.826$, and $n=0.9608$.

Besides the turbulence SGS model (Section 2.2), we did not include any additional physics to our runs, as we are interested primarily in cluster regions and length scales where extensions beyond the "adiabatic" description of the fluid play a minor role (for some investigation, see Eckert et al. 2012). We model the gas as a fluid described by an ideal equation of state, with $\gamma=5/3$.

The simulation setup has been optimised for resolving the formation of a single galaxy cluster in its cosmological environment. The simulation volume is resolved in $128^{3}$ grid cells and $128^{3}$ $N$-body particles modelling the dark matter (DM). Inside this root grid (whose resolution is also identified by the refinement level $l=0$) two additional static grids are nested around the box centre ($l=1$ and 2). Each of these nested grids is resolved in $128^{3}$ grid cells and $128^{3}$ $N$-body particles; they have a size of $128$ Mpc $h^{-1}$ ($l=1$) and $64$ Mpc $h^{-1}$ ($l=2$), respectively, with the $N$-body particle mass at $l=2$ being $8.5\times 10^{9}$ M*☉* $h^{-1}$. In an innermost volume of $38.4$ Mpc $h^{-1}$ on a side AMR is activated, according to the refinement strategies discussed below (Section 2.1). The size of this region is chosen in such a way to be, at $z=2$, about four times larger than the Lagrangian volume of the $N$-body particles that, at $z=0$, are part of the cluster. The finest allowed AMR level is $l=8$, corresponding to an effective spatial resolution of $7.8$ kpc $h^{-1}$.

2.1 Criteria for adaptive mesh refinement

In mesh-based codes AMR is a crucial tool for gaining the required dynamic range to model strongly clumped systems like galaxy clusters. However, by its very nature, AMR has issues in dealing with high refinement in extended regions of the computational volume. The choice of the refinement strategy is of great importance, also in order to avoid over-refinement on uninteresting locations: this would spoil the advantage of using AMR with respect to other alternatives (static grids, both monolithic or nested), and moreover would make the computation infeasible on modern computer architectures, where the amount of available memory per core is a significant constraint.

A sensible and widely used AMR choice in simulations of cosmological structure formation is refinement on overdensity, referring both to the DM and the baryon component. This criterion is already implemented in enzo, and a cell is flagged for refinement if the local density $\rho_{\text{i}}$, where "i" can indicate either baryons or DM, fulfils the following criterion:

[TABLE]

where $\rho_{0}=3H_{0}^{2}/8\pi G$ is the critical density, $H_{0}$ is the Hubble parameter at the present epoch, $\Omega_{\text{i}}$ are the cosmological density parameters for either baryons or DM, and the refinement factor is $N=2$. In this work the parameters $f_{\text{i}}$ for overdensity are set such that $f_{\text{b}}=f_{\text{DM}}$.

The choice of the ideal overdensity parameter depends on the kind of problem that one wants to address. The work of O’Shea et al. (2005b) has shown that the DM refinement is critical for properly resolving substructures. Our choices of $f_{\text{i}}$ are between $2$ and $8$; we present them in detail in Table LABEL:tab1sims.

Among the different choices for refinement criteria for intermittent turbulent flows (Kritsuk et al., 2006; Vazza et al., 2009), the approach chosen in this work is based on the regional variability of so-called structural invariants of the flow, i.e. on variables related to spatial derivatives of the velocity field. This method has been introduced by Schmidt et al. (2009) and applied by Iapichino et al. (2008) and Schmidt et al. (2015) in idealised subcluster simulations, and in full cosmological simulations by Iapichino & Niemeyer (2008), Paul et al. (2011), and Schmidt et al. (2014). Here for the first time it is optimised for the study of cluster outskirts.

According to this criterion, a cell is flagged for refinement if the local value of the variable $Q(\mn@boldsymbol{x},t)$ is

[TABLE]

where $\lambda_{i}$ is the maximum between the average $\langle Q\rangle$ and the standard deviation of $Q$ computed on a grid patch $i$, and $\alpha$ is a tunable parameter. The rationale of this selection procedure is somewhat similar to the method by Zhuravleva et al. (2012) for selecting density inhomogeneities, although in the present work we decide to use a different variable. Instead of density, we refine on the vorticity modulus, $Q=\omega^{2}$. Up to a factor of 2, this quantity is also known as enstrophy (e.g., Vazza et al. 2017). The vorticity is defined as the curl of the gas velocity:

[TABLE]

As for the refinement parameter, we empirically set the value of the threshold factor $\alpha=6.5$ like in Iapichino & Niemeyer (2008), for an optimal compromise between computational feasibility and refinement efficiency.

Table LABEL:tab1sims presents an overview on the performed runs. All simulations were run from $z=60$ to $z=2$ using the refinement criterion based on overdensity with thresholds $f_{\text{i}}=4$ (that is $OD4$), and then evolved further with the criterion specified in the Table. As one can see, the number of AMR cells at $z=0$ (third column) spans almost one order of magnitude for the simulation sample. The fourth column reports $f_{\text{V}}$, defined as the fraction of the computational volume which is resolved at the highest refinement level $l=8$ (normalised to the volume where AMR is allowed). This quantity (anticipated here from Fig. 3) has also a wide range of values within our sample, and together with the previous variable can be used to estimate the computational cost of the simulations presented here.

2.2 Modelling of unresolved turbulence

Our hydrodynamical code makes use of a SGS model for unresolved turbulence. This tool has been described first by Schmidt et al. (2006a) as an essential part of the flame propagation model in simulations of Type Ia supernovae (Schmidt et al., 2006b). In simulations of the evolution of the cosmic large-scale structure it has been used among others by Maier et al. (2009), Iapichino et al. (2011); Iapichino et al. (2013), Braun et al. (2014), and Schmidt et al. (2014). We notice that the adoption of turbulence SGS models has become more frequent in computational astrophysics in recent years (see e.g. Scannapieco & Brüggen 2008; Close et al. 2013), with pioneering work in SPH (Shen et al., 2010).

While the details of the model are contained in many of the works cited above (and we additionally refer the reader to Schmidt 2015 for a recent review), we recall here only the concepts that are crucial for a thorough understanding of what follows in this paper.

The turbulence SGS model is based on the decomposition of the fluid equations governing a physical system into a large-scale and a small-scale (unresolved) part, as described by the Germano (1992) formalism. This decomposition, when kept implicit in mesh-based simulations, consists merely in neglecting the effect of unresolved scales, which is equivalent to claiming that the numerical viscosity can mimic its effect. A more elaborate approach, however, is to explicitly select a filtering length scale and then to decompose any density-weighted variable of interest $f$ into a smoothed (i.e. large-scale) part $\langle f\rangle$ and a fluctuating part $f^{\prime}$, varying only on length scales larger than the filter (Favre, 1969). In this way, a filtered variable $\hat{f}$ is defined by the relation:

[TABLE]

When applied to the equation of fluid dynamics for a compressible, viscous, self-gravitating fluid, this formalism leads to new terms in the momentum and energy equations (see Schmidt 2015 for details), describing the interaction between resolved and unresolved scales. In particular, the filtered kinetic energy is expressed by

[TABLE]

where the first contribution on the right-hand side is the resolved kinetic energy. The second term contains the second-order moment of the velocity field $\hat{\tau}(v_{i},v_{j})=\langle\rho v_{i}v_{j}\rangle-\langle\rho\rangle\hat{v_{i}}\hat{v_{j}}$, and the trace of $\hat{\tau}(v_{i},v_{j})/\langle\rho\rangle$ can be interpreted, like in Germano (1992), as the square of the SGS turbulence velocity $q$, so that

[TABLE]

defines the SGS turbulence energy. Following from the definition of $e_{\text{t}}$, and given that the trace of $\hat{\tau}(v_{i},v_{j})$ is added to the thermal pressure in the filtered momentum equation (Schmidt, 2015), one can define the SGS turbulence pressure, associated with the unresolved turbulent velocity fluctuations on length scales smaller than the filter as

[TABLE]

The SGS turbulence energy is governed by the additional fluid equation,

[TABLE]

where the terms $\mathcal{D}$, $\Sigma$, $\Gamma$, $\lambda$ and $\epsilon$ at the right-hand side are source and transport terms related to the scale separation. We refer the reader to Maier et al. (2009) for their definition and description. These terms are expressed by empirical closures, and their formulation represents the SGS model. The closures used in this work are the same as in Maier et al. (2009) and Iapichino et al. (2011), with the exception of the formulation for the SGS turbulence stress tensor (contained in the definition of $\Sigma$, according to the so-called eddy viscosity closure), which follow the improved prescription by Schmidt & Federrath (2011).

An essential feature of the turbulence SGS model in its implementation within an AMR code is the consistent bookkeeping of the energy components, in particular of the energy exchange between the resolved kinetic energy and the turbulence SGS energy due to a variation of the filtering length scale, at grid refinement or derefinement. Only at that step of the algorithm, one needs to resort to the additional assumption of Kolmogorov scaling of the turbulent energy (Kolmogorov, 1941), according to the procedure described in detail by Maier et al. (2009) and Iapichino et al. (2011). As Kolmogorov scaling is appropriate for incompressible turbulence, this choice is justified on the length scales of the order of several kiloparsec, near the spatial resolution of our simulations, whereas turbulence is nearly transonic up to slightly supersonic on injection length scales (see later in Section 4.1).

3 Results

3.1 Properties of the simulated cluster

With the use of preliminary, low-resolution simulations, the initial conditions of the runs have been spatially shifted (making use of the periodic boundary conditions) in such a way to form the most massive structure at the box centre, with the nested static grids and the AMR region enclosing this location. This cluster, at $z=0$ and in run $OD4$+, has a virial mass $M_{\text{vir}}=0.95\times 10^{15}$ M${}_{\sun}\,h^{-1}$ and a virial radius $R_{\text{vir}}=2.13$ Mpc $h^{-1}$; these values differ for the other runs at a level of only a few per cent. Throughout this work we will keep a definition of cluster virial radius by setting, like e.g. in Bryan & Norman (1998), the virial overdensity as $\Delta^{c}_{\text{vir}}(z=0)\approx 102$, where the overdensity $\Delta$ is defined as the ratio between density and $\rho_{c}(z)$, the critical density of the Universe at redshift $z$. For sake of comparison with previous studies, $R_{\text{vir}}\approx 1.36\,R_{200}$ (Reiprich et al., 2013), where the radius $R_{\Delta}$ at overdensity $\Delta$ is defined by the equation

[TABLE]

here $M(<R_{\Delta})$ is the cluster mass enclosed within a sphere of radius $R_{\Delta}$.

The analysis of the mass accretion history of this cluster, performed with the hop algorithm (Eisenstein & Hut, 1998), clearly shows the cluster undergoes a major merger between $z=0.515$ and $0.427$, with a mass ratio close to 1. The merger plane is almost perpendicular to the $x$-axis of the computational box, therefore in Figure 1 the merger is optimally shown by slices on the $yz$-plane. In this time sequence of gas density slices the two subclusters can be still distinguished at $z=0.427$.

Figure 2 shows a better view of the merger event through temperature slices. In those panels, the launched merger shock is located ahead of the high-temperature, shock-heated region, propagating outside. Around $z=0.210$, the shock crosses the virial radius, then propagates into the outskirts, where it remains visible until $z=0$. During its propagation inside the virial radius the merger shock has a Mach number initially ($z=0.35$) between $2.5$ and $4$, and at later times ($z=0.275$) between $4$ and $6$, the range of values depending on the different shock regions. Interestingly, a visual inspection of Figure 1 shows that, after the major merger, the cluster accretes smaller clumps (approaching from below until $z=0.275$). These multiple minor mergers lead to additional mass growth of $1.4\times 10^{14}$ M${}_{\sun}\,h^{-1}$ between $z=0.350$ and $0.275$. This process launches a second shock, weaker than the previous one and hardly visible inside the virial radius in Figure 2 at $z=0.15$, with a Mach number $\mathcal{M}\lesssim 3$. Such a complex merger scenario is ideal for studying the stirring of turbulence in the ICM and somewhat complementary to the works of Iapichino & Niemeyer (2008) and Maier et al. (2009), based on the analysis of a relaxed galaxy cluster.

3.2 Comparison of mesh refinement strategies

The simulations presented in this work differ from each other only in their refinement strategy and resulting ability in resolving the formation of cosmological structure. A first insight on the different AMR effectiveness has been provided by the number of computational cells in Table LABEL:tab1sims. Besides the total number of grid cells at $z=0$, an interesting indicator of the refinement effectiveness is the volume fraction refined at an AMR level $l\geq n$ (Figure 3). This is normalised to the volume of the part of the computational domain were AMR is allowed (a cube with size of $38.4$ Mpc $h^{-1}$), so that this volume is equal to unity for $l=2$. One can see that the AMR criterion in run $OD8$ has the typical performance of simulations of cosmological large-scale structure, with every level occupying a volume approximately one order of magnitude smaller than its coarser one. Smaller overdensity thresholds result in refinement of larger volumes: the run $OD2$ refines a volume that is a factor of 37 larger than in run $OD8$ at the finest level $l=8$. The refinement on vorticity is effective especially on large AMR levels: the run $OD4$+ refines more than $OD2$ from $l>5$, and provides an additional factor of 5.3 more resolved volume at $l=8$ than the latter.

We emphasise that this strategy is especially successful in the refinement of the cluster outskirts. In order to demonstrate this, the cluster radial profiles of the mass-weighted AMR level at $z=0$ are shown in Figure 4. The profiles show how standard refinement criteria based on overdensity degrade their performance with distance from the cluster centre, and that even very permissive thresholds like in $OD2$ are of relatively little help outside $R_{\text{vir}}$. The refinement level in run $OD4$+ is smaller than that of $OD2$ at $r<0.6\ R_{\text{vir}}$, comparable for $0.6<\ r/R_{\text{vir}}<0.9$, and better for $r>0.9\ R_{\text{vir}}$. Not only the value of the AMR level is larger for the run $OD4$+, but we notice also that the slopes of the radial profiles for the other criteria are so low that only an unfeasibly small overdensity threshold (if any at all) would match the resolution of run $OD4$+ in the outskirts. At $r=2\ R_{\text{vir}}$ the refinement level of run $OD4$+ is larger by $\Delta l=1.4$ and $0.8$ compared to runs $OD8$ and $OD2$, respectively, corresponding to a better resolution by a factor of $2^{\Delta l}=2.64$ and $1.74$ for the two cases. In the best-resolved run, the mass-weighted AMR level within $r=2\,R_{\text{vir}}$ is better than $l=5$, corresponding to a spatial resolution of $62.5$ kpc $h^{-1}$.

3.3 Resolving turbulence stirring in the cluster outskirts

The main objective of AMR is to improve the resolution of the simulated physical system through a finer grid meshing. On the other hand, it is difficult to determine how effective an AMR strategy is in correctly capturing a turbulent flow, when simulating complex systems like galaxy clusters. Since the scope of the present work is to resolve turbulent flows, it is important that the proposed AMR strategy results in a good spatial coverage of the stirring agents, not only in the ICM but also in the cluster outskirts. The potential of the AMR based on regional variability of vorticity induced by subcluster motion has been shown by Iapichino et al. (2008) and Schmidt et al. (2015) in idealised simulations, and by Iapichino & Niemeyer (2008) in full cosmological simulations. Moreover, stirring driven by propagating merger shocks has been addressed by Paul et al. (2011). Rather than repeating these analyses, here we provide a different example, and focus on the resolution of turbulence driven by infalling cosmic filaments onto the ICM.

Cosmic filaments interconnect galaxy groups and clusters and contribute significantly to their mass accretion. Their role as turbulence drivers has been recently emphasised by Zinger et al. (2016), and put into relation with other properties of galaxy clusters like their classification as relaxed objects or the presence of a cool core. At the interface where filaments impinge on the ICM the gas is shocked, further contributing to thermalisation and flow stirring. Moreover, the interaction between propagating shocks and filaments drives hydrodynamical shearing instabilities, which have been related to radio-bright notches at the edges of some observed radio relics (Paul et al., 2011).

In Figure 5 one can see, in the density slices, several examples of filaments, like the network located from the lower right corner to the upper left, or the thin structure coming from the upper right corner. Because these filaments are only moderately overdense, their spatial resolution depends crucially on the employed refinement criterion. The AMR used in run $OD4$+ outperforms the other strategies, as clearly visualised in the AMR level slices on the left-hand side, not only in terms of the morphology of the filaments (density slices, second column of Figure 5) but also for the small-scale vorticity (third column of Figure 5). A filamentary pattern typical of a turbulent flow extends well beyond the virial radius (up to about $2\,R_{\text{vir}}$ is visualised in the slices, and analysed in the following). Idealised simulations of supersonic turbulence show very similar patterns of enhanced vorticity at the locations of shocks (see e.g., figure 2 in Federrath et al. 2010). We can therefore conclude that the refinement criterion based on vorticity is adequate for resolving turbulent driving.

In the last two columns of Figure 5, two different cell-based diagnostics of the turbulent flow are compared, namely the vorticity modulus of the flow and the SGS turbulence energy. The first one is a local indicator of spatial fluctuations of velocity, typical of turbulence, and is therefore a widely-used quantity for the study of the properties of flow, either in its basic definition or by isolating single terms of it (e.g. Miniati 2014; Vazza et al. 2017). A somewhat similar behaviour can be seen in the evolution of the SGS turbulence energy $e_{\text{t}}$, which is computed cell-wise by the turbulence SGS model employed in this work. These two quantities are indeed linked, as showed in Iapichino et al. (2011) by the correlation of the terms expressing the resolved and SGS turbulent pressure (their figure 9). The most important conceptual difference between $\omega$ and $e_{\text{t}}$ is the dependence of the latter on grid resolution: in its definition, the cell size acts as the filtering length scale introduced in Section 2.2. The key consequence of this definition is that, in multi-resolution and/or AMR grids, cells at different resolution convey different information (i.e referring to different cutoff scales) about the subgrid scale turbulence. The meaningful use of this variable involves therefore some additional step for correctly interpreting the results. This may seem cumbersome at first, but we will show (Section 3.4) that the analysis of $e_{\text{t}}$ provides similar results to the ones coming from other flow diagnostics, with the advantage of being computed locally on the grid without requiring data post-processing. The reason for it is that the properties of large-scale turbulence are imprinted onto the smallest scales through the turbulence cascade.

In considering the slices of $e_{\text{t}}$ in Figure 5, we notice that the SGS turbulence energy tracks well the large-scale structure surrounding the cluster, in particular some filaments down to within the virial radius. This happens because $e_{\text{t}}$ correlates with shocks, in particular with the merger shock and the external ones, surrounding the cluster and the filaments. It is expected to find a larger $e_{\text{t}}$ in the post-shock regions (Paul et al., 2011; Iapichino & Brüggen, 2012), because shocks inject vorticity through the baroclinic mechanism (e.g., Ryu et al. 2008; Vazza et al. 2017). Moreover, it has been demonstrated that the turbulence SGS model with the closure prescription by Schmidt & Federrath (2011), implemented in the simulation code used for this work, is reliable in computing the SGS energy in compressible flows, like the ones in the vicinity of shocks.

To complete the comparisons performed in this Section, the last column of Figure 5 presents an interesting consequence of the changing resolution level across the simulations, namely the variation in the SGS turbulence energy. Not surprisingly, in general $e_{\text{t}}$ is smaller in run $OD4$+, where the resolution is the highest, and therefore the unresolved part of the turbulent cascade is the smallest with respect to the other runs. For the same reason, the relatively large values in the outskirts of run $OD8$ mean that the flow in those regions is severely underresolved. These effects must be kept in mind for a correct interpretation of the results of the turbulence SGS model, as will be discussed in the following.

3.4 Evolution of the kinetic energy on resolved and subgrid scales

In the previous Section the SGS turbulence energy $e_{\text{t}}$ has been compared with another indicator of the flow, the vorticity modulus. Here we compare $e_{\text{t}}$ and other diagnostics of the gas energy content, by taking averages over control volumes, and studying the different ways they evolve during a major merger.

The analysis presented in this Section makes use of a somewhat arbitrary distinction between the core region of the galaxy cluster and its outskirts. This definition is based on the cluster centre (defined as the location of highest DM density) and virial radius at every data output: the cluster "core" is enclosed in a sphere with radius of $0.5\,R_{\text{vir}}$ around the cluster centre, while the outskirts are defined by the spherical shell ranging from $0.5$ to $2\,R_{\text{vir}}$. This definition can be problematic during the merger phases where spherical symmetry is not a good approximation. However, tests with different values of boundary radii showed that the results are not significantly affected.

As for the quantities included in the comparison, besides $e_{\text{t}}$ we will make use of the specific internal energy $e_{\text{int}}$, and of the kinetic energy $e_{\text{kin}}$. The internal energy is a useful comparison term in our analysis, because in virialised objects like galaxy clusters there is a strong link between development of bulk flow and turbulence and gas thermalisation during mergers. In getting a definition of kinetic energy, one should keep in mind that the velocity of the centre of mass of a cluster in a cosmological box can be of several hundred km s*-1*, which is of the same order of magnitude as the typical flow velocities. To get a cleaner indicator, the centre of mass velocity (computed by the hop halo finder in the yt toolkit) is subtracted from the velocity components, before computing the specific kinetic energy.

We consider the evolution of energy in the core (Figure 6$a$) first. From the time evolution of $e_{\text{int}}$ one can recognise the merger sequence described in Section 3.1. After early fluctuations due to the initial cluster growth, the internal energy in the cluster core shows in fact a late double peak, related to the propagation of merger shocks released by the major and second merger, respectively. The double peak is also apparent in the time evolution of the SGS turbulence energy. The steep decrease of $e_{\text{t}}$ right after $z=2$ is purely numerical, and caused by the activation of the AMR based on vorticity at that time. Different from the previous two quantities, the kinetic energy $e_{\text{kin}}$ shows an increase that starts earlier than the double peak (already at $z=0.5$). We will infer later that this is caused, at least in its first part, by the flow of filament gas and subclumps at merger stage, rather than by the merger shock propagation itself.

The effect of the major merger is then visible in the outskirts (Figure 6$b$) with some delay, due to the propagation of the merger shock outwards; the second minor merger has nearly no impact outside the cluster core. Also here, the peaks of $e_{\text{int}}$ and $e_{\text{t}}$ at $z=0.275$ coincide, but the peak of the kinetic energy is broader. From this analysis we can draw the conclusion that the kinetic energy is sensitive both to the stirring caused by the cluster merger and by the subsequent merger shock propagation. The SGS turbulence energy, being the quantity most sensitive to smallest length scales, seems affected mostly by turbulence injected at shocks. Since shocks also convert kinetic into internal energy, the correlation between $e_{\text{t}}$ and $e_{\text{int}}$ is readily explained.

It has been verified that the correlation between shocks and turbulence injection is not biased by numerical effects. In Figure 6 we also show the time evolution of $e_{\text{t,filt}}$, computed from $e_{\text{t}}$ by filtering out the cells belonging to the numerically smeared shocks (detected with the finder from Skillman et al. 2008). The evolution of $e_{\text{t,filt}}$, both in the core and outskirts, is basically identical to that of $e_{\text{t}}$, and the same is true for the other energy components (not shown here). This result does not come unexpected, because the performance of the turbulence SGS model at shocks in compressive flows has been extensively tested by Schmidt & Federrath (2011).

3.5 Turbulent velocities on different scales

The numerical values from Figure 6 are better readable, and more directly comparable with observational data and predictions for cluster turbulence, when they are expressed in terms of velocity, rather than energy. This is done in Table LABEL:turbvel, where the velocities computed from the mass-weighted values of the energies are reported. We select the values at the redshift with largest $e_{\text{t}}$ and at $z=0$, as representative of a time with actively stirred and one with decaying turbulence, respectively. The averages of $q$ have values which are at least an order of magnitude lower than those of $v$. This is expected, because the velocity fluctuations probed by $q$ are on much smaller scales than those probed by $v$. Indeed, the ratio $R_{\text{l}}$ of length scales between the turbulence injection length scale and the effective spatial resolution of the simulation is approximately of the order of $100$. Assuming a velocity scaling $v\sim l^{n}$, a velocity ratio

[TABLE]

is expected, with $n=1/3$ or $1/2$ for the incompressible, Kolmogorov or compressible, Burgers turbulence, respectively (Kolmogorov 1941; Frisch 1995; Federrath 2013; cf. also Section 3.7). Given the flow properties in the cluster core and the outskirts, the two cases can be seen as extremes in our problem, with the former more relevant to the smaller length scales. For these two cases, Equation (10) provides $v/q\simeq 5$ and $10$, respectively, to be compared with a value around $15$ (see Table LABEL:turbvel). The agreement with the simulation data within a factor of $2$ to $3$ is satisfactory, considering that this just as a rough estimate, and that large-scale velocities have a component of laminar motions in addition to the turbulent one.

3.6 Evolution of vorticity and volume filling factor

In Figure 7, the evolution of enstrophy $\mathcal{E}$ in the cluster core and outskirts is shown. This quantity (see also Section 2.1) is defined from the vorticity magnitude as

[TABLE]

At low redshift the enstrophy has peaks at the same times as the internal and SGS turbulence energy (Figure 6) both in the core and the outskirts. Analogously to those quantities, it appears therefore sensitive to the stirring on the smallest resolved scales. The correlation between $\mathcal{E}$ and $e_{\text{t}}$ is expected, as demonstrated by Iapichino et al. (2011). Different than in the time evolution of $e_{\text{t}}$, the two low-redshift maxima of $\mathcal{E}$ appear of similar magnitude. Moreover, other boosts of $\mathcal{E}$ are visible at higher redshift and are related to the early cluster buildup. The magnitude of the enstrophy peaks appears to correlate only weakly with the masses of the involved merging structures.

Another use of the vorticity is in the characterisation of the turbulent flow. Complementary to its intensity, an important property of turbulence in and around galaxy clusters is its volume filling factor $f_{\omega}$. It is difficult to find an objective criterion for the definition of this quantity and here, mainly for sake of consistency with previous works, the definition of Kang et al. (2007) is adopted, namely a cell $i$ is flagged as "turbulent" if the vorticity in the cell

[TABLE]

where $t_{\text{age}}(z)$ is the age of the universe at redshift $z$ and $N$ is a free parameter. According to this equation, at redshift $z$ we define as turbulent the gas that, within the cluster lifetime (represented by $t_{\text{age}}(z)$), has a sufficiently large number $N$ of eddy turnovers. At $z=0$ this definition is equivalent to $\omega_{i}>NH_{0}$ (Miniati, 2014). For a better comparison with the cited analyses, we set $N=10$. The volume filling factor $f_{\omega}$ in a given region is the volume fraction where the condition expressed by Equation 12 is fulfilled.

In Figure 8, the time evolution of this quantity is shown for the cluster core and the outskirts (defined as in Section 3.4). The filling factor is substantial in the cluster core, where according to the definition above it is always above $90\%$. A small increase can be seen during the merger events from $z=0.515$. The volume filling factor is quite remarkable also in the outskirts, where it is mostly larger than $60\%$. In the outskirts the evolution at low redshift is decreasing, closely resembling the evolution of $e_{\text{t}}$ (Figure 6). This is likely to be related to the merger history of the cluster, with the last substantial mergers around $z=0.275$ and the subsequent decay of small-scale turbulence.

In Table LABEL:vort-diff one can see the differences of the gas velocities on different scales, when computed in cells within the "turbulent” volume or outside of it, according to the definition in Equation 12. Since most of the core volume is turbulent, we limit the comparison to the more significant case of the outskirts. The values of $v$ and $q$ in the turbulent gas are very similar to the ones in Table LABEL:turbvel. In the non-turbulent gas the values of $v$ are between $10\%$ and $20\%$ smaller than in the turbulent gas, whereas for $q$ the difference is more pronounced. It is therefore evident that $q$ is a more sensitive diagnostic of the turbulent state of the flow than $v$. This can be understood by considering that $q$ correlates well with the vorticity of the flow (Iapichino et al., 2011), and that $v$ includes a laminar velocity component (bulk flow) besides the turbulent one.

3.7 Diagnostics of resolved and unresolved turbulence: radial profiles

In Section 3.4 several variable definitions, useful for characterising the properties of turbulence in the ICM and in the cluster outer regions have been introduced. In this way, it has been verified that $e_{\text{t}}$ is a good diagnostic of turbulence on small length scales. In this Section we want to compare it with another typical indicator, often employed in simulations to characterise the gas flow, namely the mass-weighted root mean square (henceforth rms) gas velocity $v_{\text{rms}}$. This variable is computed in radial profiles around the cluster centre and defined as

[TABLE]

where $m_{i}$ is the mass contained in the cell $i$, with the summation on $i$ computed over the cells belonging to the spherical shell centred on $r$, and that on $j$ over the spatial directions; $\langle\mn@boldsymbol{v}(r)\rangle$ is the average velocity in the shell centred at radius $r$. A time sequence of profiles of $v_{\text{rms}}$ for the run $OD4$+ is presented in Figure 9$a$. For better readability, we included only two profiles in pre-merger and post-merger phase ($z=0.61$ and $z=0$, respectively) and a sequence of profiles following the major merger. Especially in the latter ones one can see boosts of $v_{\text{rms}}$, within the virial radius, which then slowly decay. The magnitude of $v_{\text{rms}}$ is of the same order of $v$, indicative of the flow on large length scales.

The profiles of $v_{\text{rms}}$ (and particularly the ones corresponding to relatively quiet late phases of the evolution) grow with increasing distance from the cluster centre by about a factor of two. This behaviour has been often misinterpreted as a sign of more intense turbulent flow in the cluster outskirts. Indeed, as also observed by Valdarnini (2011), the meaning of the radial profiles of $v_{\text{rms}}$ is ambiguous, because no filtering scale for turbulence is explicitly included in Equation (13), and this is basically equivalent to the size of the spherical shell used for the computation. The outermost shells are therefore dominated by the large-scale (order of Mpc) bulk flows and do not convey information on the turbulent gas motions on smaller scales.

It is interesting to distinguish which fraction of $v_{\text{rms}}$ is radial, namely caused by infall of gas driven by gravitational accretion (e.g. Klessen & Hennebelle 2010; Federrath et al. 2011b), and which remaining component is most genuinely contributing to turbulence motions. For example, Nelson et al. (2014) show that the radial component of the rms gas velocity grows in the outskirts (see also Iapichino & Niemeyer 2008, for a radial velocity profile in a relaxed cluster). This could explain the growth of $v_{\text{rms}}$ in the cluster outer region, as due mainly to accretion. This indication should be read with a cautionary remark: during a major merger, like the crucial phases described in this work, a strict spherical symmetry is hard to hold, and the idea of isolating a radial component in $v_{\text{rms}}$ is questionable.

By its very nature, the turbulent SGS energy is a scale-dependent quantity. In some problems, it would be quite helpful to define a quantity where the dependency on the spatial resolution (more technically, on the AMR level of the grid) is properly accounted and tentatively corrected. To achieve this, Maier et al. (2009) define a "scaled" SGS turbulence energy, obtained by assuming a local Kolmogorov scaling for the turbulence in the ICM and modifying the value of the SGS turbulence energy in a grid cell depending on its AMR level. According to the assumption of Kolmogorov scaling in an incompressible turbulent flow,

[TABLE]

where $v_{1}$ and $v_{2}$ are turbulence velocities on different length scales $l_{\Delta,1}$ and $l_{\Delta,2}$, both within the inertial range of turbulence333In contrast to the Kolmogorov (incompressible) turbulence, where the scaling follows $v\sim l^{1/3}$, highly supersonic, Burgers (compressible) turbulence follows a stronger scaling with $v\sim l^{1/2}$ (Federrath, 2013).. In our numerical scheme, $v_{1}$ and $v_{2}$ can now be interpreted as turbulent velocities on subgrid scales, and $l_{\Delta,1}$ and $l_{\Delta,2}$ as grid resolutions corresponding to the levels of refinement $l_{1}$ and $l_{2}$. Using this interpretation of Equation (14), we implicitly assume that this statistical relation holds locally.

Based on Equation (14) and on the previous assumption, we see that the SGS turbulent energies on two different levels of refinement $l_{1}$ and $l_{2}$ with cell size $l_{\Delta,1}$ and $l_{\Delta,2}$, respectively, are related by:

[TABLE]

Thus, taking $l_{\Delta,\text{max}}$ as the effective spatial resolution at the finest AMR level $l_{\text{max}}$, we define

[TABLE]

as the turbulence SGS energy that can be computed from $e_{\text{t}}$ based on Equation (15), if one uses a Kolmogorov scaling from the generic length scale $l_{\Delta}$ to $l_{\Delta,\text{max}}<l_{\Delta}$. Here one makes use of the relation between the coarse grid resolution $l_{\Delta,\text{coarse}}$ (at $l=0$) and the resolution $l_{\Delta}$ at the generic AMR level $l$: $l_{\Delta}=l_{\Delta,\text{coarse}}/2^{l}$. In this relation it is assumed that, by setup definition, the resolution changes by a factor of 2 within neighbouring levels of refinement.

To provide a comparison term, in Figure 9$b$ we present the profiles of the velocity $q_{\text{scaled}}$, computed from the scaled SGS turbulence energy $e_{\text{t,scaled}}$ as $q_{\text{scaled}}=\sqrt{2e_{\text{t,scaled}}}$. In this definition, only turbulence on small (in this case, sub-grid) scales contributes, and the bulk flow on larger length scales is filtered out. One can clearly see that the profiles of $q_{\text{scaled}}$ are relatively flat, up to large distances from the cluster centre. Moreover, this variable is quite useful to highlight merger-induced turbulence, like the stirring caused by the merger shock, visible as the peak at $r=1.2$ Mpc $h^{-1}$ at $z=0.275$, moving outwards in the two subsequent times. The displacement and the decrease of the peak in $q_{\text{scaled}}$ in different profiles is related to the decrease in time of $e_{\text{t}}$ seen in the outskirts, in Figure 6. One can therefore conclude from this figure that, fixing a given length scale (in our case, of the order of the effective spatial resolution of our simulation sample) the SGS turbulence energy does not systematically increase in the cluster outskirts, as also shown by Valdarnini (2011), Vazza et al. (2012) and Schmidt et al. (2014) for other diagnostics of turbulence on small length scales.

4 Discussion

4.1 Turbulence in galaxy clusters

A merger event injects turbulence in the ICM. One of the main results of this work is to show that this stirring mechanism dominates not only in the cluster core, but also in the outer regions up to several virial radii. Although at first it appears intuitive, this is in contrast with earlier, less resolved simulation results by Iapichino et al. (2011), and to the theoretical analysis of the outskirts at low redshift by Cavaliere et al. (2011), which considers the effect of gas accretion onto the clusters, but not their merger activity. In both those studies, the energy in turbulent motions in the cluster peripheries increases with time to the current epoch. The same increase of turbulence in the outskirts, tracked by vorticity, is also observed in the recent simulation of Vazza et al. (2017). It is hard to a attempt a comparison on this point with that work, which refers to a single cluster undergoing a major merger but with smaller mass (and thus different energy budget involved) than this study. The analysis of a larger simulation sample would be beneficial on this problem.

We have seen in Section 3.4 that turbulence diagnostics on different length scales do not have the same time evolution. In Figure 6 the values of $e_{\text{kin}}$ start increasing earlier than $e_{\text{t}}$ as a consequence of the major merger. We interpreted it as an effect of stirring on large scales, which only in a later phase is imprinted on small (subgrid) scales. Interestingly, in the core $e_{\text{kin}}$ starts increasing even before the major merger (according to the hop tool, the merging structure is identified as a single cluster at $z=0.427$). Moreover, a series of minor merger events follow the major one (Section 3.1). The resulting picture is in agreement with our knowledge of the cosmic large-scale structure (cf. Vazza et al. 2011b; Miniati 2015): a cluster accreting along a filament is accompanied by minor clumps and by the filament gas itself. The dissipation of turbulence in clusters is therefore a more complex process than the decay in forced, idealised simulations: besides the main driver of turbulence (in our case, the major merger), there are further minor stirring agents anticipating and following it.

A corollary of the complex process described above is that it is difficult to provide a unique and consistent definition of the turbulence dissipation time in clusters. The duration of the turbulence decay (i.e. how long the cluster appears turbulent) depends on multiple driving events. Especially from Figure 7 one can visually estimate that the dissipation time of a single driving event is less than $2$ Gyr. However, during a merger the turbulent phase can be much longer. As an example in use of the tools presented in this work, we estimate here the turbulence decay time $t_{\text{decay}}$ in the core after the double merger as the time it takes for $e_{\text{t}}$ to decrease back to the value it had at $z=0.427$ (when the major merger takes place). This occurs at $z=0.15$, resulting in $t_{\text{decay}}=2.7$ Gyr in the cluster core. In the outskirts the estimate is complicated by the fact that the merger does not occur in that region. A visual inspection of Figure 6 suggests that an appropriate choice for $t_{\text{decay}}$ is the time between $z=0.35$ (corresponding to the peak of $e_{\text{kin}}$) and $z=0$, namely $4.0$ Gyr. Just to provide the comparison with a typical timescale related to turbulence, the turbulence eddy turnover time is $t_{\text{eddy}}=L/v$, where $L$ is the integral length scale (in the cluster case, of the order of $500$ kpc $h^{-1}$) and $v$ is the typical velocity on the scale $L$, around $1000$ km s*-1* from Table LABEL:turbvel. The eddy turnover time in the ICM is therefore of the order of $0.7$ Gyr, and can be better compared with the dissipation time of a single driving event, estimated above.

Paul et al. (2011) noticed that the relatively long decay time for turbulence in clusters can be problematic to reconcile with the statistics of observed radio halos, in particular with the sharp bimodality between clusters hosting halos or not (Brunetti et al., 2007, 2009). As a possible solution, the acceleration efficiency in the turbulent re-acceleration model (Brunetti & Lazarian, 2011) should be a steep function of the turbulent kinetic energy. Recently, Cassano et al. (2016) performed an explorative study on statistics of radio halos and fraction of merging galaxy clusters, and put constraints on the timescale of merger-induced disturbance which are similar to the turbulence decay timescale derived in this work.

Although the average values of turbulent velocity are similar in the core and the outskirts at all length scales (as seen in Table LABEL:turbvel and discussed in Section 3.7), the kinetic energy content in the latter regions is larger compared to the internal energy, because these regions have colder gas than the core. This is also clearly shown in terms of the Mach number of the turbulent flow on large length scales, defined here for simplicity as

[TABLE]

where $v$ is computed from $e_{\text{kin}}$ and reported in Table LABEL:turbvel, $c_{\text{s}}$ is the sound speed, $\gamma=5/3$ and the value of $e_{\text{int}}$ is the mass-weighted average computed as in Figure 6. The Mach number of the flow is smaller in the cluster core than in the outskirts. The flow is mildly subsonic in the core ($\mathcal{M}_{\text{turb}}=0.52$ at $z=0$, $0.75$ at $z=0.350$), and more compressible in the outskirts ($\mathcal{M}_{\text{turb}}=1.30$ at $z=0$, $1.32$ at $z=0.275$). Both the velocity values and the Mach numbers are very similar to the recent analysis on thermal Sunyaev-Zel’dovich (SZ) fluctuations in the Coma cluster, performed by Khatri & Gaspari (2016). Recent results (Zinger et al., 2016) have drawn attention to a point often only hinted (e.g. Paul et al. 2011; Iapichino & Brüggen 2012), namely the role of gas and substructures, accreting along cosmic filaments, in driving turbulence (see also Klessen & Hennebelle 2010, for a more general discussion of accretion as ubiquitous driver of turbulence on a wide range of scales). In a detailed study of this process, a refinement strategy like the one used in this work is mandatory, if one wants to correctly capture the evolution of mildly overdense objects and their role in the cluster energy budget. The interactions between the filaments and the outgoing merger shocks can shape the latter in a way that can reproduce complex observed morphologies, like for example the Toothbrush relic (van Weeren et al., 2012). Without any attempt to compare with those data and their interpretation, the straight appearance of the upper merger shock at $z=0.275$ and $0.21$ (Figure 2) is suggestive.

4.2 Flow in the outskirts and connection to radio relics

As mentioned in the Introduction, there is recent evidence that the acceleration model in relics is in tension with the simple DSA, and that one should invoke a pre-existing CR population, either produced by earlier stirring events or ejected by active galaxies (Shimwell et al., 2015; Stroe et al., 2016). Other works like Fujita et al. (2015) study a somewhat complementary scenario, namely turbulence driven in the region downstream the propagating merger shock (cf. Paul et al. 2011), and subsequent re-acceleration of CRs. We notice however that simple theoretical arguments (Iapichino & Brüggen, 2012) show that the flow in the post-shock region cannot sustain fully developed turbulence.

In principle, the propagation of shocks can expose the upstream conditions of the ICM, and its inhomogeneities. The volume filling factor of the turbulent flow and its turbulent energy could then be linked, through appropriate modelling, with the properties and distribution of a possible seed CR population, that can be re-accelerated during merger events. Turbulence might play a similar role also for the amplification of an upstream magnetic field, further noticing that the driving mechanism (consisting of a mixture of solenoidal or compressive modes) can greatly affect the dynamo action (Federrath et al., 2011a; Schober et al., 2013). Concerning the spatial distribution of the turbulent flow, in Figure 8 it was seen that its volume filling factor $f_{\omega}$ (here defined by a threshold in vorticity equivalent to about $10$ turnovers during the cluster lifetime, see Equation 12) is larger than $60\%$ in the outskirts at $z=0$, and approaches unity in the core. These findings compare well with the high volume-filling fraction of turbulent flows resolved in the simulations by Miniati (2014) and Vazza et al. (2010) and, on the other hand, show that a merger-driven shock has a non-negligible probability of interaction with a medium having a low turbulence energy content.

This whole argument suggests that the fraction of turbulent flow in the outskirts and the level of turbulence would be interesting additions in the models of radio relics for future simulations. These points have never been explicitly considered in previous studies (Hoeft et al., 2008; Bonafede et al., 2012; Skillman et al., 2013).

4.3 Caveats and limitations

After having shown the usefulness of our simulation approach, one should however observe that not in all problems a high resolution level in the cluster outer regions is strictly required. As shown by Maier et al. (2009), many general properties like radial profiles of thermodynamical variables do not critically depend on the resolution strategies tested in this work. The same is true for variables averaged on relatively large volumes, like the kinetic energy definition introduced in Section 3.4. On the other hand, if one is interested in turbulence injected in environments with low overdensity like post-shock regions in the outskirts or cosmic filaments, resolving small-scale flow requires an AMR approach similar to ours, or resorting to other strategies like sufficiently fine static grids (Miniati, 2014).

Moreover, in this work we did not take into account additional physics like radiative cooling and AGN feedback. While the latter is certainly a stirring agent on length scales of the order of $10$ kpc (Vazza et al., 2012; Mukherjee et al., 2016), cooling tends to form denser substructures, which are more effective in producing turbulent wakes in their motion (e.g. Iapichino et al. 2013). Both effects are worth being considered in future work.

Another limitation of this work is its focus on one single galaxy cluster realisation. While a large cluster sample would be certainly beneficial, we tried to critically guide our discussion to general properties of major mergers, rather than to the peculiarities of the simulated case. Some variance in the presented results can be retrieved by the comparison with other single clusters presented in the literature (Iapichino & Niemeyer, 2008; Maier et al., 2009; Miniati, 2014, 2015; Vazza et al., 2017).

A recent work by Zhuravleva et al. (2014) revives an idea which was proposed earlier (e.g. Dennis & Chandran 2005), namely the role of turbulent heating to offset radiative cooling in galaxy clusters hosting a cooling flow. Our simulations neither model radiative cooling nor AGN outflows, but it should be stressed again that the SGS turbulence energy is a suitable variable to evaluate the role of turbulence on small length scales, and that future studies on this problem could profit from it.

5 Summary and conclusions

In this work we present the analysis of a suite of cosmological grid-based simulations, following the evolution of a galaxy cluster undergoing a major merger. In the simulations, a subgrid scale model for the computation of unresolved turbulent energy is used. The runs differ from each other in their AMR strategy; in particular, the properties of a run using refinement based on regional variability of vorticity are compared with the widely used standard method of refinement based on overdensity. Here we focus on studying the stirring of turbulence in the cluster outskirts and provide a refinement criterion that adequately captures the turbulence in the outskirts. Our main conclusions are summarised as follows:

The refinement based on vorticity is suitable for reaching a good level of resolution outside the cluster core, when compared with strategies based only on gas and DM overdensity. The computational volume refined at the maximum AMR level in run $OD4$+ is a factor of $5.3$ larger than in run $OD2$, performed using the most permissive refinement threshold based on DM and gas overdensity (see Table LABEL:tab1sims for the naming convention of the simulations). Similarly, the volume-weighted average spatial resolution at $r=2\,R_{\text{vir}}$ in run $OD4$+ is a factor of $1.74$ better than in run $OD2$. The simulation employing refinement on vorticity is able to better resolve underdense structures at the cluster periphery, like cosmic filaments, and the small-scale stirring associated with the gas inflow onto the ICM. This performance of mesh refinement cannot be obtained with any feasible combination of thresholds of the overdensity criteria only (cf. Figure 4). 2. 2.

The turbulence SGS model provides a quantity, the SGS turbulence specific energy $e_{\text{t}}$, which is a useful indicator of small-scale turbulence and its time evolution during a merger event, in combination with the kinetic energy used as a large-scale diagnostic. This variable is a viable alternative to other approaches of characterising turbulence in galaxy clusters, like the ones in Sections 3.4, 3.6 and 3.7, and to other filtering methods (Lau et al., 2009; Maier et al., 2009; Vazza et al., 2011b; Vazza et al., 2012; Vazza et al., 2017). Cluster radial profiles of this quantity show a good correlation with boosts associated with the propagation of a merger shock, and at large distances from the cluster centre are not affected by large-scale bulk flows like the rms gas velocity $v_{\text{rms}}$. 3. 3.

The evolution of the energy budget of the cluster studied is dominated by the major merger event, not only in the cluster core, but well beyond the virial radius (our analysis has been extended to $2\,R_{\text{vir}}$, or $4.3$ Mpc $h^{-1}$ at $z=0$), where the merger shock propagates, heating the gas and injecting turbulence in the flow. The timescale for the decay of turbulence is of the order of several $10^{9}$ years, because of the complex features of the structure buildup, involving multiple submergers at different scales. The evolution of the turbulence energy in the outskirts does not show any increase at low redshift caused by accretion of pristine gas, as previously suggested (Iapichino et al., 2011; Cavaliere et al., 2011). A firmer conclusion on this point will require the study of more simulations with clusters of different dynamical evolution. 4. 4.

On length scales of $10$ kpc $h^{-1}$, of the order of the resolution scale of the simulation, and of the scale best probed by our turbulence SGS model, the turbulence velocity (represented by $q_{\text{scaled}}$, Section 3.7) is similar both in the cluster core and in the outskirts (Figure 9$b$), but the flow is more compressible in the outskirts, because of the radially decreasing thermal energy profile. 5. 5.

The volume filling factor of flow with vorticity larger than $10/t_{\text{age}}(z)$ is around $60\%$ at low redshift in the cluster outskirts. It is speculated that the volume fraction and the energy content of turbulent flow can be relevant for the theory of radio relics: useful estimates of magnetic field and pre-accelerated CR populations can be derived, through adequate modelling, from the preshock conditions.

acknowledgements

L.I. thanks W. Schmidt and J. Niemeyer for their primary contribution in the development of the turbulence SGS model. The numerical simulations were carried out on SuperMUC of the Leibniz Supercomputing Centre in Garching (Germany), in the framework of the project pr95he. C.F. acknowledges funding provided by the Australian Research Council’s Discovery Projects (grants DP150104329 and DP170100603). C.F. further acknowledges the Jülich Supercomputing Centre (grant hhd20), the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. R.S.K. thanks for support from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013) via the ERC Advanced Grant "STARLIGHT: Formation of the First Stars" under the project number 339177. R.S.K. furthermore acknowledges funding from the Deutsche Forschungsgemeinschaft via SFB 881, "The Milky Way System" (sub-projects B1, B2 and B8) and from the SPP 1573 "Physics of the Interstellar Medium". The enzo code (http://enzo-project.org) is the product of a collaborative effort of scientists at many universities and US national laboratories. Most of the data analysis and visualisation was performed using the yt toolkit (Turk et al., 2011).