COSMOLOGICAL SIMULATIONS OF ISOTROPIC CONDUCTION IN GALAXY CLUSTERS

The Enzo code couples an N-body particle-mesh solver (Efstathiou et al. 1985; Hockney & Eastwood 1988) used to follow the evolution of a collisionless dark matter component with an Eulerian AMR method for ideal gas dynamics by Berger & Colella (1989), which allows high dynamic range in gravitational physics and hydrodynamics in an expanding universe. This AMR method (referred to as structured AMR) utilizes an adaptive hierarchy of grid patches at varying levels of resolution. Each rectangular grid patch (referred to as a "grid") covers some region of space in its parent grid which requires higher resolution, and can itself become the parent grid to an even more highly resolved child grid. Enzo's implementation of structured AMR places no fundamental restrictions on the number of grids at a given level of refinement, or on the number of levels of refinement. However, owing to limited computational resources, it is practical to institute a maximum level of refinement, ℓ_max. Additionally, the Enzo AMR implementation allows arbitrary integer ratios of parent and child grid resolution, though in general for cosmological simulations (including the work described in this paper) a refinement ratio of two is used.

Multiple hydrodynamic methods are implemented in Enzo, the primary two being the piecewise parabolic method (PPM) of Colella & Woodward (1984) and the method from the ZEUS code (Stone & Norman 1992a, 1992b). PPM is more accurate than the ZEUS method, but is prone to instability in conditions with steep internal energy gradients, which are common in simulations of cosmological structure formation when radiative cooling is used. In these situations, PPM will occasionally predict negative internal energies, causing the simulation to crash. For this reason, we elect to use the ZEUS method for all the simulations presented here.

We use the metallicity-dependent radiative cooling method described in Smith et al. (2008) and Smith et al. (2011). This method solves the non-equilibrium chemistry and cooling for atomic H and He (Abel et al. 1997; Anninos et al. 1997) and calculates the cooling and heating from metals by interpolating from multi-dimensional tables created with the photo-ionization software, Cloudy⁷ (Ferland et al. 1998). Both the primordial and metal cooling also take into account heating from a time-dependent, spatially uniform metagalactic UV background (Haardt & Madau 2001).

To model the effects of star formation and feedback, we use a modified version of the algorithm presented by Cen & Ostriker (1992). A grid cell is capable of forming stars if the following criteria are met: the baryon overdensity is above some threshold (here 1000), the velocity divergence is negative (i.e., the flow is converging), the gas mass in the cell is greater than the Jeans mass, and the cooling time is less than the self-gravitational dynamical time. Because conduction can potentially balance radiative cooling, instead of the classical definition of the cooling time, $e/\dot{e}$, we use a modified cooling time that includes the change in energy from conduction, defined as

$$\begin{equation} t_{{\rm cool}} = \frac{e}{\dot{e}_{{\rm rad}} + \dot{e}_{{\rm cond}}}, \end{equation} \tag{ 1 }$$

where $\dot{e}_{{\rm rad}}$ is the cooling rate and $\dot{e}_{{\rm cond}}$ is the conduction rate. In the original implementation, if all of the above criteria are satisfied, then a star particle representing a large, coeval group of stars with the following mass is formed:

$$\begin{equation} m_{*} = f_{*} \ m_{\rm cell} \ \frac{\Delta t}{t_{\rm dyn}}, \end{equation} \tag{ 2 }$$

where f_* is an efficiency parameter, m_cell is the baryon mass in the cell, t_dyn is the dynamical time of the combined baryon and dark matter fluid, and Δt is the timestep taken by the grid containing the cell in question. The factor of t_dyn is added to provide a connection between the physical conditions of the gas and the timescale over which star formation and feedback occur. When a star particle is created, feedback in the form of thermal energy, gas, and metals is returned to the grid at a rate given by

$$\begin{equation} \Delta m_{\rm sf} = m_{*} \ \frac{\Delta t}{t_{\rm dyn}} \ \frac{(t - t_{*})}{t_{\rm dyn}} \ e^{-(t - t_{*}) / t_{\rm dyn}}, \end{equation} \tag{ 3 }$$

where t_* is the creation time for the particle, so that Δm_sf rises linearly over one dynamical time, then falls off exponentially after that. This feature ensures that the feedback response unfolds over a dynamical time, even if star formation in the simulation is formally instantaneous. We also use a distributed feedback model, in which feedback material is distributed evenly over a 3 × 3 × 3 cube of cells centered on the star particle's location. This was shown by Smith et al. (2011) to transport hot, metal-enriched gas out of galaxies more effectively, and to minimize overcooling issues that are common to simulations of structure formation. This method was also shown by Skory et al. (2013) to provide a better match to observed cluster properties than injecting all feedback into a single grid cell.

In practice, Enzo simulations often implement an additional restriction, allowing the creation of a star particle only if m_* is above some minimum mass. This prevents the formation of a large number of low-mass star particles whose presence can significantly slow down a simulation. However, the use of this algorithm in concert with conduction must be done with care. As we discuss in Section 2.2, explicit modeling of thermal conduction can require timesteps that are significantly shorter than the hydrodynamic Courant condition, which has the effect of considerably reducing the value of m_*. If the value of m_* is only slightly less than the minimum particle size, this condition will likely be satisfied in a relatively short time as condensation proceeds, causing the gas density to increase and its dynamical time to decrease. When the conduction algorithm is active, the resulting time delay can allow cold, dense gas to persist while in thermal contact with the hot ICM. Since radiative cooling is proportional to the square of the density, heat conducted into high-density gas can be easily radiated away. This artificial time delay therefore produces a spurious heat sink in the ICM that can potentially boost the rates of cooling and star formation. In order to avoid this unphysical behavior, we remove the factor of Δt/t_dyn from Equation (2) and adopt a constant value of 10 Myr for t_dyn for use in Equation (3). We have tested the effects of these changes by comparing two simulations run without conduction, using both the original star-formation implementation and this modification, and find the stellar masses at any given time differ by less than 1%.

We implement the equations of isotropic heat conduction in a manner similar to that of Jubelgas et al. (2004), where the heat flux, j, for a temperature field, T, is given by

$$\begin{equation} j = -\kappa \nabla T, \end{equation} \tag{ 4 }$$

where κ is the conductivity coefficient. In an ionized plasma, heat transport is mediated by Coulomb interactions between electrons. In this scenario, referred to as Spitzer conduction (Spitzer 1962), the conductivity coefficient is given by

$$\begin{equation} \kappa _{{\rm sp}} = 1.31 n_{e} \lambda _{e} k \left(\frac{k_BT_{e}}{m_{e}} \right)^{1/2}, \end{equation} \tag{ 5 }$$

where m_e, n_e, λ_e, and T_e are the electron mass, number density, mean free path, and temperature, and k_B is the Boltzmann constant. The electron mean free path is

$$\begin{equation} \lambda _{e} = \frac{3^{3/2}(k_BT)^{2}}{4\pi ^{1/2}e^{4}n_{e} \ln \Lambda }, \end{equation} \tag{ 6 }$$

where e is the electron charge and ln Λ is the Coulomb logarithm, defined to be the ratio of the maximum and minimum impact parameters over which Coulomb collisions yield significant momentum change in the interacting particles, which are electrons in this case. The Coulomb logarithm for electron–electron collisions is

$$\begin{equation} \ln \Lambda = 23.5 - \ln \ \left(n_{e}^{1/2} T_{e}^{-5/4}\right) - \left(\frac{10^{-5} + (\ln \ (T_{e}) - 2)^{2}}{16} \right)^{1/2} \end{equation} \tag{ 7 }$$

(Huba 2011, p. 34). Because of its relative insensitivity to n_e and T_e, we follow Jubelgas et al. (2004) and Sarazin (1988) and adopt a constant value of ln Λ = 37.8, corresponding to n_e = 1 cm⁻³ and T_e = 10⁶ K. For values more appropriate to the ICM (n_e = 10⁻³ cm⁻³, T_e = 10⁷ K), ln Λ = 43.6. At a constant density, ln Λ varies by ∼10%–20% over the range of temperatures relevant here. The Spitzer conductivity then reduces to

$$\begin{equation} \kappa _{{\rm sp}} = 4.9 \times 10^{-7} T^{5/2}\ {\rm erg\ s^{-1}\ cm^{-1}\ K^{-1}}. \end{equation} \tag{ 8 }$$

Since ln Λ increases with T and κ_sp∝1/ln Λ, including a direct calculation of ln Λ would serve to somewhat soften dependence of κ_sp on T. In low-density plasmas with large temperature gradients, the characteristic length scale of the temperature gradient, ℓ_T ≡ T/|∇T|, can be smaller than the electron mean free path, at which point Equation (5) no longer applies. Instead, the maximum allowable heat flux is described by a saturation term (Cowie & McKee 1977), given as

$$\begin{equation} j_{{\rm sat}} \simeq 0.4 n_e k_B T \left(\frac{2 k T}{\pi m_e} \right)^{1/2}. \end{equation} \tag{ 9 }$$

To smoothly connect the saturated and unsaturated regimes, we use an effective conductivity (Sarazin 1988) given by

$$\begin{equation} \kappa _{{\rm eff}} = \frac{\kappa _{{\rm sp}}}{1 + 4.2 \lambda _e / \ell _T}. \end{equation} \tag{ 10 }$$

The rate of change of the internal energy, u, due to conduction is then

$$\begin{equation} \frac{du}{dt} = - \frac{1}{\rho } \nabla \cdot j. \end{equation} \tag{ 11 }$$

Our numerical implementation closely follows that of Parrish & Stone (2005). We use an explicit, first-order, forward time, centered space algorithm, which is the most straightforward to implement in an AMR code. Equations (8)–(11) are evaluated for each grid cell individually. Equation (11) is solved by summing the heat fluxes from all grid cell faces, and calculating the electron density and temperature on the grid cell face as the arithmetic mean of the cell and its neighbor sharing that face. The time step stability criterion for an explicit solution of the conduction equation is

$$\begin{equation} dt < 0.5 \frac{\Delta x^2}{\alpha }, \end{equation} \tag{ 12 }$$

where Δx is the grid cell size and α is the thermal diffusivity, defined as

$$\begin{equation} \alpha = \frac{\kappa }{n_e k_B}. \end{equation} \tag{ 13 }$$

Equation (12) has the potential to be considerably more constraining than the hydrodynamical Courant condition, which is proportional to T^−1/2Δx, and so the conduction timestep in conditions typical of the ICM is often much shorter than the hydrodynamic timestep. The conduction timestep is calculated on a per-level basis and is taken to be the minimum of all such values on a given level, ensuring that Equation (11) is applied in a stable manner for all grid cells. Enzo pads each AMR grid patch with three rows of ghost zones from neighboring grids. This allows the conduction routine to take three sub-cycled time steps for every hydrodynamic step, and thus allows the minimum grid timestep to be a factor of three larger than Equation (12). After three steps, temperature information from the outermost ghost zone has propagated to the edge of the active region of the grid, so performing additional conduction cycles would yield inconsistent solutions with neighboring grids. By default, Enzo rebuilds the adaptive mesh hierarchy for a given refinement level after each timestep, which is computationally expensive. However, because conduction does not explicitly change the density field (the value of which is the only quantity that determines mesh refinement), we modify this behavior such that the hierarchy is only rebuilt after an amount of time has passed equivalent to what the minimum timestep would be if conduction were not enabled.

This implementation models the case of isotropic Spitzer conduction, where heat flows unimpeded along temperature gradients. If magnetic fields are present and strong enough that the electron gyroradius is small compared to the physical scales of interest (which is likely to be true in the ICM), then heat is restricted to flow primarily along magnetic field lines. Therefore, the heat flux becomes the dot product of the temperature gradient with the magnetic field direction (i.e., j = −κbb · ∇T, where b is the unit vector pointing in the direction of the magnetic field). Diffusion perpendicular to magnetic field lines is generally considered to be negligible. We do not consider magnetic fields in this work, but instead approximate the suppression of conduction by magnetic fields by adding a suppression factor, f_sp, varying from 0 to 1, to the heat flux calculation. In reality, the strength and orientation of magnetic fields in galaxy clusters is poorly understood. Hence, the degree to which conduction is suppressed from its maximum efficiency is not known. The extreme limit of tangled magnetic fields is equivalent to isotropic conduction with f_sp = 1/3. A number of works have shown that the presence of magnetic fields alongside conduction can create a variety of instabilities that greatly influence the effective level of isotropic heat transport (e.g., Parrish & Stone 2005, 2007; Parrish et al. 2009; Ruszkowski & Oh 2010; McCourt et al. 2011, 2012). For this reason, we simulate a large range of values of f_sp, from strongly suppressed (f_sp = 0.01) to fully unsuppressed (f_sp = 1).

The presence of conduction leads to higher gas temperatures both inside and outside of the conduction-free cluster's characteristic temperature peak at r/r₂₀₀ ∼ 0.07. Within r/r₂₀₀ ∼ 0.07, gas temperatures in conductive clusters are greater than those in the non-conductive control simulation, with values of f_sp ≳ 0.1 leading to nearly isothermal cores and central temperatures ∼40% greater than those in non-conductive cluster simulations. The temperature profiles of our conductive clusters are in reasonably good agreement with those of Ruszkowski et al. (2011) for their cluster with isotropic conduction and f_sp = 1/3, although the core of their cluster is less isothermal than ours, with cooler temperatures near the center. This could be because the simulations of Ruszkowski et al. (2011) did not include a prescription for star formation and feedback in cluster galaxies, which can increase the core temperature both by consuming the cold gas and through thermal and mechanical feedback. Agreement with the temperature profiles of the simulated clusters from Dolag et al. (2004) is not nearly as good. Those clusters show a significant reduction in both the peak and central temperatures with conduction present, which we do not observe. The theoretical temperature profiles of McCourt et al. (2013) for the 10^14.5 M_☉ clusters show an elevation in the core temperatures in conductive clusters in rough agreement with our results. We also see a marginal inward movement of the location of the temperature peak as they predict, although not nearly to the same degree. The temperature peak for our clusters is also at a smaller radius to begin with.

Despite its ability to maintain approximate isothermality in the cluster cores, conduction at even maximal efficiency is unable to avert a cooling catastrophe at the very center of the halo. This cooling catastrophe and the resulting condensation and star-formation activity produce sharp peaks at r/r₂₀₀ ≲ 0.03 in all the gas density profiles in Figure 4, as well as a slight temperature increase at the same location in the more highly conductive clusters. The elevated central temperatures at the centers of conductive clusters reduce the central gas density relative to the control clusters by ∼20%–30%. Material that would have fallen into the center is displaced to larger radii, as can be seen by the elevated density at 0.04 ≲ r/r₂₀₀ ≲ 0.6 in the more conductive runs. Gas density is enhanced primarily in the range 0.04 ≲ r/r₂₀₀ ≲ 0.2, but the enhancement appears to saturate at ∼20% for f_sp = 0.33. For f_sp = 1, the density is actually lower than for f_sp = 0.33 in this range, but is then higher for 0.2 ≲ r/r₂₀₀ ≲ 0.6. It is unclear what causes the outward transport (or prevention of inflow) to stall just inside 0.2 r₂₀₀, but it seems that it can be overcome for some value of f_sp above 0.33, likely closer to 1.

Further evidence of conduction-driven inflation of the core can be seen in Figure 5, where we plot averaged profiles of M/r, including the individual contributions from dark matter, stars, and gas. Despite the lack of perfect monotonicity in the density and temperature profiles, the gaseous component of the potential decreases monotonically with increasing f_sp, indicating that the clusters are indeed responding to the presence of conduction by puffing up their cores and redistributing gas to larger radii. As is the case for most simulated galaxy clusters, the stellar component is extremely centrally concentrated, dominating the gravitational potential inside 0.04 r₂₀₀. The stellar and dark matter components of the potential are largely unaffected by conduction.

In the right panel of Figure 4, we plot normalized entropy profiles, where the normalization term, K₂₀₀ (Voit 2005; Voit et al. 2005), is given by

$$\begin{equation} K_{200} = k_BT_{200} \bar{n_{e}}^{-2/3}, \end{equation} \tag{ 14 }$$

where T₂₀₀ is

$$\begin{equation} k_BT_{200} = \frac{{G M}_{200} \mu m_{p}}{2 r_{200}}, \end{equation} \tag{ 15 }$$

and $\bar{n_{e}}$ is the average electron number density within r₂₀₀, assuming a fully ionized plasma of primordial composition. The decrease in density and increase in temperature in the cluster cores with conduction yield entropy profiles that are enhanced by 40%–70%, but still show the steep decline in the very center that is typical of simulated clusters. This enhancement decreases out to r/r₂₀₀ ∼ 0.07, where the mean entropy values are approximately equivalent for all values of f_sp. The entropy profiles for the clusters with conduction then dip below the control sample by roughly 10% out to r/r₂₀₀ ∼ 0.2. In the range 0.2 ≲ r/r₂₀₀ ≲ 0.6, the combined temperature and density enhancements appear to be in perfect balance, producing nearly identical entropy profiles across the entire cluster sample with a very small variance.

Voit et al. (2005) find that for a sample of non-radiative, non-conducting clusters simulated with Enzo, the entropy profiles in the range 0.2 ≲ r/r₂₀₀ ≲ 1 are best fit by the power law

$$\begin{equation} K(r) / K_{200} = 1.51 (r / r_{200})^{1.24}, \end{equation} \tag{ 16 }$$

which we overplot on Figure 4 with a black dashed line. This power law matches our sample well in normalization, but has a slightly steeper slope. We find that our sample is best matched by

$$\begin{equation} K(r) / K_{200} = 1.50 (r / r_{200})^{1.09}. \end{equation} \tag{ 17 }$$

For comparison, we also plot the predicted entropy profile for a pure cooling model (Voit & Bryan 2001; Voit et al. 2002) of a 5 keV cluster. This model assumes gas and dark matter density distributions that follow an Navarro–Frenk–White profile (Navarro et al. 1997) with concentration c = 6. The gas is initially in hydrostatic equilibrium and allowed to cool for a Hubble time. The pure cooling model agrees quite well with our clusters in the range 0.2 ≲ r/r₂₀₀ ≲ 1, and especially well in the range 0.2 ≲ r/r₂₀₀ ≲ 0.5. Outside of 0.5 r₂₀₀, the slope of the pure-cooling model is slightly too steep and is closer in slope to the Voit et al. (2005) fit.

The outskirts of galaxy clusters, while having lower temperatures and hence lower conductivities, are not subject to intermittent heat injection from stellar feedback, and thus offer an intriguing laboratory for studying the effects of conduction. We identify two distinct regions in the cluster outskirts where conduction appears to have influence, at ∼r₂₀₀ and at ∼3r₂₀₀. At r/r₂₀₀ ≳ 0.6, the density excess seen in the cluster interiors turns into a deficit for all values of conduction simultaneously, as can be seen in Figure 4. From this point out to well past the virial radius, there exists a perfectly monotonic trend of lower gas densities for higher values of f_sp. For the maximum value of f_sp, the average gas density at the virial radius is 10% lower than that of the clusters simulated without conduction. In fact, the density is measurably lower for all values of f_sp ⩾ 0.1 out to a few virial radii. The temperature just inside r₂₀₀ is marginally higher, while the temperature just outside r₂₀₀ is reduced.

Conduction transports heat outward in these regions because the gravitational potential there produces a declining gas temperature gradient. This heat transfer causes the entropy of the outer gas to increase. However, because gas in the cluster outskirts is not pressure-confined like the gas in the core, it is free to expand outward and decrease in both density and temperature, while its temperature gradient remains determined by the gravitational potential. This happens because the timescale for conduction in the outer regions is substantially greater than the sound-crossing time. Consequently, outward conduction of heat causes the outer gas to expand without much change in the temperature gradient. Beyond r₂₀₀, however, we see a slight steepening in the temperature profile as predicted for conductive clusters by McCourt et al. (2013), presumably because inflation of the ICM due to conduction pushes gas near the virial radius farther from the cluster center without adding much thermal energy. The fact that the entropy profiles of conductive clusters are lower than those of non-conductive clusters beyond the virial radius supports this interpretation.

Further out, at ∼3r₂₀₀, there is another systematic decrease in both density and temperature. Interestingly, as was pointed out by Skillman et al. (2008), this is the typical location for a galaxy cluster's accretion shocks, which are responsible for heating gas up to the virial temperature. This raises a critical question: How can conduction affect a large halo's accretion shocks?

To test this question in a more controlled environment than a cosmological simulation, we performed a pair of one-dimensional simulations designed to mimic the conditions of an accretion shock around a galaxy cluster. Following the shock properties characterized by Skillman et al. (2008), we initiated a Mach 100 standing shockwave with preshock conditions of n = 10⁻⁴ cm⁻³ and T = 10⁴ K, comparable to gas that is falling directly onto a galaxy cluster (i.e., through spherical accretion of the surrounding intergalactic medium, rather than being accreted via filaments). We ran simulations with f_sp = 0 and 1, with profiles shown in Figure 6. Almost immediately, a conduction front in the f_sp = 1 simulation races ahead of the original shock. The conduction front, referred to as the foreshock (Tsai et al. 2002; Fadeyev et al. 2002; Tsai et al. 2007), continues to advance with ever-decreasing speed and settles into a near-steady state within 50–100 Myr, roughly 30 kpc upstream of its initial location. The final position of the foreshock can be related to the conduction timescale (Tsai et al. 2007; Parrish et al. 2009) given by

$$\begin{equation} t_{{\rm cond}} = \frac{L^{2}}{\alpha }, \end{equation} \tag{ 18 }$$

where L is a length scale and α is the thermal diffusivity (Equation (13)). The average velocity of the conduction front over the length, L, is then inversely proportional to L. One can then calculate the value of L at which the average velocity of the conduction front is equal to the upstream velocity in the frame of the original shock, here ∼1100 km s⁻¹. Using the average of the upstream and downstream temperature and density to calculate α, we arrive at a value of L ∼ 28 kpc, which agrees very well with Figure 7. This behavior has been seen in previous works, simulating solar coronae and stellar winds (e.g., Tsai et al. 2002, 2007; Fadeyev et al. 2002). The transfer of heat upstream causes a slight reduction in the pressure at the location of the original shock front, pushing the true shock downstream. As seen in Figure 7, the motion of the true shock proceeds much more slowly, requiring roughly 1 Gyr to reach an equilibrium position about 15 kpc downstream.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The lower right panel of Figure 6 shows the Mach numbers determined using the shock-finding algorithm described in Skillman et al. (2008). Using this method, we calculate the Mach number separately using the temperature (blue lines) and density (red lines) jumps. Since conduction has led to a physical separation of the two jumps, these two methods will give different results. The true shock, defined as the position of the velocity discontinuity and hence where the kinetic energy is converted to thermal energy, has a Mach number of just under 3 as calculated from the density jump. The foreshock has a temperature jump corresponding to a Mach number of roughly 70. While the foreshock is not the true shock, measuring its equivalent Mach number is important since the primary role of an accretion shock is to heat infalling gas to the virial temperature of the halo.

These results are independent of resolution. The simulations shown in Figure 6 are for a grid 128 cells across, and we observe nearly identical behavior down to a resolution of 16 cells. For this configuration, the results are strongly dependent on f_sp. At f_sp = 0.67, the distance between the separated shocks is approximately half of that at f_sp = 1 and the Mach number of the foreshock is only reduced to just under 90. However, the Mach number of the true shock, based on its density jump, is reduced to nearly the same value as the case of f_sp = 1. It appears that if conduction is strong enough to separate the foreshock form the true shock, the resulting Mach number of the true shock, unlike the foreshock, is relatively insensitive to the value of f_sp. For f_sp ≲ 0.4, the results are indistinguishable from those without conduction. However, we find that it is possible to produce a significant shockwave alteration for lower values of f_sp simply by lowering the initial Mach number and increasing the preshock temperature (resulting in gas in a thermodynamic regime comparable to gas that is being accreted from cosmological filaments). As long as the postshock gas is able to reach temperatures in the range of 10⁷ K, where the Spitzer conductivity becomes considerable, conduction is capable of bifurcating the shock front.

Figure 8 shows the average Mach number calculated from temperature jumps as a function of gas density in our simulations of the most massive halo for all shocks in the subvolume in which refinement is allowed. Conduction reduces the Mach numbers of the strongest shocks, which occur preferentially in the lowest density gas, by roughly 10% from f_sp = 0 to f_sp = 1. We therefore conclude that the density, temperature, and entropy deficits observed beyond the virial radius in our conductive clusters are indeed due to shock alterations similar to those seen in our idealized, one-dimensional simulations of conductive shock fronts. However, as was seen in the one-dimensional shock tube simulations, identification of the true shock is not as straightforward when conduction is present. Therefore, we caution the reader that the precise Mach numbers of the accretion shocks may be somewhat different.

Zoom In Zoom Out Reset image size

However, the physics of the actual accretion shocks around real galaxy clusters is undoubtedly more complex. In particular, it is important to note that the electron mean free path in that gas is several times greater than the intershock distance in our one-dimensional simulations, and furthermore that real accretion shocks are likely to be collisionless and magnetically mediated. Nonetheless, the general qualitative point these simulations illustrate is interesting: heating of preshock gas by a hot electron precursor has the potential to alter the expected relationships between the sizes and locations of the density and temperature jumps in accretion and merger shocks. The displacement of the foreshock and true shock is small (10s of kpc) compared to the typical location of accretion shocks around galaxy clusters (r ∼ 10 Mpc), but observation of such an effect for merger shocks may be possible.

Progress in understanding the effects of conduction on accretion shocks will require modeling the gas as a fully ionized plasma (Zel'dovich 1957; Shafranov 1957), which is beyond the scope of this work. The high Mach number of an accretion shock, combined with a postshock temperature high enough for significant heat flux, create conditions similar to a supercritical radiative shock, as described by Lowrie & Rauenzahn (2007). When simulated with a two-fluid approach, the combination of preheating the preshock medium via conduction in the electrons, electron–ion coupling, and compression of the ion fluid in the postshock region can produce a small region, known as a Zel'dovich spike, where the ion temperature is slightly higher than the equilibrium postshock temperature (Lowrie & Rauenzahn 2007; Lowrie & Edwards 2008; Masser et al. 2011). Our single-fluid simulations, despite their limitations, produce shocks quite similar to the non-equilibrium results of Lowrie & Edwards (2008), where a diffusion term proportional to T^7/2 (similar to the T^5/2 of Spitzer conductivity) is used. Nevertheless, a more detailed study of the characteristics of accretion shocks employing a two-fluid MHD treatment along with physically motivated conduction and cooling rates seems warranted.