COLD FRONTS AND GAS SLOSHING IN GALAXY CLUSTERS WITH ANISOTROPIC THERMAL CONDUCTION

X-ray observations of galaxy clusters often show a bright, dense, centrally peaked core comprised of the intracluster medium (ICM) much cooler than the cluster average. The X-ray radiative cooling time, inversely proportional to the square of the gas density, is very short in these dense cores, on the order of several hundred Myr in the very central regions (Sarazin 1986). This gave rise to an early "cooling flow" scenario, in which the gas cools catastrophically in the cluster center and is continuously replaced by the inflow of the surrounding cooling gas. However, Chandra and XMM-Newton have failed to detect the implied large quantities of very cool gas in the cores; in particular, the characteristic emission lines such as Fe xvii emitted by the plasma with T < 1 keV are not observed (Peterson & Fabian 2006). Since X-ray cooling is directly observed, a mechanism is required to heat the ICM in the core to compensate for the cooling. The two most studied heating mechanisms are feedback from active galactic nuclei (AGNs) and thermal conduction. In this paper, we will be discussing the latter mechanism in connection with the phenomenon of gas sloshing.

The role of thermal conduction in the context of heating the cool cores of clusters of galaxies has been studied by many authors (see, e.g., Binney & Cowie 1981; Narayan & Medvedev 2001; Ruszkowski & Begelman 2002; Zakamska & Narayan 2003; Guo et al. 2008). Since in the ICM the Larmor radius of the electrons spiraling along the B ∼ 0.1–1 μG magnetic field lines is many orders of magnitude smaller than the electron mean free path, the conduction is strongly anisotropic, occurring essentially exclusively along the field lines. Therefore, most of these works assumed an isotropically tangled magnetic field geometry and parameterized the heat flux as an effective isotropic thermal conductivity given as a fraction, f_Sp, of the ideal (field-free) Spitzer (1962) heat flux.

However, theoretical and numerical work over the last decade has revealed that the stability and thermal properties of the ICM are possibly more complicated than previously thought if anisotropic thermal conduction is modeled explicitly in the MHD regime. Early analysis by Balbus (2000, 2001) showed that when magnetized atmospheres with anisotropic thermal conduction are considered, their stability depends on the temperature gradient rather than the entropy gradient. He demonstrated that in regions where the temperature increases in the direction of gravity, the plasma is susceptible to the magnetothermal instability (hereafter MTI).

Further analytical work by Quataert (2008) and numerical work by Parrish & Quataert (2008) and Parrish et al. (2008, 2009) has shown that the cluster atmosphere is unstable even when the temperature decreases in the direction of gravity, as in the cluster cool cores (heat-flux buoyancy instability, HBI). The effect of the HBI is to realign the magnetic fields azimuthally $(\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}} \approx 0)$, and thereby significantly reduce the effective radial heat flux Q_r. This renders thermal conduction almost useless against preventing a "cooling catastrophe." Subsequent work has shown that relatively low levels of turbulence can rearrange the field lines in a more isotropic configuration, permitting heat conduction to the core (Parrish et al. 2010; Ruszkowski & Oh 2010, 2011). More recently, and of direct relevance to this work, Lecoanet et al. (2012) simulated the dynamics of Rayleigh–Taylor (RT) stable and unstable contact discontinuities with anisotropic thermal conduction. Of particular relevance is their finding that the temperature gradients of RT-stable contact discontinuities (such as sloshing cold fronts) can be smoothed by anisotropic thermal conduction if the field lines are at least partially misaligned with the surface of the discontinuity. The smoothing of the temperature gradient of the RT-stable contact discontinuities in these simulations was limited by the HBI, which reoriented magnetic field lines initially perpendicular to the front surface in a more parallel direction, suppressing conduction across the surface at later times.

So far, the numerical simulations of the evolution of the cluster atmospheres with anisotropic thermal conduction have mostly been limited to isolated, spherically symmetric clusters initially in hydrostatic equilibrium. These works have been essential in order to determine the character of the instabilities that are present when such physical effects are taken into account. A recent exception to this is the paper by Ruszkowski et al. (2011), which simulated the formation of a galaxy cluster from cosmological initial conditions.

The previous works assumed that the heat conduction along the field lines is at the Spitzer value. However, the effective conductivity along the field lines may be reduced due to strong curvature of the field lines and stochastic magnetic mirrors on small scales (e.g., Chandran et al. 1999; Malyshkin & Kulsrud 2001; Narayan & Medvedev 2001). Markevitch et al. (2003a) attempted to constrain heat conduction in the ICM from the observed temperature differences in the merging, non-cool-core cluster A754. They determined that the effective conduction in the direction of the observed temperature gradients is an order of magnitude less than Spitzer, which, under the assumption of an isotropically tangled field, implies suppression along the field lines. Finally, Schekochihin et al. (2005) suggested that plasma microscale instabilities may be very important in the high-β plasma of the ICM, and consequently its thermal conductivity may be many orders of magnitude smaller than the Spitzer value (Schekochihin et al. 2008).

In this work we consider a scenario where the initial cluster core is relaxed but undergoes a period of gas "sloshing." Many cool-core systems exhibit edges in X-ray surface brightness approximately concentric with respect to brightness peak of the cluster (e.g., Mazzotta et al. 2001; Markevitch et al. 2001, 2003b). X-ray spectra of these regions have revealed that in most cases the brighter (and therefore denser) side of the edge is the colder side, and hence these jumps in gas density have been dubbed "cold fronts" (for a detailed review, see Markevitch & Vikhlinin 2007). These features have been observed in the majority of cool-core systems (Ghizzardi et al. 2010). The presence of one or more cold fronts is seen as an indication of the motion of gas in a cluster. For relaxed systems, the primarily spiral-shaped cold fronts are believed to arise from gas sloshing in the deep dark matter-dominated potential well. These motions are initiated when the low-entropy, cool gas of the core is displaced from the bottom of the dark matter potential well, either by gravitational perturbations from infalling subclusters (Ascasibar & Markevitch 2006, hereafter AM06) or by an interaction with a shock front (Churazov et al. 2003).

In an earlier paper, ZuHone et al. (2010, hereafter ZMJ10), we have investigated the possible flow of heat to the cluster cool core as a result of the sloshing motions bringing hot gas from larger radii into contact with the cool gas of the core. We determined that if the gas viscosity is small, significant mixing between the hot and cold gas will result, raising the entropy of the core, though this heat flow was found to be insufficient to completely offset the effects of radiative cooling. In that work, we did not include thermal conduction, but only mixing of gases with different temperatures.

ZMJ10 modeled the ICM as a (possibly viscous) unmagnetized fluid. However, the ICM is known to be magnetized, and this fact has important considerations for heat transport. For example, it has been suggested that the narrow width of cold fronts in galaxy clusters can be used to constrain the effectiveness of heat conduction in the ICM (e.g., Ettori & Fabian 2000; Xiang et al. 2007). The sharp temperature gradient of the cold front should be smoothed out on short timescales if conduction is not suppressed across the front surface. Such suppression could be provided by magnetic field lines oriented parallel to the cold front surfaces, a situation believed to naturally arise due to the effects of magnetic draping and shear amplification (Asai et al. 2004, 2005, 2007; Lyutikov 2006; Keshet et al. 2010). In a previous paper (ZuHone et al. 2011, hereafter ZML11), we performed simulations of gas sloshing with magnetic fields, showing that such shear amplification and draping does occur around cold fronts in cluster cores. These simulations also showed that magnetic fields have similar effects as viscosity on the ICM, inhibiting the mixing of the hot and cold flows and reducing the increase of core entropy due to sloshing.

In this work, we build upon our earlier simulations (ZMJ10; ZML11) by including the effects of anisotropic thermal conduction in a simulation of gas sloshing in a cool-core galaxy cluster. We will show that the magnetic fields draped across the cold fronts indeed inhibit the smoothing out of the temperature gradient, but only to a degree. The crucial points to note are (1) due to fluid instabilities the magnetic field lines may not be oriented perfectly azimuthal with the front and (2) cold fronts are not completely enclosed surfaces—there are still field lines that connect the cold side of the cold fronts to regions of higher temperature in the cluster atmosphere over short distances. Therefore, if thermal conduction is important, it may provide a more efficient way than gas mixing to transfer heat between the hot and cold flows characteristic of sloshing motions.

This paper is organized as follows: in Section 2 we describe the simulations and the code. In Section 3 we describe the effects of anisotropic conduction and gas sloshing on the thermal state of the cluster core for different configurations of the magnetic field and various prescriptions for thermal conduction. Finally, in Section 4 we summarize our results and discuss future developments of this work. We assume a flat ΛCDM cosmology with h = 0.7, Ω_m = 0.3, and Ω_b = 0.02 h⁻².

In our simulations, we solve the following set of MHD equations (in Gaussian units):

$$\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 }\right) + \nabla {p} &=& \rho {\bf g} \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\\ &&- n_en_H\Lambda (T,Z) \end{eqnarray} \tag{ 3 }$$

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

where

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

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

where p_th is the gas pressure, is the gas internal energy per unit volume, and $\bf {B}$ is the magnetic field strength. For all our simulations, we assume an ideal equation of state with γ = 5/3 and primordial abundances with molecular weight μ = 0.6.

We performed our simulations using FLASH 3.3, a parallel hydrodynamics/N-body astrophysical simulation code developed at the FLASH Center for Computational Science at the University of Chicago (Fryxell et al. 2000; Dubey et al. 2009). FLASH 3.3 solves the equations of magnetohydrodynamics using a directionally unsplit staggered mesh algorithm (USM; Lee & Deane 2009). The USM algorithm used in FLASH 3.3 is based on a finite-volume, high-order Godunov scheme combined with a constrained transport method (CT), which guarantees that the evolved magnetic field satisfies the divergence-free condition (Evans & Hawley 1988). In our simulations, the order of the USM algorithm corresponds to the piecewise parabolic method (PPM) of Colella & Woodward (1984), which is ideally suited for capturing shocks and contact discontinuities (such as the cold fronts that appear in our simulations).

The heat flux due to thermal conduction is given as

$$\begin{equation} {\bf Q} = -\kappa _{\rm Sp}[f_\perp ({\bf I} - \hat{\textbf {\emph {b}}}\hat{\textbf {\emph {b}}})+f_\parallel \hat{\textbf {\emph {b}}}\hat{\textbf {\emph {b}}}] \cdot \nabla {T}, \end{equation} \tag{ 7 }$$

where we chose to represent the conductivity coefficient as a product of the Spitzer conductivity κ_Sp and the fractions f_∥ and f_⊥ of Spitzer conductivity parallel and perpendicular to the field line, respectively. $\hat{\textbf {\emph {b}}}\hat{\textbf {\emph {b}}}$ and ${\bf I} - \hat{\textbf {\emph {b}}}\hat{\textbf {\emph {b}}}$ are projection operators that single out the components of the temperature gradient ∇T parallel and perpendicular to the magnetic field line, with $\hat{\textbf {\emph {b}}}$ the unit vector in the direction of the local magnetic field. Most works assume f_∥ = 1 and f_⊥ = 0 in the conditions of the ICM (e.g., heat conducted strictly along the field lines without suppression), due to the large mean free path for collisions compared to the Larmor radius of the electron. We assume this prescription as the default, but also investigate alternative scenarios, represented by simulations with different values of f_∥ and f_⊥. For example, if the magnetic field is strongly curved on small scales, the effective conductivity along the field lines may be smaller than the Spitzer value, resulting in a smaller f_∥.

Alternatively, if a fraction of the magnetic field lines draped along the cold front surfaces reconnect with lines below these surfaces, a small amount of heat conduction perpendicular to the cold front surfaces may result. Our ideal MHD simulations do not model reconnection explicitly, as non-ideal terms including resistivity are not included (there is a small numerical resistivity associated with the finite resolution of the simulation, but this effective resistivity is negligible, see Section 3.2). However, we may account for the conduction that may occur across cold front surfaces in an approximate way by modeling it as a non-zero f_⊥, with the same temperature dependence as the parallel conduction. Since we are not modeling the reconnecting field lines and the conduction along them explicitly, our simple model is intended to provide only an order of magnitude estimate of this possible effect.

It is important to distinguish our "perpendicular" conduction across the cold front surfaces from conduction perpendicular to the field lines resulting from collisions between electrons. Because the Larmor radius is many orders of magnitude smaller than the electron mean free path, this perpendicular conduction is effectively zero for the conditions in galaxy clusters. Additionally, the temperature dependence of this conduction would be ∝T^−1/2 instead of the ∝T^5/2 dependence for conduction along the field lines (Braginskii 1965). In contrast, since our "perpendicular" conduction is intended to account for conduction that occurs parallel to reconnecting field lines, it has the same temperature dependence as the parallel conduction.

Following Cowie & McKee (1977), we included the effect of conduction saturation whenever the characteristic length scale associated with the temperature gradient is smaller than the mean free path, though in the bulk of the ICM this effect is not significant. We implemented anisotropic conduction following the approach of Sharma & Hammett (2007). More specifically, we applied a monotonized central (MC) limiter to the conductive fluxes. This method ensures that anisotropic conduction does not lead to negative temperatures in the presence of steep temperature gradients. A verification test of our implementation of the anisotropic conduction module in FLASH 3.3 was given in Ruszkowski et al. (2011).

We want to examine the effects of anisotropic thermal conduction in a sloshing galaxy cluster core in isolation (in particular, the effect on the cold fronts), and so for most of our simulations, we do not include the effects of radiative cooling. We then include radiative cooling to determine the effect of sloshing with anisotropic conduction on the thermal state of the cluster core. The X-ray cooling function is determined using the Mewe et al. (1995) model with a uniform metallicity of Z = 0.3 Z_☉.

The gravitational potential on the grid is set up as the sum of two "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 in Roediger & ZuHone (2012).

We refer the reader to ZML11 for the full details of our physical setup, but will provide the following short description. Our simulations consist of a massive ($M \approx 1.5 \times 10^{15} {\cal M}_\odot$, T_X ∼ 8 keV), cool-core cluster, initially in hydrostatic equilibrium, merging with a small (mass ratio R = 5) gasless subcluster, which sets off the sloshing of the cool core. We begin the simulation at a point in time when the cluster centers have a mutual separation of d = 3 Mpc and an impact parameter b = 500 kpc. The initial velocities of the subclusters are set up assuming that the total kinetic energy of the system is set to half of its total potential energy. For all of the simulations, we set up the main cluster within a cubical computational domain of width L = 2.4 Mpc on a side, with a finest cell size on our adaptive mesh refinement (AMR) grid of Δx = 2.34 kpc.

The magnetic field of the cluster is set up in an identical way to the simulations of ZML11, following the approach of Ruszkowski et al. (2007, 2008) and Ruszkowski & Oh (2010). A random magnetic field $\tilde{\bf B}({\bf k})$ is set up in $\bf {k}$-space on a uniform grid using independent normal random deviates for the real and imaginary components of the field. Thus, the components of the complex magnetic field in $\bf {k}$-space are set up such that

$$\begin{eqnarray} \tilde{B_x}(\bf {k}) &=& B_1[N(u_1) + iN(u_2)] \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \tilde{B_y}(\bf {k}) &=& B_2[N(u_3) + iN(u_4)] \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} \tilde{B_z}(\bf {k}) &=& B_3[N(u_5) + iN(u_6)], \end{eqnarray} \tag{ 10 }$$

where N(u) is a function of the uniformly distributed random variable u that returns Gaussian-distributed random values, and the values B_i are field amplitudes. We adopt a dependence of the magnetic field amplitude B(k) on the wavenumber $|\bf {k}|$ similar to (but not the same as) Ruszkowski et al. (2007) and Ruszkowski & Oh (2010):

$$\begin{equation} B(k) \propto k^{-11/6}{\rm exp}[-(k/k_0)^2]{\rm exp}[-k_1/k], \end{equation} \tag{ 11 }$$

where k₀ and k₁ control the exponential cutoff terms in the magnetic energy spectra. The cutoff at high wavenumber k₀ roughly corresponds to the coherence length of the magnetic field k₀ = 2π/λ₀ = 2π/(43 kpc), observed in some cluster cores. The cutoff at low wavenumber k₁ = 2π/λ₁ roughly corresponds to λ₁ ≈ r₅₀₀/2 ≈ 500 kpc, which is held fixed for all of our simulations. This field spectrum corresponds to a Kolmogorov shape for the energy spectrum (P_B(k)∝k^−5/3) with cutoffs at large and small linear scales. This field is then Fourier transformed to yield B(x), which is rescaled to have an average value of $\sqrt{8\pi {p}/\beta }$ to yield a field with a pressure that is proportional to the gas pressure, i.e., to have a spatially uniform β for the initial field. We set the initial β = 400. This procedure ensures that there will be significant power at all length scales between λ₁ and λ₀, which will result in a field that is significantly tangled down to the scale of λ₀.

As mentioned above, previous simulation works have predicted that due to the action of the HBI, an azimuthal magnetic field should arise in the cluster core, in the absence of turbulence driven by galaxy motions or AGNs. For this purpose, we also include a simulation that has a purely azimuthal initial field similar to that used in Bogdanović et al. (2009), described in spherical coordinates by

$$\begin{eqnarray} B_r &=& 0 \\ \ms B_\theta &=& 2B_0\sin {\theta }\cos {2\phi } \\ \ms B_\phi &=& -B_0\sin {2\phi }\sin {2\theta }. \end{eqnarray} \tag{ 12 }$$

Its magnitude was then scaled by the gas pressure with β = 400.

Finally, for both configurations, a field consistent with ∇ · B = 0 is generated by performing the Fourier transform and then by "cleaning" the field of divergence terms in Fourier space via the projection

$$\begin{equation} \tilde{\bf B}^{\prime }({\bf k}) = (\bf {I}-\bf {\hat{k}\hat{k}})\tilde{\bf B}({\bf k}), \end{equation} \tag{ 15 }$$

where $\hat{\bf {k}}$ is the unit vector in $\bf {k}$-space and I is the identity operator. Note that this operation does not change the power spectrum of the magnetic field fluctuations. Upon transformation back to real space, the magnetic field satisfies ∇ · B = 0. This field is then interpolated onto our AMR grid in such a way that the condition ∇ · B = 0 is maintained.

Because our sloshing simulations begin with the two clusters at a distance, it takes a period of time (∼2 Gyr) for the influence of the subcluster to perturb the cool core of the main cluster and initiate sloshing. If we had began the simulation with anisotropic thermal conduction "switched on," it would alter the thermal state of the core and the HBI and MTI would alter the cluster magnetic field spatial configuration long before sloshing had begun, since the growth times of these instabilities and the conduction timescale for the core are on the order of ∼ a few 10⁸ yr. Worse yet, if we began the simulation with radiative cooling switched on, the cool core would cool rapidly (also on a timescale of ∼ a few 10⁸ yr) and we would have a cooling catastrophe before sloshing started. To avoid these effects, we switch on conduction (and cooling in the simulations where it is included) at a time t ∼ 2 Gyr after the beginning of the simulation, shortly after the subcluster has made its pericentric passage of the cluster core and right at the beginning of the sloshing period. This enables us to examine the effects of conduction and sloshing together on a typical cool-core cluster state.

An unfortunate side effect of this evolutionary period before the core passage is that the magnetic field configuration will slowly evolve away from its initial configuration. Hence, the initially isotropically tangled field configuration with $\langle |\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}}|\rangle = 0.5$ evolves to a slightly more radial field configuration with $\langle |\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}}|\rangle \approx 0.56$, and the initially azimuthal field configuration with $\langle |\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}}|\rangle = 0$ also evolves to a slightly more radial field configuration with $\langle |\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}}|\rangle \approx 0.35$. This is not a major concern, as the former is not too far off from an isotropically tangled field configuration and the latter is not too far from the field configuration of the saturation state of the HBI as seen in simulations (Parrish et al. 2009; Bogdanović et al. 2009). In real clusters, there is always turbulence and small disturbances from galaxies and DM clumps that keep the field tangled (Parrish et al. 2010; Ruszkowski & Oh 2010, 2011).

We will refer to the simulations using the following codes. The letter "R" will indicate a cluster that is evolved in isolation, whereas "S" will indicate a cluster evolved with a gravitational interaction causing sloshing. "C" will indicate anisotropic thermal conduction is turned on, and "X" will indicate X-ray radiative cooling is turned on. Simulations that differ in the prescription for conduction and the initial field line configuration will be delimited by numbers. Table 1 lists all simulation runs with their particular parameters.

3.1. Brief Outline of the Sloshing Process

3.2. The Effect of Conduction on Cold Fronts

Figures 1–4 show slices through the center of the cluster of the gas temperature for several different epochs, for our simulations with and without anisotropic conduction (but without radiative cooling). We will first compare runs with no conduction to runs with full Spitzer conduction along the field lines (f_∥ = 1, f_⊥ = 0). Our runs S1 (top left panels) and S2 (top center panels) correspond to sloshing simulations with initially tangled fields and initially azimuthal fields, respectively. Likewise, our runs SC1 (top right panels) and SC2 (bottom left panels) correspond to sloshing simulations with anisotropic heat conduction, with initially tangled fields and initially azimuthal fields.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Including conduction has a dramatic effect on the temperature structure of the cold fronts. The adiabatic simulations without conduction (S1 and S2) develop a number of well-defined cold fronts over time, which maintain temperature differences of ∼1.5–2 across the fronts throughout their evolution. This is not the case for the simulations with anisotropic unsuppressed Spitzer conduction (simulations SC1 and SC2). Over time, the cold gas inside the fronts is heated up, so that the temperature contrasts across the fronts become smaller, roughly ∼1.3–1.5. The changing thermal state of the gas also changes the shape of the fronts in some places. In particular, due to the fact that the inner 50 kpc becomes nearly isothermal, cold fronts do not occur on the smallest scales, in conflict with observations of cool-core clusters, which typically show a number of cold fronts at both small and intermediate radii within the core. The temperature of the core in simulation SC2 increases more slowly than in simulation SC1, because the more azimuthally oriented field lines in the former suppress conduction from larger radii more efficiently than do the more tangled fields of simulation SC1.

Figures 1–4 also show results for the simulations with other values for the conduction coefficients f_⊥ and f_∥. The SC3 simulation has the same field line configuration as the SC1 simulation, but the conduction along the lines is suppressed by a factor of 10. If field lines are strongly curved on very small linear scales, the effective conduction coefficient parallel to the field lines may be smaller than the Spitzer value by an order of magnitude due to the presence of stochastic magnetic mirrors (e.g., Narayan & Medvedev 2001). We choose f_∥ = 0.1 to determine the effect of such a reduction. Alternatively, if plasma microscale instabilities play a significant role, the thermal conductivity of the ICM may be many orders of magnitude smaller (Schekochihin et al. 2008). For all practical purposes, this case is equivalent to the runs without conduction.

On the other hand, if (for example) ∼1% of the magnetic field lines parallel to the fronts reconnect, the effective fraction of Spitzer conduction perpendicular to the large-scale direction of the field could be f_⊥ ∼ 0.01. We assume this value of f_⊥ for simulation SC4, keeping the value of f_∥ = 0.1 for comparison with simulation SC3. As expected, for the SC3 simulation, far less heat is conducted to the cluster core, and the temperature contrast of the cold fronts with the surrounding medium is much closer to the case without any conduction. In the SC4 simulation, the small perpendicular conductivity, combined with a steep temperature jump across the cold front, has succeeded over time to smooth out the temperature gradient across the cold fronts.

Figure 5 compares temperature gradients across cold fronts between the S1 and the SC1 simulations at a late epoch, t = 3.25 Gyr. Four particular locations across the fronts are picked out for inspection (see Figure 3 for the profile locations, marked by the corresponding letters). The profiles show that the temperature jumps in the conducting simulation have all been reduced. For some of the profiles (e.g., (a) and possibly (d)), the width of the temperature jump in the simulation with conduction SC1 is not greater than the jump in the non-conducting simulation S1 (a few simulation cells, ∼5–10 kpc—essentially unresolved by the simulations), whereas for other profiles in simulation SC1, the temperature jumps have been smoothed out over several zone sizes (we will explore these issues further in Section 3.3). Particularly apparent in profile (d) is that conduction has modified the temperature distribution below the front, while preserving the sharp jump across the front.

Zoom In Zoom Out Reset image size

Figure 5 also highlights example temperature jumps in the SC3 and SC4 simulations (along the same profile locations, see Figure 3). The temperature difference across the fronts are similar between these two simulations and the S1 simulation, due to the high suppression of conduction, but the smoothing out of the gradient due to the conduction perpendicular to the field lines in the SC4 simulation is noticeable, with the temperature change occurring over roughly twice the number of cells (∼10 kpc) compared to the case without perpendicular conduction due to field line reconnection.

The jumps in density across the cold fronts are similarly affected by conduction. As heat flows from the hot gas above the cold fronts to the cool gas below the front, the former contracts and the latter expands in the same proportion as the change in temperature. Figure 6 shows the density jumps across the cold fronts at the same positions in Figure 3. In simulations with conduction, the density jumps are reduced and the widths of the density jumps are increased to varying degrees, in the same proportion as the temperature jumps except that the sign of the gradient is reversed. In the most extreme cases, the density jumps have been essentially erased, either in the strong parallel conduction case SC1 or the case SC4 with a small component of conduction perpendicular to the magnetic field lines. The smoothed profiles in these simulations are in disagreement with observations of cold fronts in clusters.

Zoom In Zoom Out Reset image size

The differences between the simulations S1, SC3, and SC4 indicate that numerical reconnection across the cold front surfaces and any associated heat flux due to this reconnection is negligible. The cold fronts in the simulations S1 and SC3 have essentially the same temperature and density jumps, whereas adding a small perpendicular conduction smooths out the jumps over twice the number of cells. The fact that we do not see this degree of smoothing in the SC3 simulation indicates that numerical reconnection is too small to provide any significant heat flux across the fronts.

The smaller density jumps and the modification of the density profiles below the fronts that reduce the density contrast have a significant effect on the appearance of cold fronts as seen in X-ray images, due to the strong dependence of the X-ray emission on the gas density. Figure 7 shows synthetic X-ray images for the S1, SC1, and SC3 simulations, at the epoch t = 2.75 Gyr, created by computing the bolometric MeKaL emissivity in each FLASH cell and projecting it along the z-axis (in the same manner as ZML11). Additionally, a Chandra X-ray image of A2319 is also presented for comparison (this cluster is slightly hotter than our model cluster, with average T ∼ 9–10 keV). This cluster exhibits a cold front with a very sharp jump in X-ray emission, with a width roughly the size of the Chandra point-spread function (O'Hara et al. 2004; Govoni et al. 2004).

Zoom In Zoom Out Reset image size

The characteristic spiral shape of the cold fronts is apparent in each simulated image. However, the sharp edges in emission across the front surfaces are essentially absent from the SC1 case, whereas in the SC3 and S1 simulations, the sharp jumps are clearly seen. Therefore, for our simulated hot cluster, unsuppressed Spitzer conduction appears to produce cold fronts that are far less prominent than those in observed hot clusters, and it is questionable whether or not they may be even called "cold fronts" at all. We reserve a more quantitative comparison with observed clusters of different temperatures and plasma β values and the possible constraints on conductivity for a future paper.

It also is apparent from these figures that conduction may have an effect on the smoothness of the front surfaces. For example, the cold fronts in simulations SC1 and SC2 exhibit fewer small-scale features along the front surfaces than their counterparts without conduction, S1 and S2. The widening of the cold front interfaces (seen clearly in Figure 6) will cause them to be less susceptible to the Kelvin–Helmholtz (K-H) instability (see, e.g., Churazov & Inogamov 2004).

3.3. The Geometry of the Magnetic Field Lines

The simulations from ZML11 showed magnetic field is amplified and the lines oriented parallel to the cold front surfaces, which would prevent conduction across the cold front surfaces from the hotter gas just above the fronts. The fact that there are examples of discernable temperature jumps over the same spatial range (a few simulation cells) even in the SC1 run indicates that conduction is still strongly suppressed directly across the front surfaces. What accounts for the significant changes in density and temperature across the fronts that still occur?

The cool gas below the cold fronts is only completely immune to the effects of heat conduction if there are no field lines leading from hotter regions of the cluster atmosphere into the cold regions of the fronts from other directions. The fronts seen in X-ray observations do not completely enclose the cool core. In principle, there could be magnetic field lines that reach the regions below the fronts from hotter regions on the other side of the cluster. More importantly, simulations of gas sloshing show that the resulting spiral structure has an alternating pattern of hot and cold flows, with hot flows developing below cold fronts as they expand to higher radii. Keshet (2012) showed that these latter flows are of hot gas being drawn into the core region form higher radii. These regions can be magnetically connected to the cool gas immediately below the fronts. Additionally, though the cold fronts that form in simulations are extended in directions perpendicular to the sloshing plane, the characteristic size of the fronts in this direction (along the line of sight where the observer sees the spiral structure) is smaller than in the plane. This exposes the cold regions of the fronts to magnetic field lines extending into hotter regions along directions perpendicular to the sloshing plane.

The preceding arguments lead us to expect that some field lines that pass through cold gas below the front surfaces may be connected to hotter regions over relatively short distances. Figures 8 and 9 show magnetic field lines beginning at several points in the cold region below the front surface and followed for a distance of 100 kpc. A large number of field lines are tangentially oriented to the front surfaces, as can be easily seen in both projections, shielding this cold gas from the hotter gas directly above the front. However, the figures show that there are still a number of field lines that extend into hotter regions of the cluster core along other directions. At the epochs of t = 2.5 Gyr and t = 3.0 Gyr, the coldest temperatures below the front surfaces are approximately T ∼ 4.5–5 keV, and these regions are connected to regions with temperatures of T ∼ 6–7.5 keV. In the right panel of Figure 8, field lines stretched along the y and z directions connect cold gas underneath the front surfaces to hotter gas. In Figure 9, lines connecting cool gas with hot gas that was brought in by sloshing are evident near the center of the left panel of the figure.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the simulations show that there is not always a perfect alignment of strong magnetic field layers with cold fronts. Figure 10 shows slices through the center of the simulation domain in the x–y plane (coincident with the plane of the sloshing motions) of the gas temperature and the magnetic field energy for the epochs t = 2.75 Gyr and t = 3.25 Gyr with the magnetic field vector components in the x- and y-directions overlaid. Extended regions of strong magnetic field indicate places where the magnetic field strength has been increased by shear amplification. The locations of the profiles from Figure 3 are also marked. While some locations along the cold fronts exhibit a strong field layer parallel to the front, this is not strictly true everywhere along the fronts. In particular, the t = 2.75 Gyr panel shows that the trailing edge of the cold spiral in the southwestern direction (indicated roughly by the white rectangles in Figure 10) is not draped by a high-B layer. It can be seen that the magnetic field vectors along this portion of the spiral front are not strictly aligned parallel to the front surface. Regions of ordered, strong field are interspersed with regions of weak, tangled field as shown by the vectors at these locations. At a somewhat later time, t = 3.25 Gyr, this location still lacks a strongly magnetized parallel layer. This also happens to be the location of profile (d). That particular temperature jump was shown to be smoothed out significantly in Figure 5. Though at the epoch t = 3.25 Gyr there does appear to be a small, parallel, highly magnetized layer at that particular front location, this layer was not present previously and the temperature jump likely has been smoothed out before the formation of this layer.

Zoom In Zoom Out Reset image size

Cold fronts that are not draped by strong magnetic fields will be more susceptible to the K-H instability, which will distort the front surface and lead to gas mixing between hot and cold phases. This is because the orientation of the draping fields tends to coincide with the direction of the shearing motions. Perturbations of the field lines perpendicular to the front surfaces at these locations will lead to smoothing out of the fronts as outlined in Lecoanet et al. (2012), though the extent of this smoothing will be limited by the HBI, which will reorient magnetic field lines parallel to the front surface. Lecoanet et al. (2012) estimated that the scale height of the temperature gradient for such a front would smooth to a height of h_T ∼ χ^2/3/g^1/3, where χ is the Spitzer thermal diffusivity and g is the gravitational acceleration. Taking example values from our simulation for a cold front at r ≈ 120 kpc, we find χ ∼ 5 × 10³⁰ cm² s⁻¹ and g ∼ 3 × 10⁻⁸ cm s⁻². This yields an approximate scale height of the temperature gradient across the front of h_T ∼ 30 kpc, which is similar to the widths of some of the temperature gradients across the fronts in Figure 5.

In simulations of thermal instabilities in magnetized galaxy clusters, a helpful quantity to examine has been the average angle between the magnetic field and radial unit vectors, given by

$$\begin{equation} \langle \theta _B\rangle = \cos ^{-1}\langle |\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}}|\rangle. \end{equation} \tag{ 16 }$$

Isotropically oriented tangled field lines will have an average $\langle |\hat{\textbf {\emph {b}}}\cdot \hat{\textbf {\emph {{r}}}}|\rangle = 0.5$ and hence the corresponding angle is 〈θ_B〉 ≈ 60°. The HBI drives field lines in a more azimuthal direction (θ_B ∼ 70°–80°) in the presence of a positive radial temperature gradient (as in a cool core). Sloshing motions in general preserve the overall positive temperature gradient of the cool core, implying that the HBI should still be active. However, they will also drag the field lines around into a more azimuthally oriented configuration as they are draped around the cold fronts.

Figure 11 shows the evolution of the volume-averaged angle of the magnetic field with respect to the radial direction within a radius of 100 kpc over time for the S1, S2, SC1, and SC2 simulations. The simulations with initially tangled and azimuthal fields begin with 〈θ_B〉 ≈ 56° and 〈θ_B〉 ≈ 70°, respectively. We find that largely the evolution of this "magnetic angle" is driven by the sloshing motions themselves, since the evolution of θ_B is essentially identical between the simulations with and without conduction that have the same initial magnetic field setup. An exception to this general rule are the simulations with initially azimuthal magnetic field lines, S2 and SC2. When conduction is included, the final θ_B indicates a more azimuthal direction than in the case where conduction is absent. This is possibly due to the fact that in this case the core is more magnetically isolated due to the azimuthally oriented field lines, and the positive temperature gradient is better preserved in the core than in the tangled field line case. In such a scenario, the HBI would drive the field lines in a more azimuthal direction once the conduction is turned on.

Zoom In Zoom Out Reset image size

3.4. The Evolution of the Entropy of the Cluster Core

Given that our simulations show that the magnetic field geometry of the sloshing cool core still permits heat conduction to the cluster core, it is important to determine what effect this has on the thermal state of the core. The evolution of average entropy within the inner 50 kpc of the cluster center over time is shown in Figure 12. In the simulations without conduction (S1 and S2), the average entropy of the cluster core increases due to sloshing bringing hot gas from the cluster outskirts with the cool gas of the core into contact, which then mix to a certain degree. However, the magnetic field suppresses such mixing (ZML11), and the increase in entropy levels off roughly 0.5 Gyr after the simulation begins. ZML11 showed that the difference in initial field geometry has little effect on the entropy increase in this case.

Zoom In Zoom Out Reset image size

In contrast, in the simulations with anisotropic Spitzer conduction (SC1 and SC2), heat conducts along the field lines to the cluster center, and the entropy increase of the core continues unabated for the duration of the simulation. As expected, the effect is stronger in the case with initially random field lines than in the case with initially azimuthal field lines. The average specific entropy increase in these two simulations is roughly ΔS ∼ 20–30 keV cm² higher than in the adiabatic simulations.

A similar situation exists for the simulations with weaker conduction coefficients. The simulation with conductivity that is 1/10 of the Spitzer value strictly along the field lines (SC3) has a correspondingly smaller increase in entropy over the SC1 simulation, with ΔS ∼ 5 keV cm², close to the entropy increase due to sloshing alone, but the entropy is still increasing at the end of the simulation. The inclusion of even the small, 1/100 Spitzer conduction perpendicular to the field lines (SC4) results in a significantly higher entropy increase over simulation S1 than the SC3 case, with ΔS ∼ 10 keV cm²—this is because in the sloshing core, temperature gradients across the magnetic field lines are much higher than along the lines. The result of all of the simulations with some form of anisotropic conduction is that heat can be conducted efficiently to the core, even in the presence of azimuthal cold fronts.

3.5. Simulations with Conduction and Radiative Cooling

We now include radiative cooling in our simulations. The runs with radiative cooling SCX1 and SCX2 correspond to the non-cooling runs SC1 and SC2. The cooling is "turned on" at the same moment as the anisotropic conduction, at the epoch t = 2 Gyr. As control cases to isolate the effects of sloshing and conduction, we also ran the following simulations. RX is an initially relaxed cluster left in isolation without heat conduction. Due to the absence of any heating mechanism, it experiences a cooling catastrophe in the core within a few hundred Myr. RCX is a similar run with an isolated cluster, except it also includes Spitzer conduction along the field lines. Thermal conduction partially offsets the cooling of the core in this run, but the HBI quickly reorients the magnetic field lines in the core in the azimuthal direction and a cooling catastrophe still occurs. The last run, SX, is identical to simulation S1 (without conduction and with sloshing), but with radiative cooling included.

In ZMJ10, we showed that radiative cooling may be slowed by sloshing itself, due to the mixing of hot gas from the outskirts with the cold gas for the core. However, the inclusion of even weak magnetic fields reduces this mixing and inhibits the consequent heating of the core, depending on the strength of the magnetic field (ZML11). From our results in ZMJ10, we expect that regardless of the strength of the magnetic field, the heat delivered to the core from sloshing alone will not be sufficient to prevent a cooling catastrophe at the innermost regions of the cool core. The left panels of Figure 13 show slices of temperature through the center of the cluster for selected epochs for the SX simulation. In this run, the core temperature drops below 1 keV within a radius of r ∼ 20 kpc. Additionally, the core density peaks at high values of n_e ∼ 1–5 cm⁻³, in conflict with observations of real clusters. Therefore, the goal of the runs with conduction and cooling is to determine if sloshing with conduction together is able to prevent such a cooling catastrophe.

Zoom In Zoom Out Reset image size

The center and right panels of Figure 13 shows slices of temperature through the center of the cluster for selected epochs for the sloshing simulations including cooling and conduction. Anisotropic conduction, even in the case where it is most efficient, is unable on its own to suppress the runaway cooling in the very center of the cluster, with central temperatures falling to the temperature floor of our simulation (T ∼ 0.01 keV), far below the temperatures observed in real clusters. However, the temperature is significantly raised outside the very central region of the cool core, at radii r ≳ 10 kpc.

The two leftmost panels of Figure 14 show the evolution of the gas mass with short cooling times within the inner 50 kpc of the cluster center over time. We were only able to run the RX simulation out to an epoch of t ∼ 3.6 Gyr, as the extremely high densities, low temperatures, and magnetic pressures in this simulation made the simulation too unstable to continue past this point. Finally, we have also included a cooling run SXNoB, identical to simulation SX, but without magnetic fields, to compare the rate of increase of cool gas when the mixing of hot and cool gases is not suppressed by such fields. This the same simulation as in ZMJ10, but the setup (resolution and model for the gravitational potential) is identical to our present work, to enable a direct comparison.

Zoom In Zoom Out Reset image size

The top panels of Figure 14 show the evolution of the gas mass with short cooling times within the inner 50 kpc of the cluster center over time. We were only able to run the RX simulation out to an epoch of t ∼ 3.6 Gyr, as the extremely high densities, low temperatures, and magnetic pressures in this simulation made the simulation too unstable to continue past this point. Finally, we have also included a cooling run SXNoB, identical to simulation SX, but without magnetic fields, to compare the rate of increase of cool gas when the mixing of hot and cool gases is not suppressed by such fields. This the same simulation as in ZMJ10, but the setup (resolution and model for the gravitational potential) is identical to our present work, to enable a direct comparison.

Over the course of ∼2 Gyr, all of the simulations have a significant increase of the mass of cool gas in their cores. However, both sloshing and conduction reduce the buildup of cool gas. The top left panel of Figure 14 shows the buildup of gas mass with t_cool < 0.1 Gyr. At the onset of cooling, there is no gas in the core with a cooling time this short. Within ∼0.25 Gyr, in the relaxed cluster simulations (RX and RCX), a gas mass of ∼10⁹–10¹⁰ M_☉ with such cooling time has built up in the core. This is somewhat slowed by conduction in the RCX simulation, and by t ∼ 4 Gyr the mass of cool gas is lower than that of the RX simulation by a factor of a few. In the simulations with sloshing (SX and SCX1), the initial buildup of cool gas does not begin until ∼0.5 Gyr after the onset of cooling. However, in the SX simulation, the mass of cool gas eventually catches up with that of the RX simulation. Out of all of these simulations, the run with both sloshing and conduction (SCX1) is most effective at suppressing the cooling of the gas for the long term, with nearly an order of magnitude less mass of cool gas in the core than in the RX simulation.

In the case of simulation SXNoB, no gas below t_cool = 0.1 Gyr appears until ∼1 Gyr after the onset of cooling. Gas mixing, which is much more efficient in the absence of a magnetic field (or perhaps in the presence of field variations on microscopic scales that would strongly reduce field tension), provides a source of heat to the core gas and prevents the catastrophic cooling of gas until this later stage. However, the mass of cool gas increases quickly after this point and is eventually comparable with that of simulation RCX. This indicates that in the absence of magnetic field tension, strong gas mixing due to sloshing is a comparable source of heat to anisotropic conduction.

The top right panel of Figure 14 shows the buildup of gas mass with t_cool < 1.0 Gyr. In this case, all simulations already start with the same mass of gas with cooling timescales this short. However, the general trend for each of these simulations is the same as for the case for the gas mass with t_cool < 0.1 Gyr. The runs with conduction (or gas mixing) and sloshing exhibit a slower buildup of cool gas.

The bottom panels of Figure 14 show slices through the z = 0 plane of the temperature at the epoch t = 3.5 Gyr, with contours marking the surfaces with specific cooling times of 0.1, 0.5, and 1.0 Gyr, for the simulations RX and SCX1. It is clear that the effect of sloshing and conduction together is to reduce the volume of the region with very low cooling times and temperatures. Without conduction and sloshing, a r ∼ 10 kpc region with t_cool < 0.1 Gyr develops, whereas with both of these effects included, this region is restricted to r ∼ 3 kpc.

Overall, this result indicates that sloshing, even with anisotropic thermal conduction, is not enough to completely stabilize a cool core against a cooling catastrophe, and that for the central regions of clusters, AGN feedback, or other mechanisms are necessary to regulate the temperature of the core. However, sloshing and conduction are able to suppress or significantly delay catastrophic cooling in the cluster cooling cores (r ∼ 50–100 kpc) outside the very central regions of radii r ≲ 3 kpc.

Gas sloshing, a prevalent phenomenon in relaxed galaxy clusters, is evidenced by the nearly ubiquitous presence of spiral-shaped cold fronts in cluster cool cores. The temperature jumps across these fronts imply that either thermal conduction is intrinsically weak in galaxy clusters or the magnetic field is oriented parallel to the front surfaces, restricting conduction across the fronts. Most previous works have assumed the latter, and simulations have confirmed that such magnetic layers form around cold fronts in clusters.

Sloshing may also have an important role in determining the thermal evolution of the cluster core. It brings the cold gas of the core into contact with hot gas, possibly facilitating a transfer of heat between these phases, either via mixing of the hot and cold gas or by heat conduction. This raises the possibility that sloshing may be partially responsible for preventing a "cooling catastrophe" in the cores of galaxy clusters. However, the ICM is magnetized, which, as we showed in ZML11, partially suppresses mixing.

To determine whether or not anisotropic heat conduction in the presence of magnetic fields has any effect on sloshing cold fronts and the thermal state of the cool cores, we have performed MHD simulations with various configurations of the magnetic field and the conduction coefficients. We determine that, in spite of the fact that the cold fronts indeed are mostly draped by a strong parallel magnetic field, they are not completely thermally isolated from hotter gas. The shear amplification and stretching of magnetic field lines occurs only in a thin layer at some of the cold front surfaces, which only prevents heat conduction to the cold gas below the front from the hot gas directly above the front. Heat can still flow around the cold front, along field lines connected to the cold gas from other regions, either from the hot gas at higher radii perpendicular to the sloshing plane or from the hot flows that are brought inside the cool core by sloshing. Thus, while the sharpness of the temperature gradient at the cold front surfaces may be preserved at some locations along the fronts, the cool gas below the fronts can still be heated from other directions. Also, at some points along the cold front surfaces, the sharpness of the temperature gradient is not preserved because a stabilizing parallel magnetic field layer fails to form. Conduction also leads to the reduction of the density jumps and the widening of the cold front interfaces, which can aid in suppressing K-H instabilities. Importantly, for the hot simulated cluster presented in this work, we find that if conduction is unsuppressed along the field lines, the density jumps can be smoothed to such a degree that the cold fronts are no longer apparent in synthetic X-ray observations, in profound disagreement with observed clusters. Due to the strong temperature dependence of Spitzer conduction, we cannot use our model cluster, which is hot (T_X ∼ 8 keV), to draw conclusions about the smearing out of cold fronts due to conduction in cooler clusters. In a future paper, we will attempt more quantitative comparisons with observed cold fronts in other clusters in order to more generally constrain conduction.

We find that a conduction coefficient parallel to the field lines an order of magnitude smaller than the Spitzer value has a nearly insignificant effect on the temperature of the core and the density and temperature jumps at the cold front surfaces, producing cold fronts that are in agreement with observations. On the other hand, adding even a small perpendicular conduction (0.01 of Spitzer), that may, for example, mimic the effect of field line reconnection, has an effect on the core specific entropy similar to the case with unsuppressed Spitzer conduction along the field lines (because the temperature gradients across the field lines are higher), and slightly smoothes out the temperature gradients of the fronts over the course of a few Gyr.

In the absence of radiative cooling, gas mixing and the conduction of heat between the hot and cold phases in the cluster core raises the temperature and the entropy of the core, and if left unchecked, would eliminate the cool core entirely, leaving behind a high-entropy isothermal core. When radiative cooling is taken into account, we find that, even for full Spitzer conduction along the field lines, a cooling catastrophe still occurs in the cluster center, though the effects of conduction and sloshing together can significantly suppress the buildup of cool gas in the core and restrict the cooling catastrophe to a small (r ≲ 5 kpc) volume in the very center. The fact that conduction is unable to suppress cooling completely in the cluster center is somewhat inevitable, as for isobaric perturbations the cooling rate increases with lower temperature and the heating rate from conduction decreases (Kunz et al. 2011). AGN feedback or other mechanisms are still necessary to prevent gas from cooling in these central regions. However, our results show that sloshing and conduction together can produce a significant amount of heat to offset some cooling in the cluster core.

We note that we have only considered one merger scenario and one initial β for the magnetic field. Though different choices in these parameter spaces are not expected to change our main conclusions, from our previous works we may make reasonable guesses about the effect of changing our initial conditions. Certain merger setups result in more mixing of hot and cold gases than others (ZMJ10), whereas a lower initial β results in stronger magnetic field layers at the cold front surfaces and hence more suppression of K-H instabilities (ZML11). Most importantly for our considerations, the strong temperature dependence of Spitzer conduction makes it particularly strong for our T ∼ 8 keV simulated cluster. In order to use observations of cold fronts in real clusters to put constraints on suppression of conduction, a wider parameter space of initial cluster temperatures and plasma β values would need to be explored. Such a study, including a quantitative comparison with a sample of observed cold fronts, will be presented in a future work.

J.A.Z. thanks Ian Parrish and Mikhail Medvedev for useful discussions and advice. Calculations were performed using the computational resources of the National Institute for Computational Sciences at the University of Tennessee and the Advanced Supercomputing Division at NASA/Ames Research Center. Analysis of the simulation data was carried out using the AMR analysis and visualization toolset yt (Turk et al. 2011), which is available for download at http://yt-project.org. J.A.Z. is supported by the NASA Postdoctoral Program. M.R. acknowledges NSF grant 1008454. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago.