2D Magnetohydrodynamics simulations of induced plasma dynamics in the near-core region of a galaxy cluster Free

Abstract

1 INTRODUCTION

Clusters of galaxies are the largest structures in the universe. The majority of their baryonic mass is in the form of a hot plasma filling the space between galaxies, the intracluster medium (ICM). The cooling time of the ICM near the centre of the clusters due to radiative losses via X-ray emission is shorter than the Hubble time for most clusters (Sarazin 1986; Edge, Stewart & Fabian 1992; Peres et al. 1998; Sanderson, Ponman & O’Sullivan 2006). In the absence of heating, the short cooling times (as short as ∼0.1 Gyr) are expected to lead to cooling of the ICM plasma and its inflow towards the centre of the cluster at a rate that was estimated to be as large as ∼10³ M_⊙ yr⁻¹ (Cowie & Binney 1977; Fabian & Nulsen 1977; Mathews & Bregman 1978; for a review, see Fabian 1994). However, high-resolution X-ray spectroscopy from XMM–Newton and Chandra have failed to find evidence of such ‘cooling flows’ (for a review see Peterson & Fabian 2006), strongly suggesting that one or more heating mechanisms must be balancing radiative cooling in galaxy clusters.

A wide variety of heating mechanisms have been investigated, including (i) heat transport from the outer regions of the cluster to the central cooling gas by conduction (e.g. Binney & Cowie 1981; Tucker & Rosner 1983; Bertschinger & Meiksin 1986; Bregman & David 1988; Gaetz 1989; Rosner & Tucker 1989; David, Hughes & Tucker 1992; Pistinner & Shaviv 1996; Narayan & Medvedev 2001; Voigt et al. 2002; Fabian, Voigt & Morris 2002; Zakamska & Narayan 2003; Kim & Narayan 2003a; Markevitch et al. 2003; Dolag et al. 2004; Pope et al. 2006; Xiang et al. 2007; Balbus & Reynolds 2008; Bogdanović et al. 2009; Parrish, Quataert & Sharma 2009; Parrish, Quataert & Sharma 2010); (ii) energy injection from active galactic nuclei (AGNs) via jets, bubbles, radiation and sound waves (e.g. Binney & Tabor 1995; Ciotti & Ostriker 2001; Brüggen & Kaiser 2002; Omma et al. 2004; Ruszkowski, Brüggen & Begelman 2004a,b; Voit & Donahue 2005; Fabian et al. 2005; Nusser, Silk & Babul 2006; Mathews, Faltenbacher & Brighenti 2006; Chandran & Rasera 2007; McCarthy et al. 2008; Brüggen & Scannapieco 2009); (iii) heating due to galaxy motions via dynamical friction (e.g. Schipper 1974; Miller 1986; Just et al. 1990; Fabian 2003; El-Zant, Kim & Kamionkowski 2004; Kim, El-Zant & Kamionkowski 2005; Dekel & Birnboim 2008); (iv) turbulent mixing (e.g. Loewenstein & Fabian 1990; Deiss & Just 1996; Ricker & Sarazin 2001; Cho et al. 2003; Kim & Narayan 2003b; Dennis & Chandran 2005; ZuHone, Markevitch & Johnson 2010); (v) magnetic field reconnection (Soker & Sarazin 1990); vi) viscous heating regulated by pressure anisotropy (e.g. Kunz et al. 2010) or (vii) a combination of these (e.g. Brighenti & Mathews 2002; Ruszkowski & Begelman 2002; Brighenti & Mathews 2003; Voigt & Fabian 2004; Fujita & Suzuki 2005; Pope et al. 2006; Conroy & Ostriker 2008; Guo, Oh & Ruszkowski 2008).

Due to the extremely high Hall parameter of ICM plasmas, the impact of the magnetic field on the thermal balance of the ICM has been of particular interest; extensive multidimensional magnetohydrodynamics (MHD) simulations have been employed in the last several years to assess anisotropic effects. Parrish & Quataert (2008) performed 2D and 3D simulations with the MHD code athena (Gardiner & Stone 2008; Stone et al. 2008) focusing on the evolution and saturation of the heat-flux-driven buoyancy instability (HBI; Quataert 2008) for various magnetic field strengths. Bogdanović et al. (2009) also performed 3D simulations with athena. They carried out pre-simulations, whereby starting with an isothermal ICM they allowed radiative cooling without thermal conduction to obtain a cooling core state, and then used it as the ‘initial condition’ for their self-consistent MHD simulations. For the initial magnetic field, the authors’ numerical experiments included a case of rotational field lines. They concluded that thermal conduction alone cannot be the solution to the cooling flow problem. Moreover, they hypothesized that the nature of the AGN feedback is one of ‘stirring’ rather than heating, that is, AGN activity periodically disrupting the azimuthal field structure allowing thermal conduction to sporadically heat the core. The level of effectiveness of ‘stirring’ however remains an open question. For example, cooling-core clusters like the A2199 contain prominent metal gradients that could prevent violent processes from stirring the ICM. In a similar study, Parrish et al. (2009) performed a wide range of numerical studies of the A2199 cluster with tangled and radial initial magnetic field topologies and strengths of 1 nG to 1μG, extending the parameter study of Bogdanović et al. (2009). They demonstrated that the effective reduction of the anisotropic heat flux was less than 10 per cent of the Spitzer value due to a re-orientation of the magnetic field lines relative to the radial direction induced by the HBI. This reduction was found to occur within several times the HBI growth time (∼0.125thinsp;Gyr for A2199), leading to a cooling catastrophe in less than ∼3 Gyr. The authors found similar trends at different initial magnetic field strengths and topologies. Starting with lower initial core densities and isothermal initial temperature the cooling catastrophe was significantly delayed to ∼7 Gyr. The cooling catastrophe was also delayed when a central heating source was assumed via an idealized model of the heating luminosity within a radius of 20 kpc. The authors concluded (at least for A2199) that in the absence of AGN feedback and/or strong magnetic fields near the galaxy cluster core, thermal conduction alone cannot be responsible for balancing the radiative losses in the ICM, in agreement with Bogdanović et al. (2009). Similar results about the reduction and evolution of the effective heat flux were obtained by Ruszkowski & Oh (2010) using the 3D MHD code flash.1

There is a large body of numerical work on various aspects of AGN–ICM interactions that has spanned purely hydrodynamic models (e.g. Churazov et al. 2001; Brüggen & Kaiser 2001, 2002; Reynolds, Heinz & Begelman 2002; Omma et al. 2004; Zanni et al. 2005; Vernaleo & Reynolds 2006), as well as MHD models (e.g. Ruszkowski et al. 2004a; Robinson et al. 2004; Jones & De Young 2005; De Young 2010). Yet, a conclusive picture regarding the cooling flow problem based on numerical simulations appears to remain elusive. For example, Vernaleo & Reynolds (2006) simulated direct injection of a supersonic jet into the ICM atmosphere. However, this model failed to maintain long-term thermal balance in the cluster, leading to hypotheses of heating by the dissipation of sound waves or the decay of global ICM modes (e.g. Omma et al. 2004). Also, De Young (2010) argue based on comparisons with experimental data that AGN outflows are strongly decelerated, generating turbulent sonic or subsonic flows due to their interaction with the surrounding medium, which support such hypotheses.

In this paper we present results from a series of 2D numerical experiments aimed at assessing the driving physics and related sensitivities in a galaxy cluster. We have chosen to simulate A2199, a cooling-core cluster that has been simulated by other authors (e.g. Zakamska & Narayan 2003; Parrish et al. 2009; Ruszkowski & Oh 2010), and for which observed data for the electron number density and temperature exist at a relatively close distance (R∼ 1 kpc) from the cluster centre. The density at the core is an important parameter in MHD parameter studies that aim to assess mechanisms associated with the sustainment (or collapse) of the thermal balance in a cluster, since the radiative cooling rate depends strongly on it (∝n²_e). Because A2199 is a cluster that is often modelled by MHD simulations, we perform here a series of idealized numerical experiments in 2D axisymmetric geometry to assess the sensitivity of the collapse rate on the near-core density, using two different spatial profiles that bounded the observational data (Johnstone et al. 2002).

In 3D MHD simulations, it is common to prescribe initial profiles that are in hydrostatic equilibrium. Often, the observational data guide solutions to the hydrostatic equation and then these solutions are prescribed as the ‘initial conditions’ of the ICM. MHD simulations then proceed to seek of mechanisms that maintain thermal balance in the cluster. In this study we perform a series of numerical experiments that do not begin with solutions in hydrostatic equilibrium, to assess dynamical effects in the energy balance of the system. Our simulations include a range of initial magnetic field profiles (similar to those used by Parrish et al. 2009 and Bogdanović et al. 2009) and quantify the self-consistent evolution of the ICM in 2D planar geometry.

The paper is organized as follows. Section 2 describes the computational methodology and governing laws used to perform the numerical simulations. Initial and boundary conditions for the cluster are also given in this section. Section 3 presents our simulation results in 2D axisymmetric (r, z) and in 2D planar (x, y) physical domains. We summarize the results and provide our conclusions in Sections 4 and 5, respectively.

2 METHOD

The numerical simulations have been performed with the Multiblock Arbitrary Coordinate Hydromagnetic (mach) code (Peterkin 1998). machx (where x = 2 or 3) are time-dependent, resistive, compressible MHD solvers of the Arbitrary Lagrangian Eulerian variety that were developed by the Center for Plasma Theory and Computation at the Air Force Research Laboratory and NumerEx Inc., to support applications with non-ideal MHD physics in complex geometries. mach2 and mach3 are the two-and-one-half-dimensional [2(1/2)D] and 3D versions of the solvers, respectively. Over the last three decades, they have been upgraded significantly through the contributions of many scientists in the United States. Although the majority of the work with mach today is in the area of pulsed power, their application base has expanded to multiple areas in engineering. However, to the best of our knowledge, they have not yet been used to study problems in astrophysics. In order to introduce these codes to the astrophysical community, we provide a more detailed description of the codes and their applications in Appendix A.

mach2 is the version used in the present study. A feature of mach2 that has been particularly useful in our studies is its modularity with the terms in the MHD equations, which allowed us to perform instructive comparisons between the various simulation physics in the evolution of the ICM. The numerical experiments have been carried out using a reduced set of the full MHD conservation laws of mach (see Appendix A). Simulations in only two dimensions allow a wide range of highly resolved cases at reduced computational time, but at the expense of excluding processes in the evolution of the ICM that are inherently 3D. Nevertheless, the range of physics cases examined has provided us with a useful set of results, especially in regards to the ‘initial state’ of the ICM and the potential importance of low-speed hydrodynamics near the centre of the cluster. These 2D results now also form the basis of our future simulations with mach3.

2.1 Governing equations and numerical approach

The ICM simulations have been performed using the following reduced set of the full mach conservation laws, representing ideal MHD physics in the presence of anisotropic thermal conduction and radiative cooling:

1

2

3

4

In equations (1)–(4), ρ is the mass density of the plasma fluid, v is the fluid velocity field, B is the magnetic induction field, ε is the internal specific energy, q is the conductive heat flux, p is the thermal pressure and g is the gravitational acceleration.

Analytic models or semi-empirical tabular equations of state (EOS) and transport coefficients may be used in mach by accessing the SESAME library2 (see Appenix A for more information about the library) during each computational cycle. For the cluster simulations, we have used an ideal-gas EOS with a specific heat ratio γ = 5/3,

5

Also, an analytic model of the anisotropic thermal conductivity coefficients has been used, which incorporates the full Braginskii coefficients (Braginskii 1965) for the parallel (κ_e//) and perpendicular (κ_e⊥) conductivities

6

where is the thermal conductivity tensor, is the unit vector in the direction of the magnetic field and is the delta tensor; is the magnetic field unit dyad. The Braginskii coefficients for are provided in Appendix A. The conduction heat flux is given by

7

In the relevant range of temperatures (∼1–10 keV) the plasma in the cluster cools primarily by thermal bremsstrahlung at a rate given by (in units of J m⁻³ s⁻¹) (e.g. Rybicki & Lightman 1979) with denoting the velocity-averaged Gaunt factor and n_e is the electron number density. mach2 uses a mean Planck opacity χ_p (with dimensions of area cross-section per unit mass, so a photon mean-free-path is 1/χ_p). The opacity values are taken from the SESAME library. In the optically thin limit, radiation energy is removed from a mach2 computational element at a volumetric cooling rate of

8

where a is Stefan’s constant and c_ℓ is the speed of light in vacuum. A comparison of the cooling rates Φ_Tb and Φ_eR is provided in Appendix A.

The gravitational acceleration is obtained from Poisson’s equation

9

where ρ_DM is the mass density of the dark matter. Newton’s gravitational constant is denoted by G. The present study accounts only for the dark matter density, which is prescribed based on a Navarro–Frenk–White (NFW) profile with a soft core (Navarro, Frenk & White 1997). The density profile may be expressed in non-dimensional form as follows:

10

where ⁠, and with R denoting the radius measured from the cluster origin (r=z = 0 or x=y = 0). The normalization mass is taken to be M₀ = 3.8 × 10¹⁴ M_⊙ from Zakamska & Narayan (2003) and it has been kept fixed throughout all simulations presented herein. The standard scale radius for A2199, R_s = 390 kpc, is taken also from Zakamska & Narayan (2003) who inferred its value from observations of the outer plasma temperature. The value of the core radius R_c is uncertain. For the 10 clusters studied by Zakamska & Narayn, the authors’‘best-fitting’ solutions correspond to values of 0 and 1/20. For A2199 they chose R_c = 0, whereas Ruszkowski & Oh (2010) chose R_c = 20 kpc. We found that our main conclusions were not affected by this value and present our results for R_c = 20 kpc. The solution of equation (9) for the gravitational acceleration is spherically symmetric, ⁠, where is the position unit vector defined relative to the cluster origin. In non-dimensional form, the analytical solution to equation (9) is given by

11

where ⁠. Based on estimated values of the NFW characteristic radii and normalized mass for A2199, we find ρ_DM≈ 1.3 × 10⁻²⁴ g cm⁻³ at R = 0.95 kpc. At the same location the observed electron number density is 0.115 cm⁻³ and is the highest value we have used in our initial conditions. This gives ρ≈ 0.2 × 10⁻²⁴ g cm⁻³ for the plasma assuming that the molecular weight per electron is μ_e = 1.18. The latter assumed that the hydrogen and helium fractions are X = 0.7 and Y = 0.28, respectively (Zakamska & Narayan 2003). Thus, our assumption in equation (9) about the dominance of the dark matter density in the determination of the gravitational acceleration appears acceptable.

Equations (1)–(4) are solved in a fractional time-split manner. The Lagrangian time advance of the mass density, velocity and magnetic vector fields are computed implicitly on a computational domain that, unless otherwise noted, is discretized by a uniform mesh that consists of square computational elements (or ‘cells’). Simulations in both 2D axisymmetric (r, z) and 2D planar (x, y) coordinate systems have been performed but no mesh movement and/or adaptation has been invoked. Spatial discretizations are obtained using finite volume differencing. The radiation cooling algorithm advances the internal energy explicitly, whereas thermal conduction is performed implicitly. Because we have assumed that the gravitational acceleration is dependent only on the dark matter density, g is prescribed using equation (11) and is held fixed throughout each simulation.

2.2 Initial and boundary conditions

2.2.1 Initial conditions for the plasma density and temperature

Several combinations of initial conditions for the scalars and vectors have been used throughout this study, all associated with the A2199 galaxy cluster. Specifically, for the number density and temperature of the plasma, we employ two different profiles to illustrate the sensitivity of the ICM simulations on the assumed initial conditions. The first profile is a fit of the Zakamska & Narayan ‘best’ solution to an idealized set of four, time-independent, ordinary differential equations that account for thermal conduction and radiation cooling only, in hydrostatic equilibrium. The authors sought the solution that best fits the Chandra observations (Johnstone et al. 2002) by adjusting the thermal conduction heat flux factor f, defined here as

12

where ⁠. The Spitzer–Braginskii thermal conductivity is κ_e//=κ_sp = 3.16k_B²n_eT_eτ_e/m_e, where τ_e is the (classical) e–i collision time. The authors showed that for five out of the 10 clusters studied a pure conduction model balancing radiative cooling, i.e.

13

was plausible if the factor f was reduced from unity by more than 50 per cent. Because this best-fitting solution yields an electron number density near the core which is about 35 per cent lower than the observed value, we shall label it Low Core Density or ‘LCD’. The LCD profiles are plotted in Fig. 1 against Chandra data (Johnstone et al. 2002, also in Zakamska & Narayan 2003).

Figure 1

Chandra data (Johnstone et al. 2002, also in Zakamska & Narayan 2003) and initial profiles used in the A2199 galaxy cluster simulations.

Open in new tabDownload slide

Some MHD simulations in the literature have employed initial density profiles that underestimate the observed electron number density by factors of 5–6 near the cluster core [e.g. see Parrish et al. (2009) and Ruszkowski & Oh (2010)]. To assess the importance of the core conditions on the evolution of the ICM, we constructed a second set of initial conditions in a manner similar to that followed by Ruszkowski & Oh (2010), who assumed hydrostatic equilibrium and a polytropic EOS with a prescribed profile for the entropy S (e.g. Cavagnolo et al. 2009). We write the entropy as a function of radius in non-dimensional form:

14

where the barred quantities are defined as ⁠, ⁠, ⁠, β≡ (T₂₀₀/T_c)(ρ_c/ρ₂₀₀)^{γ− 1} and c=R₂₀₀/R_s is the concentration parameter. Here, variables carrying subscript ‘c’ represent values at the cluster centre and those with a subscript ‘200’ denote values at a radius where the mean density of the dark matter is 200 times the critical density of the universe at the redshift of the cluster ρ_crit(z).

The hydrostatic equation is given by

15

Assuming spherical symmetry, equation (15) may be written in terms of the non-dimensional density and entropy as follows:

16

and, given ⁠, yields the density as a function of radius. The gravitational acceleration is taken from equation (11). The constant on the right of equation (15) is C≡ 2GM₀μm_u/R_sk_BT_c (μ denotes the mean molecular weight, m_u is the atomic mass unit and k_B is Boltzmann’s constant). Equation (15) has been solved numerically to produce a second set of cluster profiles by assuming the following parameters: T_c = 1.6 keV, T₂₀₀ = 4 keV, ρ_c = 23.5 × 10⁻²⁶ g cm⁻³, ρ₂₀₀ = 0.054 × 10⁻²⁶ g cm⁻³, c = 4.378, γ = 5/3 and μ = 0.62. The profiles are plotted in Fig. 1 and shall be denoted by High Core Density or ‘HCD’ hereinafter. Although neither the LCD nor the HCD profiles represent accurately the density observations, they were implemented in this manner deliberately as limiting solutions that bounded the full range of the data, within the constraints of the hydrostatic equilibrium condition, while illustrating the sensitivity of the simulation results on the density.

2.2.2 Initial conditions for the magnetic field

An extensive body of work exists on the magnetic field structure in galaxy clusters covering both cooling-core and non-cooling-core clusters. Faraday rotation measures towards extended radio sources in and behind galaxy clusters have demonstrated the turbulent nature of the magnetic field there (e.g. see Murgia et al. 2004; Vogt & Enßlin 2005; Laing et al. 2006; Guidetti et al. 2008; Bonafede et al. 2010; Govoni et al. 2010; Guidetti et al. 2010; Vacca et al. 2010). In general, it has been found that a power-law distribution of their fluctuation spectrum and some radially declining profile of their average strength represent well these magnetic field structures. Nearly all hydrodynamical simulations within a cosmological context predict around 10 per cent of turbulence within the ICM (e.g. see Bryan & Norman 1998; Inogamov & Sunyaev 2003; Gardini et al. 2004; Vazza et al. 2006; Iapichino et al. 2008; Paul et al. 2010), and are in agreement with observations (e.g. see Churazov et al. 2004; Schuecker et al. 2004). Consequently, since the magnetic field is frozen within such turbulent ICM, cosmological MHD simulations predict that the magnetic field has a turbulent spectrum (specific examples can be found in Dolag, Bartelmann & Lesch 2002; Brüggen & Hoeft 2005) with no exception when studying a sample of clusters or during the time evolution of a single cluster. The reason for this is that clusters forming in a cosmological context contain accreting matter, part of which forms by small-scale-structure mergers. Therefore, clusters may be populated by thousands of substructures in the form of galaxies or former groups of galaxies that move through the cluster with speeds that are typically on the order of a 1000 km s⁻¹, permanently perturbing the underlying potential. It has been shown that such gravitational perturbations alone can produce significant motions within the core of clusters (e.g. see Ascasibar & Markevitch 2006; Roediger et al. 2010; ZuHone, Markevitch & Johnson 2010). Some very small fraction of these substructures contains also gaseous atmospheres or discs capable of steering the ICM directly.

To assess the impact of initial magnetic field topologies on the evolution of the effective thermal conduction in the ICM, we employed here the following three initial magnetic field configurations: (1) a turbulent field (using random number deviates), (2) a purely rotational field and (3) a purely radial field. The intent of the last two highly idealized configurations was to assess the sensitivity of the anisotropic heat flux reductions, and have served as the two limiting cases of our parameter studies on magnetic fields. In all three cases, the out-of-plane component of the magnetic field, B_θ in the 2D axisymmetric cases and B_z in the planar cases, was set to zero.

The turbulent field was constructed based on the approach followed by Ruszkowski et al. (2007) (see also Ruszkowski & Oh 2010). First, we set up a random field in k-space, using independent normal random deviates N_a(0, 1) (with 0 and 1 being the values of the mean and standard deviation of the random distribution, respectively) for the real and imaginary components:

17

where the wavenumber magnitude is | k| = (k_x²+k_y²)^1/2 and subscript ‘a’ is a number used to identify different random deviates. Ruszkowski et al. (2007) proposed the following dependence of the magnetic field amplitude on the wavenumber | k|:

18

where k₀ and k₁ control the exponential cut-off terms in the magnetic energy spectra. The former controls the spectra at the large wave numbers and is given by k₀ = 2 π/λ₀ with λ₀∝ 43 h⁻¹ and h denoting the normalized Hubble parameter (0.5–0.8). The second, k₁, controls the cut-off of the spectra at the small wavelengths; the case we present in this paper uses k₀ = 0.095 and k₁ = 0.05. Once the complex vector field in k-space has been generated, we perform vector divergence cleaning. Using conventional Einstein tensor summation notation where δ_ij is the Kronecker delta, this operation may be expressed by (e.g. see Balsara 1998)

19

Finally, we perform (complex) inverse Fourier transformation in 2D to obtain the r-space components [B_x(r), B_y(r)].

Because the precise distribution of magnetic fields in galaxy clusters is uncertain, there is no strong evidence in favour of the assumed spectrum in equation (18). Therefore, we extend the range of our numerical experiments to include the following two (extreme) cases: a rotational profile and a radial profile as given by equations (20) and (21), respectively, where R₀ is a scaling radius.

20

21

These two additional cases have been used exclusively in the 2D planar simulations and the results are presented in Section 3.2. For both cases the out-of-plane component B_z was set to zero. The initial magnetic field profiles for all three cases are illustrated in Fig. 2 for a peak magnetic field value of 1 μG. Hereinafter, we shall denote the maximum strength of the magnetic field in our simulations by B.

Figure 2

Initial magnetic field profiles used in the 2D planar simulations. The x–y plots show contours of magnetic pressure (in Pa) overlaid by magnetic field unit vectors. The magnitudes of the magnetic field for the three cases used in the simulations – rotational, radial and turbulent – are compared in the lower left plot along x at y = 0. In all planar (x, y) simulations, the computational mesh was uniform consisting of 128² square cells.

Open in new tabDownload slide

2.2.3 Boundary conditions

Boundary conditions in mach2 are defined by employing ‘ghost’ cells and vertices. The simulations in planar geometry were performed using a computational domain that spanned −L≤x≤L and −L≤y≤L with L=500 kpc. At all outer boundaries in this geometry, we have imposed a Dirichlet condition for the temperature and the specified value is set equal to the initial value as given in Fig. 1. Radiative cooling is allowed at these boundaries. For each time-step the velocity at the outer boundaries is set to zero (‘no slip’), the density is set to the initial value (1) and the magnetic field is copied from the real cell adjacent to the boundary to the ghost cell.

For the axisymmetric (r, z) simulations, the computational region spans 0 ≤r≤L and 0 ≤z≤L. At the boundaries (L, z) and (r, L) the conditions are the same as those at the outer boundaries of the planar geometry. At (0, z) and (r, 0) we impose thermally insulating conditions. Only one case in this computational domain employs a magnetic field and it is kept frozen everywhere throughout the simulation.

3 NUMERICAL EXPERIMENTS

A series of numerical experiments were performed in 2D axisymmetric (r, z) and 2D planar (x, y) geometry with the two initial plasma profiles described in Section 2.2.1, LCD and HCD, and three initial magnetic field profiles described in Section 2.2.2 for a wide range of combinations and simulation physics. Hereinafter, it will be implied that when we present simulation results without a magnetic field, thermal conduction has been computed isotropically, whereas for simulations that include a magnetic field thermal conduction has taken into account the anisotropic effects (see equations 6 and 7).

Table 1 lists the cases carried out in 2D axisymmetric geometry. A goal of the simulations in this geometry was to assess the sensitivity of the global heat flux factor f on the plasma density near the A2199 core under idealized simulation physics. Specifically, no hydrodynamics were permitted in any of these simulations and the magnetic field, when used, was held frozen. Table 2 outlines our cases in 2D planar geometry. In all of these cases, no control of the thermal conductivity was exercised and the anisotropic heat flux was computed self-consistently. The focus of these simulations was the early dynamics in the near-core region of the cluster and the results are presented in Section 3.2.

3.1 2D simulations in axisymmetric (r, z) geometry

The main intent of the simulations in 2D axisymmetric geometry was to conduct a preliminary assessment of the sensitivity of the ICM solution to the prescribed initial conditions before invoking full MHD simulations. Due to the dependence of the cooling rate on the square of the electron density, ⁠, we expect the solution will depend strongly on the near-core gas density (e.g. see Conroy & Ostriker 2008). In fact, if one assumes that the dominant competing mechanisms establishing the thermal balance in the ICM are thermal conduction and radiation only, and all other mechanisms are captured by a single effective reduction of the thermal conductivity (represented by the quantity f), we can estimate of this sensitivity from the steady-state energy equation (13), by computing directly the heat flux factor f of the cluster. The value of the ICM density in the cluster core (R≈ 0.95 kpc), based on the Chandra observations, is n_e≈ 0.115 cm⁻³, whereas the LCD model at the same location gives n_e≈ 0.075 cm⁻³. Since the temperature predicted by the LCD and HCD models is approximately the same at this location, we use the LCD model to estimate ∇· (κ_sp∇T_e). By approximating the core region as a small finite cylindrical volume with dimensions Δr=Δz≈ L/128 = 3.906 kpc, we find

• f≈ 0.8 when n_e≈ 0.115 cm⁻³

• f≈ 0.5 when n_e≈ 0.075 cm⁻³,

which illustrates to first order the sensitivity of idealized solutions on the near-core plasma density.

The significance of this sensitivity may be better appreciated by comparing the following relevant time-scales. First, the characteristic time τ_q associated with the isothermalization of a region with a temperature gradient ∼ΔT_e/ΔR, may be estimated by setting

22

yielding

23

(with ℓ∼ΔR). Next, the characteristic time τ_Φ in which the plasma radiates away all of its internal energy is obtained in a similar manner yielding

24

If the mechanism(s) in the ICM produce an effective reduction of the thermal conductivity f≈ 0.5, then the conduction time would be more than double the radiation time, assuming the plasma density near the core is as high as the observed value. In other words, without any additional heating source(s), these estimates suggest that thermal conduction from the outer regions of the cluster cannot avert runaway cooling near the core.

Indeed, 2D axisymmetric simulations with mach2 that incorporated only thermal conduction and radiation (MHD was turned off in this series of calculations) show that thermal balance is achieved at a relatively low value of f (between 0.5 and 0.6) for the LCD profile, whereas the HCD profile leads to the collapse of the core for the same values of f. Fig. 3 (top) shows the evolution of the computed temperature at R=20 kpc for different values of f, when LCD initial profiles were used. Mesh sensitivity calculations confirmed that a 64² uniform mesh provided sufficient resolution. Fig. 3 (bottom) compares the LCD solutions to those obtained using the HCD profiles for f = 0.5 and 0.6, showing collapse of the core in less than 1 Gyr. Thus, determining a global value for f that would sustain thermal balance in the cluster may be in general instructive; quantitatively, however, it is of less significance due to the sensitivity of f on the central density.

Figure 3

Results from 2D axisymmetric simulations with thermal conduction and radiation only. The solutions are plotted for different values of the heat flux factor f, at R = 20 kpc. Top: results using the LCD initial conditions for 64² and 128² uniform (U) computational meshes. The 64² non-uniform (NU) mesh used quadratic spacing with highest resolution of 0.25 kpc at the cluster core. Bottom: comparisons of solutions for the LCD and HCD initial conditions at f = 0.5 and 0.6 (it is noted that differences in the two HCD solutions are negligible).

Open in new tabDownload slide

Due to the extremely high Hall parameters in the ICM (Ω_e > 10⁹ in A2199), the potential impact of turbulent magnetic fields in the transport of heat from the outer (hot) regions of the cluster has been widely recognized. Although the evolution and impact of MHD instabilities, specifically of the HBI (e.g. see Quataert 2008; Parrish & Quataert 2008) and of the magneto-thermal instability (MTI) (e.g. see Balbus 2000; Parrish, Stone & Lemaster 2008), remain an ongoing topic of study, we performed here an idealized numerical experiment with a turbulent initial magnetic field (see Section 2.2.2) that was kept frozen throughout the simulation. The intent was to emulate a scenario in which continuous galaxy motions sustain a nearly isotropic turbulent field. The simulation case was motivated in part by the work of Ruszkowski & Oh (2010), who suggest it is possible that turbulent velocity fields in the cluster may drive (or stir) turbulence in the magnetic field in a manner that allows for sufficiently rapid transport of heat from the outer regions of the ICM to sustain thermal balance.

Fig. 4 (top) shows contours of the initial temperature overlaid by unit vectors of the initial magnetic field as implemented in our 2D axisymmetric simulation. For this case the peak magnitude of the magnetic field did not exceed 1 μG. The evolution of the computed temperature at R=R_c is shown in Fig. 4 (bottom). Also shown for comparison is the case with no magnetic field. For all cases in this series of simulations the thermal conductivity was computed self-consistently, accounting fully for the anisotropy in the heat flux, but MHD were excluded. We found that thermal conduction alone is unable to avert a cooling catastrophe, even when the magnetic field is completely tangled. The collapse of the core for the 1-μG case was found to be insensitive to the value of peak magnetic field, which is expected due to the high values of Ω_e. In this arrangement a change in the collapse rate was obtained only when the peak magnetic field magnitude was lowered by several orders of magnitude (where Ω_e→ 1).

Figure 4

Results from 2D axisymmetric simulations with anisotropic thermal conduction, radiation and a turbulent magnetic field. The magnetic field was frozen throughout the simulation. Top: temperature contours overlaid by magnetic field unit vectors at t = 0. Bottom: evolution of the average temperature at R = 20 kpc for different values of the peak magnetic field amplitude allowed. It is noted that differences between the solutions for the 0.1 and 1 μG are negligible.

Open in new tabDownload slide

By comparison, Ruszkowski & Oh (2010) found in their 3D simulations of A2199 that a cooling catastrophe could be averted. The authors applied a driving force that resulted not only in keeping the magnetic field tangled but also in turbulent flows that introduced additional physical effects such as mixing, which we do not account for in our 2D axisymmetric simulation. We also recognize that the real topology of the magnetic field is unknown and may differ significantly from the idealized model described in Section 2.2.2. Nevertheless, in view of the sensitivity of thermal balance calculations on the plasma density it appears that hypotheses about ‘stirring’ or other dynamical activity in the ICM, proposed as mechanisms that could sustain thermal balance in the cluster, would only be strengthened if such mechanisms sustained balance in the cluster for a wide range of possible core densities. Ruszkowski & Oh (2010) employed an initial temperature profile that was in close agreement with the Chandra data and a plasma density profile with values almost six times lower than the observed data near the cluster core.

3.2 2D simulations in planar (x, y) geometry and in hydrostatic non-equilibrium

In this section we perform an idealized study in two dimensions near the core of A2199, focusing on possible mechanisms that may overcome HBI-driven collapse (Parrish et al. 2009; Ruszkowski & Oh 2010). In the 3D MHD simulation cases of Parrish et al. (2009) and Ruszkowski & Oh (2010), some initial ICM states did not achieve thermal balance but all cases imposed hydrostatic equilibrium at t = 0. Here, we deliberately violate this condition only as a means to generate and study near-core flow dynamics (<20 kpc) that may possibly be associated with AGN activity.

We deviate from hydrostatic equilibrium by prescribing the spherically symmetric initial conditions LCD and HCD (Fig. 1) in a 2D planar geometry. We note that although both these profiles correspond to hydrostatic equilibrium in three dimensions, they do not preserve such equilibria in two dimensions, thereby producing subsonic wave motion in the plasma. Here we focus only on the early dynamics of the MHD flow field, t < 2.5 Gyr, near the cluster centre (R < 20 kpc), since several simulations of A2199 have predicted a cooling catastrophe due to the HBI within this time for a wide range of magnetic field strengths and profiles. Assuming a weak magnetic field initially in the direction of the gravitational force, the HBI has a characteristic growth rate (Quataert 2008)

25

≈(0.125 Gyr)⁻¹. Thus the instability is allowed approximately 20 growth times in ∼2.5 Gyr.

We find in the 2D planar simulations with a turbulent initial magnetic field (see Fig. 2 bottom right) that the average temperature at R = 20 kpc is sustained and a cooling catastrophe does not occur, as shown in Fig. 5. For comparison, the solutions with the same simulation physics are also shown in Fig. 5 for the following cases: (1) radial initial magnetic field, (2) rotational initial magnetic field, (3) turbulent initial magnetic field profile but held frozen throughout the simulation (a case with no MHD but one that accounted for anisotropic thermal conduction and radiative cooling), (4) isotropic thermal conduction and radiation without MHD and (4) isotropic thermal conduction, radiation and hydrodynamics but no magnetic field. Similar to the cases of isotropic conduction no catastrophic cooling is found in the simulations with large-scale flows (induced by starting out of hydrostatic equilibrium), independent of the initial configuration of the magnetic field and despite a significant reduction of the heat flux factor. Fig. 6 (left) plots the evolution of the cosine of the angle between the magnetic field and radial position vectors, and the heat flux factor 〈f〉 in the case of the turbulent initial magnetic field and the LCD profile. Unless otherwise noted all quantities have been averaged over the area of a ring of radius 20 kpc and maximum thickness of Δ, which denotes the side length of a single computational cell (recall our mesh consists of square cells). Also plotted in Fig. 6 (left) are the absolute radial velocity and Mach number ⁠, where c_s is the adiabatic acoustic speed.

Figure 5

Comparisons of the average temperature at R = 20 kpc, computed in planar geometry for cases with the same (LCD) profiles but different simulation physics (see Table 2 for simulation physics nomenclature). For all cases with a magnetic field the maximum strength was 1 μG.

Open in new tabDownload slide

Figure 6

Results from 2D planar simulations with CRM simulation physics starting with the LCD initial condition for the density and temperature and a turbulent initial magnetic field with maximum magnitude of 1 μG. Left: evolution of the cosine of the angle between the magnetic field and radial position vectors, ⁠, heat flux factor 〈f〉, absolute radial velocity and Mach number ⁠, averaged over the area of a disc of radius R = 20 kpc (±Δ), Right: temperature contours (in eV) overlaid by magnetic field unit vectors at t = 1.27 Gyr showing the orientation of the field to be largely rotational near the cluster core.

Open in new tabDownload slide

The evolution of the velocity field as shown in Fig. 6 (left) defines a characteristic interval of time, 0.2 < t < 1.3 Gyr. The interval is bounded by two characteristic times. The shortest time is associated with the hydrostatic imbalance that is imposed initially in the system, producing a subsonic expansion with a characteristic speed of ∼100 km s⁻¹ (〈M〉∼ 0.1). The time it takes for the expansion wave to reach a radius of 20 kpc is ∼20 kpc/100 km s⁻¹ = 0.2 Gyr. As shown in Fig. 6 (left), upon passage of this first wave the rarefaction flow at R=20 kpc moves at an approximately constant speed, ∼10 km s⁻¹ (〈M〉∼ 0.01) for about another 1.1 Gyr. Such ‘winds’ are, in fact, observed in luminous galaxies (and are commonly labelled ‘broad-line regions', e.g. see Ostriker et al. 2010). In the ensuing discussions, we focus in this ‘quiet’ rarefaction period 0.2 < t < 1.3 Gyr as an emulator of low-speed flow activity (possibly associated with an AGN). Fig. 7 plots 2D temperature contours and unit vectors of the magnetic field at t = 0, 0.16 and 0.32 Gyr showing the evolution of the expansion.

Figure 7

Evolution of the temperature and magnetic (unit vector) field for the case of the LCD initial profile and a 1-μG turbulent initial magnetic field. The calculations employed CRM simulation physics in planar geometry.

Open in new tabDownload slide

The two-way transit time for sound waves generated near the core is ∼2L/c_s = 1.3 Gyr. At about this time the heat flux factor and cosine of the magnetic field angle at R = 20 kpc have reached their minimum values of 0.17 and 0.35 (angle = 70°), respectively (note that the angle in fact decreases within t≈ 0.2 Gyr to about ∼45° from the initial value of ∼60° due to the expanding flow). Fig. 6 (right) shows the re-orientation of the magnetic field lines due to the HBI within ∼100 kpc at 1.27 Gyr, which is in qualitative agreement with the 3D simulation results of Parrish et al. (2009) and Ruszkowski & Oh (2010). Also noted is the qualitative similarity of our solution for the heat flux factor in Fig. 6 (left) with that obtained by Parrish et al. (2009) and Ruszkowski & Oh (2010) in three dimensions before saturation (to the value of 0.07). At about 1.3 Gyr, the velocity at R=20 kpc begins to increase as the first expansion wave is reflected off the outer boundaries of our physical domain and reaches the near-core region. Beyond this time the magnetic field orientation with respect to the radial position vector, the heat flux factor and the temperature are driven by the hydrodynamics of the flow.

A closer look at the various terms in the energy equation shows the diminishing importance of thermal conduction near the core in the early energy balance of the system, with the convective and radiative contributions almost fully driving the system by ∼1.4 Gyr. Fig. 8 compares the power densities for radiative cooling Φ_eR, thermal conduction ∇·q and the convective terms −(ρv·∇ε+p∇·v), at three different times along the y = 0 axis. In all cases, we found that p∇·v is the dominant term over ρv·∇ε. Also, because the initial conditions imposed zero velocity, the plot at t = 3.2 × 10⁻⁵ Gyr simply reflects that the evolution of the internal energy of the system is initially driven by thermal conduction and radiation, but by ∼1.4 Gyr temperature gradients near the core are found to be almost completely smoothed out.

Figure 8

Comparison of the contributions to the ICM energy balance from the convection and p dV work, radiation, and anisotropic thermal conduction. The comparisons are for the case of the LCD initial profile in planar geometry, employing CRM simulation physics and a 1-μG turbulent magnetic field. The slice plots are along x at y = 0. The computational mesh consisted of 128² square elements yielding a minimum size of 7.8 kpc for each element.

Open in new tabDownload slide

The diminishing effect of anisotropic thermal conduction is even more evident in the case of a rotational initial magnetic field, which is found to reduce on average the heat flux to the core within the time interval of interest by more than 1 order of magnitude compared to the case of a turbulent initial field (see Fig. 6, left). The evolution of the pertinent average parameters is shown Fig. 9 (left). The rotational magnetic field remains approximately orthogonal to the radial unit vector; the angle between them does not fall below ∼78°. Yet, this has little effect on the evolution of the temperature during this time. This may seem counterintuitive at first since the characteristic times associated with the smoothing of temperature gradients in the ICM by thermal conduction are larger in the rotational-field case as the effective heat flux factor is significantly lower compared to the turbulent field [equation (25) gives τ_q∼ 12 Gyr for f = 0.02]. However, thermal conduction does not drive the temperature in these cases. The characteristic time for temperature changes induced by the pressure work on the fluid may be estimated by

26

yielding

27

For the induced conditions in our numerical experiments with the LCD profiles, τ_p≈ 0.13 Gyr (≪τ_q). Fig. 10 depicts the evolution of the temperature showing the smoothing of the temperature gradients despite the preservation of a highly rotational magnetic field within the time interval of interest. Thus, under the induced conditions p∇·v dominates over thermal conduction.

Figure 9

Evolution of ⁠, 〈f〉, and ⁠, averaged over the area of a disc of radius R = 20 kpc (±Δ). The cases employed CRM simulation physics in planar geometry. Left: results starting with the LCD initial condition for the density and temperature and a rotational initial magnetic field with maximum magnitude of 1 μG. The velocity flow field here fluctuates between radially inward and outward directions after the first expansion at t∼ 0.1 Gyr (note that we plot the absolute value of the radial velocity). Right: results starting with the HCD initial conditions for the density and temperature and a turbulent initial magnetic field with maximum magnitude of 1 μG. The flow field is directed inward after the first expansion at t∼ 0.1 Gyr (see also Fig. 11).

Open in new tabDownload slide

Figure 10

Evolution of the temperature (in eV) and magnetic (unit vector) field for the case of the LCD initial profile and a rotational initial magnetic field with peak magnitude of a 1 μG. CRM simulation physics were used in planar geometry.

Open in new tabDownload slide

The 2D axisymmetric simulations in Section 3.1 with only radiation and thermal conduction demonstrate a strong sensitivity of the solution on the assumed initial core density, with the HCD profile yielding a collapse of the core at values of f for which the LCD profile yielded no cooling. The results of the cosine of the magnetic field angle, heat flux factor, radial velocity and Mach number are depicted in Fig. 9 (right) for the HCD profiles. The hydrodynamics of the flow here are significantly different compared to the LCD case. The induced hydrostatic non-equilibrium state for the ICM encompasses higher energy density in the system but at a reduced pressure gradient. The resulting expansion is weaker compared to the LCD case (10 km s⁻¹ wave front speed in HCD versus 100 km s⁻¹ in LCD) and the gravitational force is able to turn the flow at R = 20 kpc inward relatively fast, within about 0.1–0.2 Gyr, and the plasma beyond this time continues to flow radially inward. This is in contrast to the LCD case with a rotational magnetic field where the velocity field fluctuates between inward and outward directions (note that we plot the absolute value of the radial velocity). The temperature contours at t = 0.16 Gyr are shown in Fig. 11. Thus, the response of the temperature in Fig. 12 and of the non-dimensional quantities in Fig. 9 (right) responded considerably faster than in the LCD case. The solutions with a peak strength for the initial magnetic field of 1 nG show only minor differences with the 1-μG cases.

Figure 11

Temperature (in eV) and velocity (unit vector) field for the case of the HCD initial profile and a turbulent initial magnetic field with peak magnitude of a 1 μG and CRM simulation physics in planar geometry.

Open in new tabDownload slide

Figure 12

Comparisons between the average temperature for the LCD and HCD profiles starting with a turbulent initial magnetic field with maximum magnitude of 1 nG and 1 μG. All cases employed CRM simulation physics in planar geometry.

Open in new tabDownload slide

4 SUMMARY AND DISCUSSION

The works of Zakamska & Narayan (2003), Conroy & Ostriker (2008) and Parrish et al. (2009) have been particularly influential in formulating the different cases presented in this paper. Our 2D axisymmetric simulations reproduce closely the effective heat flux factor of Zakamska & Narayan (2003) under similar assumptions for A2199, when thermal conduction and (bremsstrahlung) radiation are the only physics allowed in the energy balance of the cluster, and when the electron number density and temperature profiles closely resemble the authors’‘best-fitting’ solution; we have termed this solution as the ‘LCD profile’.

Because the LCD profile underestimates the Chandra density observations near the core by approximately 35 per cent (albeit in excellent agreement with the data beyond ∼10 kpc), we performed a second set of simulations with an initial profile that was in agreement with the observed electron density near the cluster centre but overestimated the data at larger distances. Differences in the temperature profiles between the two cases were less than 5 per cent for R < 200 kpc. Because the ratio of characteristic time for thermal conduction over the time for radiative cooling is strongly dependent on the electron density, τ_q/τ_Φ∝n_e² (assuming everything else fixed), the goal of these simulations was to determine what effect such an increase would have on the thermal balance in the cluster. By comparison to the LCD case, we found that the ICM collapsed catastrophically in the HCD case for the same values of the heat flux factor f. Our conclusions from this series of idealized simulations regarding the importance of the core density in the evolution of the ICM are in agreement with the conclusions of Conroy & Ostriker (2008), who illustrated the dependence of the ICM collapse rate on the value of the plasma density at the centre. A turbulent initial magnetic field that was held fixed throughout the simulation also led to the collapse of the ICM for both initial density profiles (LCD and HCD). In this series of simulations the anisotropy of the conduction heat flux was accounted for by invoking the approximations of Braginskii (1965) for the thermal conductivity coefficients. The magnetic field was generated by a method proposed by Ruszkowski & Oh (2010), and was held frozen to emulate an extreme case in which galaxy motion maintains the turbulent field in the absence of any flow dynamics in the ICM.

Since it remains largely infeasible to perform fully consistent cosmological simulations of clusters that account for all pertinent physics and scalelengths, most often in MHD simulations the thermodynamic state and magnetic field topology are prescribed as the ‘initial conditions’ of the ICM cluster and the simulations then proceed to seek of mechanisms that maintain thermal balance in the cluster. Such simulations do not always begin with an ICM in thermal balance but all studies we have reviewed imposed hydrostatic equilibrium at t = 0. Referring to Fig. 1, it is interesting to note that the trend of increasing density near the core (R≤ 10 kpc) based on the observed data appears to be one that crosses the two hydrostatic-equilibrium solutions, namely the LCD and HCD profiles. This suggests that, at least near the core of A2199, deviations from such equilibria may be at play. In fact, Sanders & Fabian (2006) hypothesize the presence of a relatively weak (M∼ 1.5) isothermal shock in this region. The MHD simulations in 2D planar geometry we performed here were intended to serve as a preliminary series of numerical experiments to assess possible deviations from hydrostatic equilibrium in connection with transient AGN activity and/or remnant dynamical activity from the early formation of the galaxy cluster. To span a wide range of magnetic field effects, we employed initial strengths and topologies similar to those imposed by Parrish et al. (2009) and Bogdanović et al. (2009). The 2D results now also form the basis for near-future 3D simulations with mach3.

We were particularly interested in the initial response of the plasma near the core (t < 1.3 Gyr, R = 20 kpc) because radiative cooling and the HBI have been shown to lead to a runaway catastrophe in A2199 in this region in less than 3 Gyr (Parrish et al. 2009; Ruszkowski & Oh 2010). Moreover, it has been proposed that mass and momentum transport from low-speed winds (∼10 000 km s⁻¹) associated with AGN outbursts, can have a dramatic effect on the growth of the central supermassive black hole (Ostriker et al. 2010) and, in turn, on the AGN–ICM coupling. We deliberately imposed in the initial thermodynamic state of the ICM a deviation from hydrostatic equilibrium, thereby producing subsonic expansion flow upon initiation of the simulation. By inducing a non-zero velocity field early in the simulation, we allowed for (thermal) pressure work to contribute to the energy balance of the system. In fact, under the induced conditions, pressure work eventually overwhelmed the deleterious effects of the magnetic field, which in the idealized hydrostatic case would have shut off the main supplier of heat to the region (in the absence of an AGN sources), namely thermal conduction. It is noted that only the scalar (thermal) pressure was included in this work; its importance in the presence of a subsonic velocity field and a magnetic field in this region suggests that the full pressure tensor may have to be accounted for as proposed by Kunz et al. (2010). In all three cases studied with different initial magnetic field, and for both initial density profiles (HCD and LCD), we obtained similar results by 2D MHD simulations about the sustainment of the temperature near the core.

5 CONCLUSIONS

In this paper we presented results from a series of numerical experiments in two dimensions for the A2199 galaxy cluster. The simulations employed a relatively wide combination of initial conditions and simulation physics. It is acknowledged that the 2D physical domains did not permit us to address a wide range of possible physics in clusters that are inherently 3D, and the propositions presented herein shall be re-evaluated in three dimensions. Nevertheless, the number of cases performed in this study yielded constructive insight that now also forms the basis of our near-future simulations with mach3.

We chose A2199 largely because observed data for this cluster exist at a relatively small radius from its centre (R∼ 1 kpc). As a consequence, our idealized 2D axisymmetric simulations aimed to quantify the sensitivity of the ICM thermal balance on the plasma density profiles with less ambiguity compared to other clusters for which such data do not exist. The computed differences in the heat flux factor for the two profiles we used were not insignificant. Therefore, it appears it would be instructive that MHD parametric studies, performed to assess the sensitivity of proposed dominant physics in the ICM and its sustainment/collapse (e.g. HBI, MTI, AGN heating, turbulence stirring), extended their range to include higher values of the plasma density near the cluster centre.

The results of the 2D planar simulations demonstrate that imbalances in the system’s hydrostatic equilibrium, inducing relatively weak flow dynamics, may overcome the effects of the magnetic field on the diffusion of heat from the outer regions of the cluster. This marginalizes the significance of anisotropic thermal conduction over transient flow events in the ICM. Such events may possibly be associated with AGN pulsations at the cluster centre, with characteristic dynamical times that are hundreds of times shorter than the Hubble time.

The authors wish to thank Jeremiah Ostriker and Matt Kunz for their useful comments. We are also grateful to Moustafa T. Chahine for his support through the Innovative Spontaneous Concepts Research and Technology Program. This work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, and made extensive use of the NASA Astrophysics Data System and arXiv.org preprint server.

REFERENCES

Appendix

APPENDIX A: THE FULL MHD EQUATIONS OF mach

A1 Conservation equations

The mach codes solve equations (A1)–(A6), the time-dependent, compressible, single-fluid, multitemperature resistive MHD conservation laws, on a mesh that is composed of arbitrary hexahedral cells.

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

In equations (A1)–(A6), ρ is the mass density of the fluid, v is the fluid velocity, B is the magnetic induction, ε_e/i is the electron/ion specific internal energy, ε_R is the radiation energy density, and is related to the radiation temperature T_R via Stefan’s constant a: ε_R=aT⁴_R. Stefan’s constant is related to the Stefan–Boltzmann constant σ=ac/4, where c is the speed of light in vacuum. T_e/i is the electron/ion temperature, p_e/i is the electron/ion thermal pressure, Q is the artificial numerical compressional viscosity, J is the current density, n_e is the electron number density, C_V is the specific heat, τ_ei is the electron–ion equilibration time, e is the magnitude of the charge of an electron and μ₀ is the permeability of free space. The electrical resistivity tensor is and is the electron/ion thermal conductivity tensor.

Energy is conserved separately for electrons and ions. Models for radiation emission (thin limit), equilibrium conduction (thick limit) as well as non-equilibrium (flux-limited) radiation conduction also exist. The SESAME tables serve as tabular EOS for more than 100 materials. The library is created and maintained by two groups at the Los Alamos National Laboratory (EOS and conductivity tables are the work of T-1; opacity tables are created and maintained by T-4). Although tabular values for the electrical and thermal conductivities are often available, the default models of our mach version use the Spitzer–Braginskii conductivities parallel and perpendicular to the magnetic field separately for electrons and ions. Most of the physical processes are time-advanced in an implicit fashion that allows the Courant condition to be exceeded without exciting numerical instabilities. The spatial derivatives are obtained by integrating over finite volumes. The clean-up procedure applied to B to maintain its divergence at an acceptably small level (solenoidal condition) uses a successive over-relaxation solver to iterate on the in-plane components towards ∇·B = 0 by solving for the scalar potential φ (elliptic correction since a Poisson’s equation is solved), and is based on the method proposed by Brackbill & Barnes (1980).

In the last three decades, mach2 has been employed to study a variety of plasma problems some of which include plasma opening switches (Buff et al. 1987; Degnan et al. 1987), gas and solid density z-pinch implosion physics (Degnan et al. 1995; Shumlak, Hussey & Peterkin 1995), very high power plasma source designs (Mikellides, Turchi & Mikellides 2002b), magnetic nozzles (Mikellides et al. 2002a) and a variety of plasma thrusters (Mikellides, Turchi & Roderick 2000; Turchi et al. 2000; Mikellides & Villarreal 2007).

mach is a three-temperature code, that is, the temperature of the radiation field T_R may be different from T_e/i. In this case, the radiation conduction equation (A5) is solved where both Planck (optically thin limit) χ_p and Rosseland (optically thick limit) χ_r mean opacities must be supplied. The radiation coupling term is in general Φ_eR=acρχ_p(T⁴_e−T⁴_R). Fig. A1 compares the radiation emission rates based on the mach2 tabular opacities and the formula by Rybicki & Lightman (1979) with a Gaunt factor of 1.5 (note the factor ranges 1.1–1.5).

Figure A1

Comparison of the radiation emission rates based on the mach2 tabular opacities and the formula by Rybicki & Lightman (1979) with a Gaunt factor of 1.5.

Open in new tabDownload slide

A2 Anisotropic thermal conduction flux in mach2

In general, the conductive heat flux for electrons in a magnetic field (in the absence the thermoelectric effect) consists of parallel, perpendicular and transverse components as follows:

(A7)

where

(A8)

and ⁠, ⁠. The conductivity coefficients are given by

(A9)

and are the approximations of Braginskii (1965) for the case of a fully ionized plasma, with κ_e0=k_Bp_eτ_e/m_e. Thermal diffusion is solved for implicitly in mach with an iterative Jacobi method. The rate of convergence of the Jacobi iteration may be increased by multigrid methods (e.g. Brandt 1973), which we have employed in the present mach2 simulations.