THE EFFECT OF ANISOTROPIC VISCOSITY ON COLD FRONTS IN GALAXY CLUSTERS

Our simulations solve the following set of MHD equations:

$$\begin{eqnarray} \frac{\partial {\rho }}{\partial {t}} + \nabla \cdot (\rho {\bf v}) &=& 0, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \frac{\partial {(\rho {\bf v})}}{\partial {t}} + \nabla \cdot \left(\rho {\bf vv} - \frac{\bf BB}{4\pi } + p{\mathsf I}\right) &=& \rho {\bf g} - \nabla \cdot { \Pi }, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} && \frac{\partial {E}}{\partial {t}} + \nabla \cdot \left[{\bf v}(E+p) - \frac{{\bf B}({\bf v \cdot B})}{4\pi }\right] = \rho {\bf g \cdot v} - \nabla \cdot {\bf Q} \nonumber \\ &&\quad\quad - \nabla \cdot ({ \Pi } \cdot {\bf v}), \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \frac{\partial {\bf B}}{\partial {t}} + \nabla \cdot ({\bf vB} - {\bf Bv}) &=& 0, \end{eqnarray} \tag{ 4 }$$

where

$$\begin{equation} p = p_{\rm th} + \frac{B^2}{8\pi } \end{equation} \tag{ 5 }$$

is the total pressure,

$$\begin{equation} E = \frac{\rho {v^2}}{2} + \epsilon + \frac{B^2}{8\pi } \end{equation} \tag{ 6 }$$

is the total energy per unit volume, ρ is the gas density, v is the fluid velocity, p_th is the gas pressure, is the gas internal energy per unit volume, B is the magnetic field vector, and $\mathsf {I}$ is the unit dyad. We assume an ideal equation of state with γ = 5/3, equal electron and ion temperatures, and primordial abundances with molecular weight $\bar{A} = 0.6$.

We model the viscous transport of momentum as either isotropic or anisotropic. The viscous flux in the isotropic case is given by

$$\begin{equation} {\Pi }_{\rm iso} = -\mu \nabla {\bf v}, \end{equation} \tag{ 7 }$$

where μ is the dynamic viscosity coefficient. To illustrate the physics behind anisotropic (Braginskii) viscosity, we will follow a more detailed derivation. Equation (2) may be written in an alternative form:

$$\begin{equation} \frac{\partial {(\rho {\bf v})}}{\partial {t}} + \nabla \cdot \left(\rho {\bf vv} + {\mathsf P}\right) = \rho {\bf g}, \end{equation} \tag{ 8 }$$

where ${\mathsf P}$ is the total (thermal+magnetic) pressure tensor:

$$\begin{equation} {\mathsf P} = \left(p_\perp +\frac{B^2}{8\pi }\right){\mathsf I}-\left(p_\perp -p_\parallel +\frac{B^2}{4\pi }\right)\hat{{\bf b}}\hat{{\bf b}}, \end{equation} \tag{ 9 }$$

p_⊥ (p_∥) is the thermal pressure perpendicular (parallel) to the magnetic field, and $\hat{{\bf b}} = {\bf B}/B$ is the unit vector in the direction of the magnetic field. The total thermal pressure satisfies

$$\begin{equation} p = \frac{2}{3} \,p_\perp +\frac{1}{3} \,p_\parallel. \end{equation} \tag{ 10 }$$

Differences in these two components of the thermal pressure arise from the conservation of the first and second adiabatic invariants for each particle on timescales much greater than the inverse of the ion gyrofrequency, $\Omega _g^{-1}$ (Chew et al. 1956). When the ion–ion collision frequency ν_ii is much larger than the rates of change of all fields, an equation for the pressure anisotropy can be obtained by balancing its production by adiabatic invariance with its relaxation via collisions (cf. Schekochihin et al. 2005):

$$\begin{equation} p_\perp -p_\parallel = 0.960\,\frac{p_{\rm i}}{\nu _{\rm ii}}\frac{d}{dt}\ln {\frac{B^3}{\rho ^2}}, \end{equation} \tag{ 11 }$$

where p_i is the thermal pressure of the ions. The ion dynamic viscosity coefficient for the ICM is given by Spitzer (1962), Braginskii (1965), and Sarazin (1988)

$$\begin{eqnarray} \mu &=& 0.960\,\frac{n_{\rm i}k_{\rm B}T}{\nu _{\rm ii}} \\ \nonumber &\approx & 2.2 \times 10^{-15}\frac{T^{5/2}}{\ln \Lambda _{\rm i}}\, {\rm g\, cm^{-1}\,s^{-1}}, \end{eqnarray} \tag{ 12 }$$

where the temperature T is in Kelvin and ln Λ_i is the ion Coulomb logarithm, which is a weak function of ρ and T; for simplicity we follow Roediger et al. (2013a) by approximating ln Λ_i = 40, appropriate for conditions in the ICM. The kinematic coefficient of viscosity ν = μ/ρ. Using Equations (1) and (4) to replace the time derivatives of density and magnetic field strength with velocity gradients, Equation (11) may be written as

$$\begin{equation} p_\perp -p_\parallel = 3\mu \left(\hat{{\bf b}}\hat{{\bf b}}-\frac{1}{3}{\mathsf I}\right):\nabla {\bf v}. \end{equation} \tag{ 13 }$$

It can then be shown that Equations (2) and (8) are equivalent, with the resulting viscous flux for the anisotropic (Braginskii) case being

$$\begin{equation} { \Pi }_{\rm aniso} = -3\mu \left(\hat{{\bf b}}\hat{{\bf b}}-\frac{1}{3}{\mathsf I}\right)\left(\hat{{\bf b}}\hat{{\bf b}}-\frac{1}{3}{\mathsf I}\right):\nabla {\bf v}. \end{equation} \tag{ 14 }$$

These two different fluxes (isotropic, Equation (7), and anisotropic, Equation (14)) are implemented in our different simulations. For the isotropic cases, we model viscosities less than the Spitzer value by including a multiplicative suppression factor f.

We also include anisotropic thermal conduction in some of our runs. The heat flux due to thermal conduction by electrons is given by

$$\begin{equation} {\bf Q} = -\kappa \hat{{\bf b}}\hat{{\bf b}} \cdot \nabla {T}, \end{equation} \tag{ 15 }$$

where the conductivity coefficient (Spitzer 1962; Braginskii 1965; Sarazin 1988)

$$\begin{eqnarray} \kappa &=& 3.2 k_{\rm B} \frac{n_{\rm e} k_{\rm B} T}{m_{\rm e}\nu _{\rm ee}} \\ \nonumber &\approx & 1.84 \times 10^{-5}\frac{T^{5/2}}{\ln \Lambda _{\rm e}}\, {\rm erg\, cm^{-1}\, s^{-1}\, K^{-1}}, \end{eqnarray} \tag{ 16 }$$

where ν_ee is the electron–electron collision frequency and ln Λ_e = 40 is the electron Coulomb logarithm. In this model, the conduction proceeds at full Spitzer rate parallel to the field lines, and is zero perpendicular to the lines. The corresponding thermal diffusivity χ = κT/p. Following Cowie & McKee (1977), we include the effect of conduction saturation whenever the characteristic length scale associated with the temperature gradient is smaller than the electron mean free path (although in the bulk of the ICM this effect is insignificant).

We performed our simulations using Athena 4.1, a parallel conservative magnetohydrodynamic (MHD) astrophysical simulation code (Stone et al. 2008). The MHD algorithms in Athena are detailed in Gardiner & Stone (2005, 2008). The directionally unsplit corner transport upwind integration method and the HLLD Riemann solver are used in all of our simulations, with third-order (piecewise parabolic; Colella & Woodward 1984) reconstruction.

We included anisotropic conduction and Braginskii viscosity following the approaches of Sharma & Hammett (2007), Dong & Stone (2009), Parrish et al. (2012), and Kunz et al. (2012). The conductive and viscous fluxes are implemented via operator splitting, with monotonized central (MC) limiters applied to the fluxes to preserve monotonicity. The latter ensures that unphysical transport of energy and momentum does not occur in the presence of steep gradients.⁵ Since these diffusive processes are modeled by explicit time-stepping methods, they have very restrictive Courant-limited time steps ∝(Δx)². For this reason, we choose to accelerate the calculation via the method of super-time-stepping (STS; Alexiades et al. 1996). STS achieves this speed up by requiring stability of the calculation at a time step that is significantly larger than the Courant-restricted time step (in our case, the time step associated with the diffusive processes of viscosity and conduction). This time step is typically on the order of the hydrodynamical time step associated with the sound speed of the gas. Finally, since we are only concerned in this work with the effects of ICM microphysics on the cold fronts, we do not explicitly include the effects of radiative cooling during the simulation, though our initial condition implicitly assumes that the cool core was originally formed by such cooling.

We model two cluster merger scenarios, both drawn from previous investigations. The first is the setup from Roediger et al. (2011, 2013a), which we dub "Virgo." Virgo is modeled as a cool-core cluster with M ∼ 2 × 10¹⁴ M_☉ and T_X ∼ 2 keV. The gas-less subcluster is modeled as a Hernquist DM halo with a mass of 2 × 10¹³ M_☉. The trajectory results in a closest passage of the subcluster of ≈100 kpc at time t = 1.0 Gyr in the simulation.⁶ The moment identified by Roediger et al. (2011) as most consistent with the shape and size of the Virgo cold fronts is t = 2.7 Gyr in our simulations.

The second merger scenario is identical to our previous setup in ZML11 and Z13, which was originally derived from Ascasibar & Markevitch (2006). In this case, the cool-core cluster is more massive (M ∼ 10¹⁵ M_☉) and hotter (T ∼ 7 keV), resembling A2029 (though not exactly reproducing it). The gas-less subcluster is 5 times less massive than the main cluster and approaches it with an initial impact parameter of b = 500 kpc. We will refer to this setup as "AM06."

The gravitational potential on the grid is the sum of two collisionless "rigid bodies" corresponding to the contributions to the potential from both clusters. This approach to modeling the potential is used for simplicity and speed over solving the Poisson equation for the matter distribution, and is an adequate approximation for our purposes. It is the same approach that we used in ZML11, and is justified and explained in Roediger & Zuhone (2012).

The tangled magnetic field of the cluster is set up in an similar way to the simulations of ZML11 and Z13. We use a radial profile for the magnetic field strength of B = B₀[ρ(r)/ρ₀]^η with η = 0.5 (Bonafede et al. 2010). The average magnetic field strength in the core region corresponds to β = p_th/_B ∼ 1000–1500, which is on the weaker end of the parameter space from ZML11. This low field strength was chosen to ensure the development of K-H instabilities in order to determine the effect of viscosity on their development. For the three-dimensional tangled structure of the initial magnetic field, we follow the approach of Ruszkowski et al. (2007, 2008) and Ruszkowski & Oh (2010), and we refer the reader to those papers for the details. Figure 1 shows the initial radial profiles of density, temperature, and magnetic-field strength for both clusters.

Zoom In Zoom Out Reset image size

For all of the simulations, we set up the main cluster within a cubical computational domain of width L = 4 Mpc on a side. We employ static mesh refinement (SMR), with 256 cells on a side on the top-level domain and five smaller domains at increasing levels of refinement, each half the size of the domain just above it, centered on the cluster potential minimum. This results in a finest cell size of Δx = 0.98 kpc within a region of l³ = (250 kpc)³. Appendix A details the results of a test where we experiment with higher and lower spatial resolutions. Roediger et al. (2013b) examined the effect of varying resolution on the development of K-H instabilities in a simulation of the ICM, and determined that the gross morphologies of the K-H rolls is captured well down to resolutions of 32 cells per perturbation wavelength, a wavelength of λ ∼ 32 kpc in our simulations (for further discussion of the effects of resolution see Section 4.2).

Similar to the approach taken in Z13 for simulations with diffusive effects (viscosity and/or conduction), we switch on these effects at the pericentric passage (t = 1.0 Gyr in the case of the "Virgo" simulations and t ∼ 1.5 Gyr in the case of the "AM06" simulations). This is just before the onset of sloshing, and allows us to examine the effects of these processes on the sloshing cold fronts without any significant change in our initial condition of the main cluster that would occur if these effects (particularly conduction) were switched on.

For these two cluster models, we perform four simulations each with different types of viscosity. Also, for each model we perform two additional simulations with anisotropic thermal conduction, one of the two including Braginskii viscosity. In addition, to see the effect of viscosity separate from that of the magnetic field for the same cluster setup, we also ran unmagnetized runs with and without viscosity of both setups (detailed in Section 4.6). These are all exploratory runs, not intended to cover the full parameter space of the viscosity and magnetic field values; we will try to reproduce particular observed clusters in a future work.

We begin by visually examining the development of K-H instabilities at cold front surfaces in slices through the center of the simulation domain in the merger plane. Figures 2 and 3 show temperature slices for the two clusters at various times in the different simulations. The inviscid simulations quickly develop K-H instabilities at the cold front surfaces, which continue to grow and develop as the simulations progress, though the development of K-H instabilities is suppressed to a degree by the weak magnetic field compared to what would be obtained in an unmagnetized simulation (ZML11 and Section 4.6). These instabilities generate turbulence (to be discussed in Section 3.3) and significant small-scale structure in temperature. In the anisotropic and isotropic f = 0.1 simulations, the same instabilities develop at early times, but as the cold fronts expand, these are progressively damped out by viscosity. We mark the positions of cold fronts with different morphology in the different simulations with arrows. However, there are still some locations along the fronts where there are indications of instability, as well as some small-scale temperature fluctuations within the envelope of the cold fronts. In the case of the isotropic, f = 1 simulations, no K-H instabilities or small-scale temperature structures develop at all during the course of the simulations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 4 and 5 show slices of the magnetic-field strength for the two clusters at the same times for the different simulations. The overall evolution of the magnetic field is very similar in all the simulations, but there are some differences. The magnetized layers that occur along cold front surfaces are disrupted by K-H instabilities in the inviscid simulations, whereas viscosity results in magnetic-field layers that are smoother (see the arrows in the figures), though there is not a perfect wrapping of fields around the cold fronts in any of the simulations (we will examine the consequences of this for anisotropic thermal conduction in Section 3.4). As viscosity is increased, there is less small-scale structure in the magnetic field.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To make closer comparisons with observations of cold fronts in actual clusters, we construct synthetic Chandra observations. We generate X-ray photons from an APEC model (Smith et al. 2001) of a thermal gas with Z = 0.3 Z_☉ for each simulation cell within a radius of R = 250 kpc from the cluster potential minimum.⁷ We also apply Galactic absorption assuming N_H = 10²¹ cm². Our procedure for generating these observations is outlined in detail in Appendix B. To generate an approximate number of counts that would be expected from a mosaic of Chandra observations of each of our simulated clusters, we use the on-axis effective area curve of the ACIS-S3 chip. For the "Virgo" simulations, we use the actual redshift of the Virgo cluster, z = 0.0036, and for the "AM06" simulations, we choose a redshift of z = 0.1. The images have been blocked so that the pixels correspond to the size of the finest simulation cell size, Δx = 0.98 kpc.⁸ To determine the effect of photon statistics, we simulate two different exposures, 300 ks and 30 ks.

Figures 6 and 7 show the X-ray counts images of the two model clusters at different epochs and for different exposure times. Projection effects and Poisson noise make the interpretation of these images more difficult than for the temperature slices. It is clear that long exposures (in our case, ∼300 ks) are needed to determine the amount of K-H disruption at the cold front surfaces. In the 30 ks exposures, it is much more difficult to make distinctions between the cold fronts in the simulations due to poor statistics. For the rest of this section, we will focus on the 300 ks exposures.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the inviscid simulations, the presence of K-H instabilities along the cold front surfaces is readily apparent from the X-ray images. By contrast, in the isotropic, f = 1 simulations, the cold front arcs are completely smooth in appearance. In the isotropic, f = 0.1, and anisotropic simulations, while the cold front surfaces are mostly smooth, small indications of instability are visible. Some fronts that appear in the viscous simulations are so disrupted in the inviscid simulations that they are no longer technically cold fronts (one example is in the "AM06" simulations, indicated by the cyan arrow in Figure 7). Another indicator of cold front instability is the overall shape of the cold fronts. Unstable fronts tend to have a "boxy" shape, due to the development of the instability over the curved surface, whereas fronts stabilized by viscosity appear more round (a fact previously pointed out by Roediger et al. 2012b). This is most evident in the "AM06" case.

The characteristics of the cold front structures may be seen more clearly by looking at the surface brightness residuals instead of the surface brightness itself. Figure 8 shows the surface brightness residuals for both models at a late epoch, generated by computing an azimuthally averaged radial profile of the counts from the 300 ks observation and subtracting the corresponding image from the original image. As in Figures 6 and 7, the arrows mark the positions of cold fronts with different morphologies in the different simulations. The difference in the degree of K-H disruption can be more clearly seen. In particular, it is obvious from the right panel of Figure 8 that the cold front at the inner radius (indicated by the cyan arrow) is essentially impossible to discern even in the residual image in the inviscid case. Figures 9 and 10 show two "close-ups" of several quantities side by side in the region of cold fronts where instabilities have been suppressed due to anisotropic viscosity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Viscosity can have an effect on turbulence driven by sloshing by damping out turbulent motions on the smallest length scales. Future X-ray telescopes such as Astro-H⁹ and Athena+¹⁰ will be able to measure turbulence and bulk motions from the line shift and line widths of spectral lines (Inogamov & Sunyaev 2003; Sunyaev et al. 2003). It is instructive to determine the effect of our different viscosity models on the turbulence driven by sloshing.

A straightforward way to distinguish the turbulent velocity field between the different models is to compare the velocity power spectra. The sloshing motions themselves complicate matters, since the velocity jumps across the cold fronts impose a ∼k⁻² signature on the power spectrum at all wavenumbers k. Since the bulk of the motion in the simulation is in the sloshing motions, this signal dominates. This reduces to the more general problem of separating turbulence from bulk motions in hydrodynamic simulations, a subject of previous investigations (Vazza et al. 2012; ZuHone et al. 2013a). To carry out this decomposition, we use the procedure outlined in Vazza et al. (2012), which uses a multi-scale filtering scheme to iteratively converge on the mean velocity field around each cell, which is then subtracted from the total velocity to get the turbulent component. As noted by ZuHone et al. (2013a), carrying out the filtering procedure naively on the velocity field from sloshing cold fronts will not completely remove the effect of the bulk motions, because the velocity differences across the cold fronts will result in the local mean velocity changing discontinuously across these surfaces. Vazza et al. (2012) noted the same difficulty with shock fronts, and proposed filtering on the skew of the local velocity field as a way to avoid these difficulties. We adopt the same approach here for our simulations.

Figure 11 shows the filtered ("turbulent") velocity power spectrum as a function of wavenumber k for all of the simulations, as well as the unfiltered velocity power spectrum for the inviscid simulations (for comparison). The anisotropic simulations with the full Braginskii viscosity result in a very modest reduction of the turbulent power, whereas for the isotropic simulations with full Spitzer viscosity the turbulent power is reduced by more than an order of magnitude.

Zoom In Zoom Out Reset image size

Figures 12 and 13 show the emission-weighted projected line-of-sight mean velocity and line-of-sight velocity dispersion along the z-axis of the simulation domain for all of the simulations. For this exercise, we have not filtered out the bulk motions. The projections of the two inviscid simulations show evidence of turbulence in small-scale random velocity fluctuations. In the anisotropic simulations, turbulence is somewhat reduced, but the small-scale velocity fluctuations are still present. This is also the case for the isotropic, f = 0.1 simulations. In the isotropic, f = 1 simulations, essentially all of the small-scale motions are absent.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also find in all of our simulations that significant velocity dispersions may be measured along the line of sight that do not correspond to turbulent motions but instead to large variations in the bulk motion along the line of sight. Though we have determined the dispersion here by simply computing it from the velocity field itself, this effect will also show up in the broadening of X-ray spectral lines, which is how turbulence is to be measured. We stress that care must be taken in interpreting the results of future measurements of the width of spectral lines resulting from turbulence, since it is possible for bulk motions along the line of sight to produce similar signatures. We will examine this effect in the context of a sloshing cool core in more detail in a separate work.

Z13 found that anisotropic thermal conduction acted to eliminate the sharp cold front temperature and density jumps despite the presence of strongly magnetized layers stretched along the fronts. This may be the result of imperfections in the magnetic layers due to field-line tangling from K-H instabilities, resulting in some "leaking" of heat from the hot plasma above the cold front, or to the presence of other regions of hot plasma that are not magnetically isolated from the colder gas inside the fronts (e.g., the spiral-shaped "tongues" of hot gas in between the cold fronts; see Figures 2 and 3). If viscosity acts to prevent the development of K-H instabilities, it is possible that the magnetized layers that develop will be more aligned to the smooth cold front surfaces and hence they may be more effective at suppressing conduction than in the inviscid case. In Figures 2 and 4, we saw two cold fronts where viscosity has suppressed K-H instabilities and produced magnetic fields more aligned to the cold front surfaces. This hypothesis motivates our simulations with anisotropic thermal conduction and anisotropic viscosity.

We performed two additional simulations for each of our model clusters, one with anisotropic thermal conduction (full Spitzer rate along the field lines and zero across the lines; see Section 2.1) and one with both anisotropic thermal conduction and Braginskii viscosity. Figure 14 shows images of the projected temperature constructed at specific epochs for the "Virgo" and "AM06" simulations. In the simulations with anisotropic thermal conduction, though there are still surface brightness gradients at the positions of the cold fronts in the inviscid simulations, these once-sharp gradients have been reduced. This is demonstrated more quantitatively in Figure 15, which shows the extracted surface brightness and projected temperature profiles from the colored annular regions in Figure 14. In the simulations with conduction, the gradients are no longer sharp, with or without viscosity, and the radial profiles are essentially identical. The inclusion of Braginskii viscosity appears to have little effect on the reduction of cold front jumps by conduction.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In our exploratory set of simulations, we have bracketed the range of possibilities for K-H instabilities at sloshing cold fronts for our two examples of low-mass and massive clusters. On the one hand, the cold fronts in our inviscid simulations are very susceptible to K-H instabilities, despite the presence of magnetic fields (this is even more true if magnetic fields are very weak or absent; cf. Section 4.6; ZML11). There appear to be more instabilities in these simulations than we see at cold fronts in actual clusters (a quantitative comparison is forthcoming in our future work). On the other hand, an ICM with an isotropic Spitzer viscosity exhibits no evidence of instabilities at all, which is also contrary to some observations. Our more physically motivated simulations, with either Braginskii viscosity or a reduced isotropic Spitzer viscosity, appear to qualitatively reproduce the observations in terms of the degree of smoothness of the cold fronts. A quantitative comparison with several well-observed cold fronts will be presented in a subsequent publication. However, from our present study, we infer that it may be difficult to distinguish between an isotropic "sub-Spitzer" viscosity and anisotropic Braginskii viscosity using X-ray observations alone. The degree of suppression of K-H instabilities at cold front surfaces and the effect on turbulence are both qualitatively similar. From a purely theoretical perspective, the Braginskii viscosity is the preferred model because it is known that the ICM is magnetized, and the small Larmor radius of the ions will constrain the momentum transport accordingly.

Other avenues for constraining viscosity in the ICM exist. Dong & Stone (2009) used Braginskii–MHD simulations of buoyant bubbles in a cluster atmosphere to investigate their stability properties. They found that bubble stability was impacted most by the magnetic field geometry. In the core of the Perseus cluster, the smooth shape of some of the Hα filaments may indicate the flow around these filaments is laminar and significantly viscous (Fabian et al. 2003). Further simulation studies of these systems and their observations need to be combined with studies of cold fronts to provide better constraints on ICM viscosity.

The different prescriptions for viscosity in our simulations produce significant differences in the amount of turbulence that is produced. To more fully understand the origin of these differences, it is instructive to examine the relevant scales for turbulence in our simulations. The first is the scale of the sloshing motions themselves, L₀, which at early times is roughly L₀ ∼ 50 kpc, increasing at later times to L₀ ∼ 200 kpc. The second scale is the scale at which the turbulent cascade is dissipated by viscosity, l_diss, assuming a roughly continuous generation of turbulence by the sloshing motions as they expand. For dissipation due to Coulomb collisions, this scale depends on the Reynolds number of the plasma (e.g., Roediger et al. 2013a):

$$\begin{eqnarray} {\rm Re} &=& 7 f^{-1} \left(\frac{V_L}{300\, {\rm km \,s^{-1}}}\right)\left(\frac{L_0}{10\, {\rm kpc}}\right) \nonumber\\ && \times \left(\frac{n_{\rm th}}{3 \times 10^{-3} {\rm \,cm^{-3}}}\right)\left(\frac{T}{3\, {\rm keV}}\right), \end{eqnarray} \tag{ 17 }$$

where V_L is the characteristic turbulent eddy velocity on the length scale L₀, n_th is the thermal gas number density, T is the thermal gas temperature, and f is the viscous suppression factor. For Kolmogorov turbulence, l_diss ∼ L₀(Re)^−3/4. For our "Virgo" setup, Re ∼200 and l_diss ∼ 2 kpc, and for the "AM06" setup, Re ∼100 and l_diss ∼ 3 kpc, in the relevant region of the cold fronts (assuming f = 1 and a value of L₀ ∼ 100 kpc for the length scale of the sloshing motions). Since these values are very close to the size of our smallest cells, the numerical dissipation dominates, since its effect sets in at a length scale ∼8Δx (ZuHone et al. 2013a).

It therefore may seem surprising at first that the different simulations have such wide differences in the amount of turbulence that is produced, since all of the simulations have the same cell size, and in any case the turbulent dissipation from viscosity does not occur except at the smallest scales. However, the relevant scale for generating turbulence, as well as for disrupting the cold front surfaces, is not the dissipation scale of the turbulent cascades, but rather it is the wavelength of the smallest perturbations that will not be damped out by viscosity. In an extensive set of plane-parallel simulations of the K-H instability, Roediger et al. (2013b) determined this scale to be

$$\begin{eqnarray} l_{\rm crit} &=& {\rm 30\; kpc}\left(\frac{\rm Re_{crit}}{30}\right)f\left(\frac{V_L}{400\, {\rm km\, s}^{-1}}\right)^{-1}\nonumber \\ &&\times \left(\frac{n_{\rm th}}{10^{-3} {\rm \,cm}^{-3}}\right)^{-1}\left(\frac{T}{2.4\, {\rm keV}}\right)^{5/2}. \end{eqnarray} \tag{ 18 }$$

Taking our lead from Roediger et al. (2013b), we assume the critical Reynolds number at which K-H instabilities develop is Re_crit ∼ 30. This criterion yields roughly l_crit ∼ 40 kpc for the outermost cold fronts at late epochs ("Virgo," t = 2.7 Gyr; "AM06," t = 3.25 Gyr) in both of our models assuming full Spitzer viscosity (f = 1). For f = 0.1, l_crit is an order of magnitude smaller for each model, and a wider range of perturbations that we can resolve in our simulations can become unstable. For Braginskii viscosity, the effective f may be as high as 1/5 for an isotropically tangled field or even smaller (Nulsen & McNamara 2013), depending on the local magnetic field geometry.

If we examine the cold fronts in the early stages of sloshing (upper panels of Figures 2 and 3), we find that in all cases except, in the isotropic, f = 1 simulations, that K-H instabilities develop very quickly. These instabilities generate turbulence that eventually develops within the sloshing region in these simulations. In the isotropic, f = 1 simulation, instabilities never really develop in the first place, the flow is always quite laminar, and as such there is not much turbulence at later stages. Therefore, the most relevant scale that determines the difference between our simulations is l_crit, the scale of the smallest perturbations that become unstable. This scale is very well-resolved by our simulations if f = 1, l_crit ∼ 40 kpc. In the f = 0.1 case, l_crit ∼ 4 kpc, and this scale is resolved only marginally well (the 8Δx criterion is not fulfilled), so numerical resolution may start to have a significant effect on our results (see Section 2.3).

When the pressure anisotropy violates the approximate inequalities

$$\begin{equation} {-}\frac{B^2}{4\pi } \lesssim p_\perp - p_\parallel \lesssim \frac{B^2}{8\pi }, \end{equation} \tag{ 19 }$$

a situation which is expected to occur readily in the ICM (Schekochihin et al. 2005; Lyutikov 2007; Kunz et al. 2011), rapidly growing Larmor-scale instabilities (namely, the firehose instability on the left side of the equation and the mirror instability on the right side) are triggered and act to regulate the pressure anisotropy back to within its stability boundaries. The lack of finite-Larmor-radius effects in the Braginskii–MHD equations means that the fastest-growing firehose modes occur at arbitrarily small scales and that the mirror mode, an inherently kinetic instability, is not described properly. How these microscale instabilities regulate the pressure anisotropy and, in doing so, place constraints on the allowed transport of momentum and heat remains very much an open question, one that has received increased attention in recent years in both the astrophysical (e.g., Schekochihin & Cowley 2006; Sharma et al. 2006) and solar wind (e.g., Gary et al. 2001; Bale et al. 2009) communities. While it is perhaps premature to judge whether the qualitative evolution found by our Braginskii–MHD description of large-scale sloshing and K-H instabilities in the ICM is contingent upon formulating a rigorous description of the kinetic microphysics, we nevertheless feel obliged to speculate upon how such microphysics may impact our results.

A key finding to emerge from (collisionless) kinetic simulations of the driven firehose instability (e.g., Matteini et al. 2006; Hellinger & Trávníček 2008; Kunz et al. 2014) is that the collisionality of the plasma, supplemented by the anomalous scattering of particles off the microscale fluctuations, adjusts to maintain a marginally firehose-stable pressure anisotropy. In a weakly collisional plasma, this effectively reduces the parallel viscosity from $\sim v^2_{\rm th} / \nu _{\rm ii}$ to $\sim v^2_{\rm A} / S$, where S is the shear frequency of the viscous-scale motions (Kunz et al. 2014; Mogavero & Schekochihin 2014). If this behavior holds true in the ICM, the effective Reynolds number may be in fact larger than that given by Equation (18) and K-H instabilities may develop more easily. Moreover, with the pressure anisotropy microphysically pinned at the firehose stability threshold (the left inequality of Equation (19)), the resulting viscous stress would effectively cancel the magnetic tension—another effect that may ease the development of K-H instabilities. The situation with the driven mirror instability is a bit less clear (see Kunz et al. 2014; Riquelme et al. 2014), in that the pressure anisotropy appears to be regulated not just by anomalous particle scattering but also by an increasing population of resonant particles becoming trapped in magnetic mirrors where the pressure is naturally less anisotropic. It is this trapped population that has been suspected of dramatically reducing the effective thermal conductivity of the ICM (Schekochihin et al. 2008), an effect that would call into question the applicability of the simulations in Section 3.4 to real clusters.

Clearly, more work is needed not only to assess the effects of microphysical kinetic processes on the transport of heat and momentum, but to also converge on a set of well-posed MHD-like equations that can be profitably used to study the large-scale dynamics and thermodynamics of the ICM. For now, we take comfort in the fact that the cold fronts produced in our simulations are also the locations where the magnetic-field strength is very much increased, making the inequalities (Equation (19)) difficult to violate.

We find in this work that the inclusion of Braginskii viscosity fails to have any significant effect on the elimination of cold front temperature and density jumps by anisotropic thermal conduction, despite its effect to reduce the tangling of magnetic field lines along cold front surfaces by K-H instabilities. This is perhaps not surprising, given the value of the Prandtl number of the ICM:

$$\begin{equation} {\rm Pr} \equiv \frac{\nu }{\chi } = 0.5 \, \frac{\ln \Lambda _{\rm e}}{\ln \Lambda _{\rm i}}\left(\frac{2m_{\rm e}}{m_{\rm i}}\right)^{1/2} \simeq 0.02. \end{equation} \tag{ 20 }$$

This indicates that viscous forces operate on a much longer timescale than thermal conduction. Even if heat conduction is prevented from occurring directly across the cold front interface by a magnetized layer, the heating of the cold side of the front still occurs, since the front is still surrounded on multiple sides in three dimensions by hot gas that it is not magnetically isolated from (cf. Z13).

We have already noted the many points of comparison between our work and that of Roediger et al. (2013a). In the limit that the magnetic field is very weak, modeling the ICM using hydrodynamics with a suppressed Spitzer viscosity would be a sufficient approximation to an MHD simulation with Braginskii viscosity to reproduce the cold fronts as we see them in observations. If the magnetic field is strong and stretched along the cold front surfaces, it will provide an additional source of suppression of K-H instabilities. Without an independent measurement of the magnetic field strength and direction, it would be difficult to determine which mechanism is largely responsible for determining the shape of the fronts. We investigate this matter further in Section 4.6.

The unmagnetized, reduced Spitzer viscosity approach is also limited when it comes to turbulence, since an isotropic viscosity will damp all modes of a turbulent cascade, whereas Braginskii viscosity will damp the magnetosonic modes only. In the sloshing cluster core, the dominant motions are solenoidal, and the resulting turbulent cascade is mainly Alfvénic. In a major merger, with strong compressible turbulence, the turbulent cascade is likely to be dominated by magnetosonic waves. The details of the effect of the different prescriptions for viscosity on this latter scenario are not straightforward to determine from the conclusions of this study, and requires further work. Though in our study the viscous dissipation scale was very small in the cluster core, in a major merger with low densities and high temperatures this scale could be much larger, and so the differences between the two approaches at this scale would be more prevalent.

Suzuki et al. (2013) simulated the effect of isotropic and anisotropic viscosity on a "bullet" of cold gas propagating through a hot, magnetized medium. For the MHD simulations presented in their work, they assumed a uniform magnetic field direction but varied its direction with respect to the cold bullet's velocity (along the x-axis) between the different simulations. Not surprisingly, for isotropic Spitzer viscosity they find that K-H instabilities are completely suppressed. For anisotropic viscosity, the degree of suppression is highly dependent on the magnetic field direction. For magnetic field lines perpendicular to the direction of motion of the bullet, the cold bullet sweeps these lines up into a draping layer very quickly, with the velocity gradients perpendicular to the field lines. As a result, the K-H instability is only weakly suppressed. For magnetic field lines inclined at a 45° angle to the x-axis, more suppression occurs, due to the fact that the magnetic field is more aligned with the velocity gradient initially.

In this work, our initial conditions are more appropriate for conditions in a relaxed cluster core with an initially tangled magnetic field. It is difficult to compare our results with theirs directly, due to the fact that their setup is more akin to a cold front produced by a major merger (such as in the "Bullet Cluster") than that of sloshing motions. However, our results are consonant with theirs in the sense that Braginskii viscosity results in much less suppression of K-H instabilities than does full isotropic viscosity, except in the case where there is significant alignment of the velocity gradients with the magnetic field.

ZML11 showed that shear-amplified magnetic fields can suppress K-H instabilities in ICM cold fronts, even without viscosity, although this is highly dependent upon the strength of the initial magnetic field. The magnetic field in our simulations obviously plays a similar role, although determining whether its influence dominates over that of viscosity is complicated by the nonlinear nature of the problem. To date, most simulations of cold-front formation and evolution that included viscosity (e.g., ZuHone et al. 2010; Roediger et al. 2013a) had not included magnetic fields as well (Suzuki et al. 2013, is a notable exception). For completeness, we have performed unmagnetized versions of our "Virgo" and "AM06" simulations, both without viscosity and with an isotropic Spitzer viscosity reduced by f = 0.1. We also performed an unmagnetized simulation of the "AM06" model with an isotropic Spitzer viscosity reduced by f = 0.2. Sample temperature slices of these runs at different epochs, compared to their magnetized counterparts, are shown in Figures 16 and 17.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the unmagnetized, inviscid simulations, more small-scale perturbations of the cold-front surface develop than in any of our magnetized simulations, since there is no mechanism (aside from the comparatively small numerical viscosity) to suppress them. More coherent large-scale instabilities are able to grow as well, which are evident in the temperature slices.

In the unmagnetized simulations with f = 0.1 Spitzer isotropic viscosity, K-H instabilities are somewhat inhibited (especially the smallest-scale perturbations), but some large-scale features still persist. In the "Virgo" simulation (Figure 16), the cold fronts appear qualitatively similar to those in its magnetized f = 0.1 counterpart. This is not the case, however, in the "AM06" simulation (Figure 17), as its cold fronts appear more disrupted by K-H instabilities than its magnetized f = 0.1 counterpart. This is especially true at later epochs, on length scales of several tens of kpc. The additional "AM06" unmagnetized run with f = 0.2 shows better qualitative agreement with the magnetized f = 0.1 run. This demonstrates that even the somewhat weak initial field in these simulations (β ∼ 1000–1500; cf. ZML11) can have a non-negligible effect on the appearance of the cold fronts, even in a viscous ICM.

These differences illustrate the difficulty in discerning the precise mechanism for K-H suppression at cold front surfaces from X-ray observations. In general, some combination of the effects of viscosity and amplified magnetic fields are responsible, but these do not combine in a linear way and the separation of these two effects is complicated by the fact that, while the thermal properties of the ICM plasma are well-constrained in the vicinity of a cold front, the local magnetic-field strength is generally unknown (unless it can be constrained by other observations, e.g., radio mini-halo emission, Faraday rotation, or discontinuities in the thermal pressure—for the latter see Reiss & Keshet 2012). Unfortunately, these results indicate that there is not a general answer at this time as to whether the magnetic field or the viscosity is dominant in producing the observed smoothness of many cold fronts.