LOCAL SIMULATIONS OF INSTABILITIES DRIVEN BY COMPOSITION GRADIENTS IN THE ICM

Atmospheres composed of a plasma that is weakly collisional and weakly magnetized have stability properties that differ qualitatively from collisional atmospheres. Instabilities such as the magnetothermal instability (MTI; Balbus 2000, 2001) and the heat-flux-driven buoyancy instability (HBI; Quataert 2008) can arise when there is a gradient in the temperature either parallel or anti-parallel to the gravitational field. These instabilities, which feed off a gradient in temperature, have been extensively studied (Balbus 2000, 2001; Parrish & Stone 2005, 2007; Parrish & Quataert 2008; Parrish et al. 2008, 2009, 2010, 2012a, 2012b; Quataert 2008; Bogdanović et al. 2009; Ruszkowski & Oh 2010; Kunz 2011; McCourt et al. 2011, 2012; Kunz et al. 2012; Latter & Kunz 2012), and they are believed to be important for the understanding of the dynamical evolution of the intracluster medium (ICM) of galaxy clusters.

These studies assumed that the composition of the plasma is uniform, an assumption that might not be appropriate if heavier elements are able to sediment toward the core of the cluster (Fabian & Pringle 1977). In parallel and complementary studies, the long-term evolution of the radial distribution of elements has been studied using one-dimensional models (Fabian & Pringle 1977; Gilfanov & Syunyaev 1984; Chuzhoy & Nusser 2003; Chuzhoy & Loeb 2004; Peng & Nagai 2009; Shtykovskiy & Gilfanov 2010). The ensuing nonuniform composition has been argued to introduce biases in cluster properties as inferred from observations (Markevitch 2007; Peng & Nagai 2009).

While the studies of the MTI and HBI assumed a uniform plasma, the sedimentation models have yet to include magnetic fields. In an attempt to bridge the gap between the different approaches, and with the goal of understanding the long-term evolution of the composition of the ICM, Pessah & Chakraborty (2013) studied the stability properties of weakly collisional atmospheres with gradients in both temperature and composition. They found that gradients in composition, either parallel or anti-parallel to the gravitational field, can trigger instabilities. In a subsequent study, Berlok & Pessah (2015) carried out a comprehensive study using linear mode analysis and showed that these instabilities are expected to render the composition profiles obtained with current sedimentation models unstable, as it was illustrated using the model of Peng & Nagai (2009).

In this paper, we present the first nonlinear, two-dimensional (2D), numerical simulations of the instabilities that feed off a gradient in composition using a modified version of the MHD code Athena (Stone et al. 2008). The instabilities considered are (i) the magneto-thermo-compositional instability (MTCI), which is maximally unstable when the magnetic field is perpendicular to gravity; (ii) the heat- and particle-flux-driven buoyancy instability (HPBI), which is maximally unstable when the magnetic field is parallel to gravity; and (iii) the diffusion modes, which are maximally unstable when the magnetic field is parallel to gravity. These instabilities arise due to the weakly collisional nature of the ICM, which fundamentally changes the transport properties of a plasma. In this regime, where the gyro-radii of the particles are much smaller than the mean free path for particle collisions, the transport of heat, momentum, and particles will be primarily along the magnetic field lines.

The MTCI and HPBI will be present in isothermal atmospheres in which the composition increases with height, while diffusion modes can be present regardless of the direction of the gradient in composition (Pessah & Chakraborty 2013). The linear dispersion relation presented in Berlok & Pessah (2015) is used to compare with the linear evolution of the simulations. We find good agreement, thereby confirming both the linear theory and our numerical method. For the nonlinear evolution of the instabilities we find that the magnetic field inclination goes to roughly 45° independently of whether the magnetic field is initially horizontal (MTCI) or vertical (HPBI). This is contrary to the instabilities driven by temperature gradients, where the average magnetic field becomes almost vertical (horizontal) for an initially horizontal (vertical) magnetic field (Parrish & Stone 2005; Parrish & Quataert 2008). The simple explanation is that the MTCI and HPBI, both of which grow when the composition increases with height, can operate simultaneously. They are therefore driving the average angle in opposite directions, compromising at roughly 45°. The MTI and HBI, being dependent on temperature gradients in opposite directions, cannot grow at the same time, and so they grow unabated by their counterpart. We also find that both types of instabilities cause turbulent mixing of the helium concentration. We conclude that, in the idealized numerical settings that we employ, instabilities driven by the free energy supplied by a gradient in composition saturate by alleviating the gradient and thereby removing the source of free energy.

The rest of the paper is organized as follows: we start out by introducing the equations of kinetic MHD in Section 2 and how they can be solved numerically in Section 3. In Section 4 we demonstrate that the simulations agree with the linear theory for isothermal atmospheres, and we illustrate how the growth rates depend on some of the key parameters of the problem. We also use atmospheres with gradients in both temperature and composition, motivated by the model of Peng & Nagai (2009) and discussed in Berlok & Pessah (2015), to show that the theory and simulations also agree with both gradients present. In Section 5, we consider the nonlinear evolution of the MTCI and HPBI in isothermal atmospheres in order to determine how they saturate. We summarize and outline future work in Section 6.

We consider a fully ionized, weakly magnetized, and weakly collisional plasma consisting of a mixture of hydrogen and helium. We model such a plasma using the set of equations introduced in Pessah & Chakraborty (2013),¹

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {\boldsymbol{v}})}{\partial t}+{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;\left(\rho {\boldsymbol{v}}{\boldsymbol{v}}+{P}_{{\rm{T}}}{\mathsf{I}}-\displaystyle \frac{{B}^{2}}{4\pi }\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}\right)=-{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{\rm{\Pi }}+\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial E}{\partial t}+{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;\left[(E+{P}_{{\rm{T}}}){\boldsymbol{v}}-\displaystyle \frac{{\boldsymbol{B}}({\boldsymbol{B\cdot v}})}{4\pi }\right]\\ &&\quad =-{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{s}}}-{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;({\rm{\Pi }}\;{\boldsymbol{\cdot }}\;{\boldsymbol{v}})+\rho {\boldsymbol{g\cdot v}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}={\rm{\nabla }}\times ({\boldsymbol{v}}\times {\boldsymbol{B}}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial (c\rho )}{\partial t}+{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;(c\rho {\boldsymbol{v}})=-{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{c}}}.\end{eqnarray} \tag{ 5 }$$

In these equations ρ is the mass density, ${\boldsymbol{v}}$ is the fluid velocity, ${\boldsymbol{B}}$ is the magnetic field with direction $\hat{{\boldsymbol{b}}}=({b}_{x},0,{b}_{z}),{\boldsymbol{g}}=(0,0,-g)$ is the gravitational acceleration, and ${\mathsf{I}}$ is the identity matrix. The total pressure is ${P}_{{\rm{T}}}=P+{B}^{2}/8\pi$, where P is the thermal pressure and the total energy density, E, is

$$\begin{eqnarray}&&E=\displaystyle \frac{1}{2}\rho {v}^{2}+\displaystyle \frac{{B}^{2}}{8\pi }+\displaystyle \frac{P}{\gamma -1},\end{eqnarray} \tag{ 6 }$$

where $\gamma =5/3$ is the adiabatic index.

The composition of the plasma, c, is defined to be the ratio of the helium density to the total gas density

$$\begin{eqnarray}&&c\equiv \displaystyle \frac{{\rho }_{{\rm{He}}}}{{\rho }_{{\rm{H}}}+{\rho }_{{\rm{He}}}}=\displaystyle \frac{{\rho }_{{\rm{He}}}}{\rho },\end{eqnarray} \tag{ 7 }$$

and the associated mean molecular weight, μ, is given by

$$\begin{eqnarray}&&\mu =\displaystyle \frac{4}{8-5c},\end{eqnarray} \tag{ 8 }$$

for a completely ionized plasma consisting of helium and hydrogen. The mean molecular weight can modify the dynamics of the plasma through the equation of state

$$\begin{eqnarray}&&P=\displaystyle \frac{\rho {k}_{{\rm{B}}}T}{\mu {m}_{{\rm{H}}}},\end{eqnarray} \tag{ 9 }$$

where ${k}_{{\rm{B}}}$ is Boltzmann's constant, T is the temperature, and ${m}_{{\rm{H}}}$ is the proton mass.

We consider the plasma to be influenced by three different nonideal effects: Braginskii viscosity, which arises due to differences in pressure parallel (${p}_{\parallel }$) and perpendicular (${p}_{\perp }$) to the magnetic field, described by the viscosity tensor (Braginskii 1965)

$$\begin{eqnarray}&&{\rm{\Pi }}=-3\rho {\nu }_{\parallel }\left(\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}-\displaystyle \frac{1}{3}{\boldsymbol{I}}\right)\left(\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}-\displaystyle \frac{1}{3}{\boldsymbol{I}}\right):{\rm{\nabla }}{\boldsymbol{v}},\end{eqnarray} \tag{ 10 }$$

anisotropic heat conduction described by the heat flux (Spitzer 1962; Braginskii 1965),

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{{\rm{s}}}=-{\chi }_{\parallel }\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}T,\end{eqnarray} \tag{ 11 }$$

and anisotropic diffusion of composition described by the composition flux (Bahcall & Loeb 1990),

$$\begin{eqnarray}&&{{\boldsymbol{Q}}}_{{\rm{c}}}=-D\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}c.\end{eqnarray} \tag{ 12 }$$

The transport coefficients for Braginskii viscosity (${\nu }_{\parallel }$), heat conductivity (${\chi }_{\parallel }$), and diffusion of composition (D) depend on the temperature, density, and composition of the plasma. The dependences are given by Equations (64)–(66) in Berlok & Pessah (2015). Finally, we define the thermal velocity, ${v}_{{\rm{th}}}=\sqrt{P/\rho }$, and the plasma-β given by $\beta =8\pi P/{B}^{2}=2{v}_{{\rm{th}}}^{2}/{v}_{{\rm{A}}}^{2},$ where ${v}_{{\rm{A}}}=B/\sqrt{4\pi \rho }$ is the Alfvén velocity.²

The equations of kinetic MHD, Equations (1)–(5), are solved using a modified version of the conservative MHD code Athena (Stone et al. 2008). The algorithms used in Athena are described in Gardiner & Stone (2005) and Stone & Gardiner (2009), and a description of the implementation of anisotropic thermal conduction and Braginskii viscosity can be found in Parrish & Stone (2005) and Parrish et al. (2012b), respectively.

In order to carry out the numerical simulations of interest, we have modified Athena to include a spatially varying mean molecular weight, μ. This is done by using the inbuilt method for adding a passive scalar, defined by a spatially varying concentration, c, and then making it active by using the value of c when calculating the temperature used in the heat conduction module. Furthermore, we implemented a module that takes account of diffusion of helium by using operator splitting. This module has been built by following the same approach employed in the heat conduction module that is already present in the current publicly available version of Athena (Parrish & Stone 2005; Sharma & Hammett 2007). Our implementation allows for nonconstant values of the parameters ${\nu }_{\parallel },{\chi }_{\parallel }$, and D through user-defined functions. This feature is, however, not used in this work, as we employ a local approximation and thus treat these parameters as constants. The diffusion terms are solved explicitly, which can make the time-step constraint on viscosity, thermal conduction, and diffusion of helium very restrictive. In order to circumvent this, we use subcycling, which we limit to a maximum of 10 steps per MHD step (Kunz et al. 2012).

3.1. Plane-parallel Atmosphere with Gradients in Temperature and Composition

In this section, we introduce the two different atmospheres used as initial conditions in the simulations. The atmospheres considered are plane-parallel, i.e., all quantities are constant along a horizontal slice, perpendicular to gravity. The atmosphere is assumed to be composed of an ideal gas, characterized by the equation of state given by Equation (9), and is assumed to be in hydrostatic equilibrium, i.e.,³

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial z}=-g\rho .\end{eqnarray} \tag{ 13 }$$

3.1.1. Isothermal Atmosphere with a Composition Gradient in the Absence of Particle Diffusion

The simplest atmosphere we use is inspired by the original numerical work on the MTI (Parrish & Stone 2005). We consider an isothermal atmosphere with $T={T}_{0}$ and

$$\begin{eqnarray}&&P={P}_{0}{\left(1-\displaystyle \frac{z}{3{H}_{0}}\right)}^{3},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\mu ={\mu }_{0}{\left(1-\displaystyle \frac{z}{3{H}_{0}}\right)}^{-1}.\end{eqnarray} \tag{ 15 }$$

where ${P}_{0},{T}_{0}$, and ${\mu }_{0}$ are the values of the pressure, temperature, and mean molecular weight at z = 0, respectively, and H₀ is the scale height

$$\begin{eqnarray}&&{H}_{0}=\displaystyle \frac{{k}_{{\rm{B}}}{T}_{0}}{{\mu }_{0}{m}_{{\rm{H}}}g}.\end{eqnarray} \tag{ 16 }$$

The density can be determined using Equation (9).

This isothermal atmosphere is used for simulations of the linear regime of the MTCI in Section 4.1 and the linear regime of the HPBI in Section 4.2. It is also used for simulations of the nonlinear regime of the MTCI and HPBI in Section 5. The magnetic field can have any orientation as long as D = 0. The structure of this atmosphere is, however, not in equilibrium if $D\ne 0$ and ${b}_{z}\ne 0$. In that case, we will have to use a more sophisticated atmosphere, which we introduce next.

3.1.2. Atmosphere with Thermal and Composition Gradients

Steady state requires that the divergence of the heat and particle fluxes vanish, i.e.,

$$\begin{eqnarray}&&{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{s}}}=0,\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{c}}}=0.\end{eqnarray} \tag{ 18 }$$

Both conditions are trivially satisfied if b_z = 0 and D = 0. If, however, ${b}_{z}\ne 0$ and $D\ne 0$, these requirements can still be met by simple atmospheric models if ${\chi }_{\parallel }$ and D do not depend on z. Such an assumption is reasonable for the local simulations that we will consider, where the height of the box, L_z, satisfies the criterion ${L}_{z}\ll {H}_{0}$. When there is both a gradient in temperature T and mean molecular weight μ, the requirements that ${\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{s}}}=0$ and ${\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{c}}}=0$ can be integrated to yield

$$\begin{eqnarray}&&T(z)={T}_{0}+{s}_{{\rm{T}}}z,\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&c(z)={c}_{0}+{s}_{{\rm{c}}}z,\end{eqnarray} \tag{ 20 }$$

where ${s}_{{\rm{T}}}=({T}_{{\rm{Z}}}-{T}_{0})/{L}_{{\rm{Z}}}$ and ${s}_{{\rm{c}}}=({c}_{{\rm{Z}}}-{c}_{0})/{L}_{z}$ are the constant slopes in temperature and composition. Here T₀ (${T}_{{\rm{Z}}}$) is the temperature at the bottom (top) of the box and c₀ (${c}_{{\rm{Z}}}$) is the helium mass concentration at the bottom (top) of the box.

The pressure is found by solving Equation (13), leading to

$$\begin{eqnarray}&&P(z)={P}_{0}{\left(\displaystyle \frac{T(z)\mu (z)}{{T}_{0}{\mu }_{0}}\right)}^{\alpha },\end{eqnarray} \tag{ 21 }$$

where $\mu (z)$ is related to c(z) by Equation (8) and the constant α is given by

$$\begin{eqnarray}&&\alpha =-\displaystyle \frac{{T}_{0}}{{H}_{0}}\displaystyle \frac{4}{4{s}_{{\rm{T}}}+5{\mu }_{0}{T}_{0}{s}_{{\rm{c}}}}.\end{eqnarray} \tag{ 22 }$$

This solution for the pressure profile of the atmosphere is replaced with a simple exponential atmosphere, $P(z)={P}_{0}\mathrm{exp}(-z/{H}_{0})$, with scale height H₀ if ${s}_{{\rm{T}}}={s}_{{\rm{c}}}=0$.

We use this model atmosphere to perform simulations of modes driven by diffusion in Section 4.3. These modes are unstable when there is a vertical gradient in composition, a nonzero vertical component of the magnetic field, ${b}_{z}\ne 0$, and anisotropic diffusion of helium, D $\ne$ 0. We also use this atmosphere in Section 4.4 for simulations of the linear regime of the MTCI and the HBPI with gradients in both temperature and composition.

3.2. Boundary Conditions

The equations are made dimensionless by scaling the density with ${\rho }_{0}$, distances with H₀, and velocities with the thermal velocity ${v}_{{\rm{th,0}}}$. The magnetic field strength B₀ is found from the dimensionless parameter ${\beta }_{0}$. Here the subscript "0" denotes the value at the bottom of the computational domain, z = 0. With this convention, the unit of time is ${H}_{0}/{v}_{\mathrm{th},0}$, temperature is scaled with ${T}_{0},\mu$ is scaled with ${\mu }_{0}$, pressure, as well as energy density, is scaled with ${P}_{0}={\rho }_{0}{v}_{\mathrm{th},0}^{2}$, and the value of g is unity. As a consequence, the coefficient for anisotropic heat conduction, ${\chi }_{\parallel }$, is scaled with ${\rho }_{0}{v}_{\mathrm{th},0}^{3}{H}_{0}/{T}_{0}$, and the coefficients for Braginskii viscosity$,{\nu }_{\parallel }$, and anisotropic diffusion of composition, D, are both scaled with ${v}_{\mathrm{th},0}{H}_{0}$.

We begin by comparing the simulations with the linear theory. In order to do so, we use the quasi-periodic boundary conditions described in the previous section and in Appendix B. A Cartesian box of size $[0,{L}_{x}]\times [0,{L}_{z}]$ with ${L}_{x}={L}_{z}=0.1$ and a resolution of 64 × 64 is used in all simulations unless otherwise noted. An overview of the simulations of the linear regime can be found in Table 1.

Table 1.  Simulations of the Linear Regime Using Quasi-periodic Boundary Conditions

Simulation	(${n}_{x},{n}_{z}$)	θ	${\beta }_{0}$	${\chi }_{\parallel }$	${\nu }_{\parallel }$	D	Resolution	Figure
MTCI_chi	(1, 0)	0°	2 × 108	⋯	0	0	64 × 64	3(a)
MTCI_B	$(1,0)$	0°	⋯	3 × 10−4	0	0	64 × 64	3(b)
HPBI_nu	$(1,1)$	90°	2 × 108	10−4	⋯	0	64 × 64	4(a)
HPBI_n	(⋯, ⋯)	90°	2 × 106	10−4	0	0	256 × 256	4(b)
D-mode_D	$(1,1)$	90°	2 × 108	10−3	0	⋯	256 × 256	5(a)
D-mode_nu	$(1,1)$	90°	2 × 108	10−3	⋯	10−3	64 × 64	5(b)
MTCI_ICMa	(⋯, 0)	0°	2 × 106	$1.4\;\times \;{10}^{-2}$	$4.0\;\times \;{10}^{-4}$	0	256 × 32	6(a)
HPBI_ICMa	(⋯, ⋯)	90°	2 × 106	$2.7\;\times \;{10}^{-4}$	$4.5\;\times \;{10}^{-6}$	0	256 × 256	6(b)

Note. Each row represents a series of simulations, where the ellipses denote that the associated parameter is being varied.

^aUsing gradients in both temperature and composition.

Download table as:  ASCIITypeset image

The instabilities are excited by seeding a given mode, with components ($\delta {v}_{x},\delta {v}_{z},\delta {B}_{x},\delta {B}_{z},\delta \rho /\rho ,\delta T/T,\delta \mu /\mu$), as derived by solving the eigenvalue system associated with the dispersion relation introduced in Pessah & Chakraborty (2013) and Berlok & Pessah (2015). We set the overall mode amplitude by enforcing $\delta \rho /\rho ={10}^{-4}$, so that the velocity perturbation is subsonic (Parrish & Stone 2005). The amplitudes of the other components are fixed by the solution to the linear eigenvalue problem, which predicts that unstable modes grow exponentially as $\mathrm{exp}(\sigma t)$ while the ratio of their components remains constant in time.

We begin by considering D = 0 and the hydrostatic atmosphere given in Section 3.1.1, which has $d\mathrm{ln}\mu /d\mathrm{ln}P=-1/3$ and $d\mathrm{ln}T/d\mathrm{ln}P=0$. This atmosphere is unstable regardless of whether the magnetic field is oriented horizontally (MTCI) or vertically (HPBI), as described in Berlok & Pessah (2015).

4.1. The Magneto-thermo-compositional Instability

When the magnetic field is perpendicular to gravity, the general dispersion relation, Equation (13) in Berlok & Pessah (2015), reduces to

$$\begin{eqnarray}&&{\sigma }^{2}\approx -g\displaystyle \frac{d\mathrm{ln}(T/\mu )}{{dz}}\displaystyle \frac{{k}_{x}^{2}+{k}_{y}^{2}}{{k}^{2}},\end{eqnarray} \tag{ 23 }$$

in the limit of fast heat conduction and weak magnetic field. When $\mu$ increases with height and the atmosphere is isothermal, we have $\sigma \gt 0$. This is the instability known as the MTCI (Pessah & Chakraborty 2013). In order to excite a single MTCI mode, we use a perturbation of the form⁴ ${k}_{z}=0$ and ${k}_{x}=2\pi /{L}_{x}$. We are interested in a direct visual comparison of the spatial dependence of the perturbations in the simulations and the one expected from the linear theory. In order to illustrate this, we consider a setting with ${\chi }_{\parallel }=3\times {10}^{-4}$ and ${\beta }_{0}=2\times {10}^{8}$. In Figure 1, we show the values of the perturbations (green crosses) $\delta {v}_{z},\delta {B}_{z}/B,\delta \mu /\mu$, and $\delta T/T$ as a function of the x-coordinate. The data slices are drawn at a fixed height, $z={L}_{z}/2$ at the time t = 5 in dimensionless units. The numerical results show good agreement with the analytical results shown with blue solid lines.

Zoom In Zoom Out Reset image size

In order to calculate the growth rate of the mode, we perform an exponential fit to the time evolution of the box average of the absolute value of any of the perturbed quantities, which are shown in Figure 2. As expected from the local linear mode analysis, the amplitudes of the various components of the perturbation grow exponentially at the same rate.

Zoom In Zoom Out Reset image size

The growth rate of the MTCI depends on, among other things, the value of the heat conductivity, ${\chi }_{\parallel }$, and the initial magnetic field strength, B₀. In order to illustrate this dependence, and at the same time test our modification to the code, we perform a parameter study. In the left panel of Figure 3, we show how the growth rate increases with the value of the heat conductivity, ${\chi }_{\parallel }$. This is to be expected because the MTCI is driven by heat transfer along magnetic field lines. In the right panel of Figure 3, we show how the growth rate decreases with the value of ${\beta }_{0}^{-1}$. The explanation for this behavior is that magnetic tension tends to stabilize the MTCI (Berlok & Pessah 2015). Magnetic tension has stabilizing effects in the limit ${\omega }_{{\rm{A}}}\gg {\omega }_{{\rm{dyn}}}$, where ${\omega }_{{\rm{A}}}={k}_{\parallel }{v}_{{\rm{A}}}$ and ${\omega }_{{\rm{dyn}}}=\sqrt{g/{H}_{0}}$. In dimensionless units, this requirement can be written as $2{k}_{\parallel }^{2}\gg \beta$. From this estimate, the growth rates shown in the right panel of Figure 3 should be negligible when ${\beta }_{0}^{-1}\gg {10}^{-4}$. The simulations and the solution to the dispersion relation show that the growth rates are already inhibited by magnetic tension at lower values of ${\beta }_{0}^{-1}$. These examples were generated by running 10 simulations at a modest resolution (64 × 64). At this resolution the growth rates match to within a percent of the values expected from linear theory.

Zoom In Zoom Out Reset image size

4.2. The Heat- and Particle-flux-driven Buoyancy Instability

When the magnetic field is parallel to gravity, the general dispersion relation reduces to

$$\begin{eqnarray}&&{\sigma }^{2}\approx g\displaystyle \frac{d\mathrm{ln}(T\mu )}{{dz}}\displaystyle \frac{{k}_{x}^{2}+{k}_{y}^{2}}{{k}^{2}},\end{eqnarray} \tag{ 24 }$$

in the limit of fast heat conduction and weak magnetic field. The isothermal atmosphere where μ increases with height, which we considered in the previous section, is therefore also unstable when the magnetic field is vertical. In this case, the instability has been termed the HPBI (Pessah & Chakraborty 2013).

In this example, we include Braginskii viscosity, which inhibits the growth rate by damping perturbations perpendicular to the magnetic field. Braginskii viscosity can be important for the HPBI (Berlok & Pessah 2015). The mechanism is similar to the mechanism described by Kunz (2011) for the HBI. In order to excite a single mode of the HPBI, we use a perturbation with wavenumbers ${k}_{x}=2\pi {n}_{x}/{L}_{x}$ and ${k}_{z}=2\pi {n}_{z}/{L}_{z}$, where $n={n}_{x}={n}_{z}$ is the mode number. We show the growth rate as a function of the Braginskii viscosity coefficient, ${\nu }_{\parallel }$ in the left panel of Figure 4. These simulations used a fixed value of ${\chi }_{\parallel }={10}^{-4},n=1$, and a numerical resolution of 64 × 64. As expected, the growth rate indeed decreases with increasing value of viscosity ${\nu }_{\parallel }$.

Zoom In Zoom Out Reset image size

The second dependence we study for the HPBI is the one on the mode number, n. High wavenumbers require higher numerical resolution in order to be resolved, and we use a resolution of 256 × 256 for these simulations. For the sake of simplicity, Braginskii viscosity is not included in these simulations. The result is shown in the right panel of Figure 4. The growth rate increases for increasing wavenumber because small wavelength perturbations have a shorter timescale for heat conduction. When the wavelength is too short, magnetic field tension renders the modes stable. A naive estimate, using $2{k}_{\parallel }^{2}\gg \beta$, suggests that this should happen when $n\gg 16$, but the exact solution to the dispersion relation shows that the instability is quenched already when n = 7. Using such simulations, we can directly see the cutoff in unstable wavenumbers resulting from magnetic field tension (as in this case) or viscosity (not shown here).

4.3. Modes Driven by Diffusion

One of the interesting findings of Pessah & Chakraborty (2013) is that there are instabilities that are driven by particle diffusion. This means that even though the equilibrium is stable according to Equation (24), the fact that $D\ne 0$ makes the equilibrium unstable. In order to study these unstable modes, we assume, for simplicity, an isothermal atmosphere with an initially vertical magnetic field.

In this case, as explained in Section 3.1, an equilibrium configuration needs to fulfill ${\rm{\nabla }}\cdot {{\boldsymbol{Q}}}_{{\rm{c}}}=0$, and so we consider the atmosphere given in Section 3.1.2 as an initial condition. According to Equation (24), this configuration is unstable to the HPBI for an isothermal atmosphere when the helium concentration increases vertically. If instead the helium concentration decreases with height, the atmosphere is stable in the absence of anisotropic particle diffusion. Choosing the slope in composition to be ${s}_{{\rm{c}}}=-0.01$, we do not observe any instabilities in the simulation when D = 0. The situation changes dramatically, turning unstable when $D\ne 0$. The growth rates found in such simulations are compared with the predictions from the linear theory in Figure 5. Since the modes are driven by diffusion of helium, we expect the growth rate to increase with the value of D (left panel). The modes have a damped growth rate when Braginskii viscosity is included. We observe a decrease in the growth rate with increasing ${\nu }_{\parallel }$, in agreement with the solution to the dispersion relation (right panel).

Zoom In Zoom Out Reset image size

4.4. Gradients in Temperature and Composition

Having tested the case of isothermal atmospheres, we now consider a more general situation where both ${dT}/{dz}\ne 0$ and $d\mu /{dz}\ne 0$. In order to work with sensible values for these gradients, we consider the models in Peng & Nagai (2009), who analyzed the long-term evolution of the concentration of helium in a one-dimensional setting by solving a coupled set of Burgers' equations for a multicomponent plasma in the absence of a magnetic field. Berlok & Pessah (2015) analyzed the stability of the Peng & Nagai (2009) model by focusing on local regions, characterized by fixed temperature and composition gradients, and modeling these as a plane-parallel atmosphere.

In this section, we present local simulations with gradients in temperature and composition estimated at $r/{r}_{500}=0.02$ and $r/{r}_{500}=0.5$ with ${r}_{500}=1.63$ Mpc in the Peng & Nagai (2009) model. These are the locations that were analyzed in Sections 6.6 and 6.4 in Berlok & Pessah (2015), indicated with a C and an A in Figure 8 in that paper. These two locations correspond to the inner region where the temperature and composition increase with radius and the outer region where the temperature and composition decrease with radius. At these radii, the values for the logarithmic gradients are $d\mathrm{ln}T/d\mathrm{ln}P=-0.4$ and $d\mathrm{ln}\mu /d\mathrm{ln}P=-0.13$ at $r/{r}_{500}=0.02$ and $d\mathrm{ln}T/d\mathrm{ln}P=0.16$ and $d\mathrm{ln}\mu /d\mathrm{ln}P=0.05$ at $r/{r}_{500}=0.5$.

We use the equilibrium derived in Section 3.1.2 with values taken from the model of Peng & Nagai (2009), ${T}_{0}=6.6\;{\rm{keV}}$ (${T}_{0}=9.5$ keV) and ${c}_{0}=0.56$ (${c}_{0}=0.29$) for the inner (outer) region. The computational domain is ${L}_{x}={L}_{z}={H}_{0}/10$ where ${H}_{0}=50$ kpc for the inner region and ${L}_{x}=10{L}_{z}={H}_{0}$ where ${H}_{0}=0.35$ Mpc for the outer region. The gradients in composition and temperature are set such that the dimensionless values of $d\mathrm{ln}T/d\mathrm{ln}P$ and $d\mathrm{ln}\mu /d\mathrm{ln}P$ in the plane-parallel atmosphere agree with the values in the model of Peng & Nagai (2009). They are given by ${s}_{{\rm{c}}}=2.9$ Mpc⁻¹ (${s}_{{\rm{c}}}=-0.19$ Mpc⁻¹) and ${s}_{{\rm{T}}}=0.058\;{\rm{keV}}$ kpc⁻¹ (${s}_{{\rm{T}}}=-4.3\;{\rm{keV}}$ Mpc⁻¹) for the inner (outer) region.

The values for ${\nu }_{\parallel },{\chi }_{\parallel }$, and D are calculated from the model of Peng & Nagai (2009) as explained in the appendix of Berlok & Pessah (2015). The dimensionless values are so large that high-resolution numerical simulations become very computationally expensive. As this is a test, we have arbitrarily reduced the values by a factor of 100 in the simulations. We use a value of $\beta =2\times {10}^{6}$ for both sets of simulations and adopt a resolution of 32 × 256 (MTCI) and 256 × 256 (HPBI). Some of the details of the simulations are listed in Table 1 with the names HPBI_ICM and MTCI_ICM. The growth rates also depend on the wavenumbers, k_x and k_z. For the HPBI (in the inner region) we take $k={k}_{x}={k}_{z}$ and investigate growth rate as a function of k. For the MTCI (in the outer region) we take k_z = 0 and investigate the growth rate as a function of k_x. The results are shown in Figure 6, with the growth rates of the HPBI in the left panel and the growth rates of the MTCI in the right panel. An estimate shows that the HPBI should be suppressed by magnetic tension for $k/2\pi \gg 3.4\;{{\rm{kpc}}}^{-1}$ and the MTCI should be suppressed for $k/2\pi \gg 450\;{{\rm{Mpc}}}^{-1}$. The growth rates are in units of 50 and 280 Myr, respectively. Therefore, in physical units, the maximum growth rates in these simulations are ${\sigma }_{\mathrm{max}}=6.4\;{\mathrm{Gyr}}^{-1}$ for the HBPI and ${\sigma }_{\mathrm{max}}=1.2\;{\mathrm{Gyr}}^{-1}$ for the MTCI.

Zoom In Zoom Out Reset image size

In order to study the nonlinear evolution of the MTCI and the HPBI, we use the reflective boundaries described in Appendix B. We use the isothermal atmospheres presented in Section 3.1.1 and seed both velocity components with Gaussian noise with a standard deviation of 10⁻⁴. The simulations are run without Braginskii viscosity or anisotropic diffusion of helium, but anisotropic heat conduction is accounted for with a value of ${\chi }_{\parallel }=5\times {10}^{-4}$. We use a value of 2 × 10⁸ for the plasma-β. An overview of the simulations of the nonlinear regime can be found in Table 2.

Table 2.  Overview of the Simulations of the Nonlinear Regime Using the Reflective Boundary Conditions

Simulation	θ	${\beta }_{0}$	${\chi }_{\parallel }$	${\nu }_{\parallel }$	D	Resolution	Figure
MTCI256	0°	2 × 108	5 × 10−4	0	0	256 × 256	7(a), 8, 10
HPBI128	90°	2 × 108	5 × 10−4	0	0	128 × 128	9
HPBI256	90°	2 × 108	5 × 10−4	0	0	256 × 256	9
HPBI512	90°	2 × 108	5 × 10−4	0	0	512 × 512	7(b), 9, 10

Download table as:  ASCIITypeset image

We start out by studying the evolution of the MTCI, i.e., we consider an atmosphere threaded by a horizontal magnetic field. The subsequent evolution of the magnetic field and the plasma composition is illustrated in the upper panel of Figure 7. In this figure, it is evident that the MTCI is able to mix the helium content and to completely rearrange the initially ordered magnetic field. The resulting growths in kinetic and magnetic energy densities are shown, respectively, in the left and right panels of Figure 8. The kinetic energies associated with the two velocity components are roughly in equipartition throughout the simulation, i.e., $\langle \rho {v}_{x}^{2}\rangle \approx \langle \rho {v}_{z}^{2}\rangle$ with $\langle \rho {v}_{z}^{2}\rangle$ always larger but never exceeding $\langle \rho {v}_{x}^{2}\rangle$ by more than an order of magnitude. The exponential phase of the instability ends at $t\approx 30$. After this point in time, both the kinetic and magnetic energies saturate, with the former exceeding the latter by two orders of magnitude. In spite of the fact that $\langle {B}_{z}^{2}\rangle$ vanishes initially, by the end of the simulation the energies associated with the two magnetic field components are roughly in equipartition with $\langle {B}_{z}^{2}\rangle$ larger than $\langle {B}_{x}^{2}\rangle$ by a factor of $\approx 2$, with $\langle {B}_{x}^{2}\rangle$ having grown by a factor of $\approx 8$ with respect to its initial value.

Zoom In Zoom Out Reset image size

We now consider the evolution of the HPBI. The setup is essentially the same, but the initial magnetic field is now vertical.⁵ The evolution of the HPBI is illustrated in the lower panel of Figure 7 with a resolution of 512 × 512. The initial vertical magnetic field is rearranged by the HPBI, and, as for the MTCI, the helium content is mixed by the action of the instability. The HPBI leads to growth in the magnetic and kinetic energy densities. In order to asses whether this growth is numerically converged, we have also run simulations at resolutions of 128 × 128 and 256 × 256. We show the evolution of $\langle {B}_{x}^{2}\rangle /8\pi$ and $\langle {B}_{z}^{2}\rangle /8\pi$ for the three different numerical resolutions in Figure 9. We observe that the instability leads to exponential growth followed by saturation in both $\langle {B}_{x}^{2}\rangle /8\pi$ and $\langle {B}_{z}^{2}\rangle /8\pi$. While the growth rate increases with increasing resolution, the values in the saturated state agree quite well.

Zoom In Zoom Out Reset image size

It is also of interest to understand how the magnetic field changes from being initially vertical to having a large horizontal component because the magnetic field inclination has consequences for heat transport along the vertical direction of the box. Such studies have been done for both the MTI (Parrish & Stone 2005, 2007) and the HBI (Parrish & Quataert 2008). These studies were motivated by a need to understand the cooling flow problem of galaxy clusters (Fabian 1994), and whether magnetic fields could alleviate this problem. While the MTI could potentially increase heat transport toward the core by making the magnetic field be preferentially in the radial direction (Parrish et al. 2008), the HBI has been shown to lead to core insulation by driving the magnetic field to be perpendicular to the radial direction (Parrish & Quataert 2008; Bogdanović et al. 2009; Parrish et al. 2009), which would exacerbate the cooling flow problem.

In Figure 10, we show the average magnetic field inclination as a function of time for the simulations of the MTCI and the HPBI. The average inclination saturates to a value of approximately $\theta \approx 45^\circ$ for both the simulations. This behavior is qualitatively different from the behavior of the magnetic field inclination for the MTI and the HBI. The difference can be explained in the following way. The MTI, which is maximally unstable when the magnetic field is horizontal, has been found to drive the saturated magnetic field to be roughly vertical (Parrish & Stone 2007). The HBI, which is maximally unstable when the magnetic field is vertical, drives the magnetic field to be roughly horizontal (Parrish & Quataert 2008). These instabilities depend on gradients in temperature that have opposite directions, and so they cannot be present at the same time. On the other hand, both the MTCI and the HPBI require a mean molecular weight that increases with height, and so they can both be present at the same time. This feature of the MTCI and the HPBI was discussed in Berlok & Pessah (2015); see especially Figure 4 in that paper. The interpretation of the left panel of Figure 10 is therefore that the MTCI aims at driving the magnetic field angle toward 90° while the HPBI aims at driving the magnetic field angle toward 0°. In the end, they reach a compromise at roughly 45°.

Zoom In Zoom Out Reset image size

The helium mass concentration, c, dramatically changes, and the initial gradient is diminished by the instability as time progresses. This is illustrated in the lower panel of Figure 10 for both the MTCI and the HPBI. As explained in the introduction, gradients in composition can introduce biases in key cluster parameters. We are therefore interested in understanding whether such gradients, if initially present, will be robust. The simulations presented here are heavily idealized, among many reasons because the gas is assumed to be initially isothermal and the simulations are local. Nevertheless, these simulations serve as a proof of principle that gradients in composition can indeed be altered by turbulent mixing induced by plasma instabilities. Future work, using realistic gradients for temperature and composition as well as transport coefficients, should allow us to understand whether such mixing can occur on timescales relevant for galaxy clusters.

In this paper we have introduced a modified version of Athena (Stone et al. 2008) for performing kinetic MHD simulations of weakly collisional plasmas with nonuniform composition. We have employed this modified code to perform the first simulations of the MTCI, the HPBI, and the diffusion modes introduced in Pessah & Chakraborty (2013). The set of simulations, aimed at investigating the linear evolution of these instabilities, served as a test for both the modification to Athena and the local linear mode analysis in Pessah & Chakraborty (2013) and Berlok & Pessah (2015).

The simulations of weakly collisional, isothermal atmospheres with a gradient in helium presented in Section 5 showed that the plasma instabilities, feeding off gradients in composition, can induce turbulent mixing of the helium content. This conclusion is valid for compositions that increase in the direction anti-parallel to gravity, regardless of whether the initial magnetic field is parallel or perpendicular to the direction of gravity. In the saturated state, the magnetic field components in the x and z directions have roughly the same average energy, but the energies are a factor of 10 higher for the HPBI than for the MTCI. The kinetic energy components are also roughly in equipartition. In both cases, the instabilities saturate by driving the average magnetic field inclination to roughly 45°. This effect seems to open the possibility of alleviating the core insulation observed in previous homogeneous simulations of the HBI, provided that the global cluster dynamics were to allow for an increase in the mean molecular weight with radius in the inner region, as envisioned by current (one-dimensional, unmagnetized) helium sedimentation models (Peng & Nagai 2009).

The simulations of the nonlinear regime of the MTCI and the HBPI presented in this paper considered an isothermal atmosphere as the equilibrium background. It would be an improvement to use the model of Peng & Nagai (2009) to determine the gradients in both temperature and composition. This would provide insight into the saturation of instabilities in potentially more realistic scenarios where the dynamical evolution is determined by the simultaneous effects of both gradients. Before proceeding with this endeavor, there are, however, a few issues that would be desirable to address, as we detail below.

It was found in Berlok & Pessah (2015) that the HPBI, if present in the inner regions of the ICM model of Peng & Nagai (2009), will have its fastest growth rates at wavelengths that are longer than the scale height of the atmosphere. This conclusion is similar to what was found for the HBI in Kunz (2011). Neither local linear theory nor local simulations will therefore capture the physics of the HPBI in the inner region of the ICM. This implies that both a quasi-global theory and simulations are needed in order to study the influence of the possible gradient in composition on the dynamics of the inner region of the ICM.

Local simulations of the MTI have been shown to underestimate the turbulence (McCourt et al. 2011), and boundary effects can also modify the conclusions from local simulations. The solution to this problem for the MTI has been to sandwich the unstable region between stable layers, thereby isolating it from the boundaries (Parrish & Stone 2005, 2007; Kunz et al. 2012). A similar approach seems reasonable for the MTCI.

Other complications stem from the fact that pressure anisotropies, shown to be important for the evolution of the MTI and HBI (Kunz 2011; Kunz et al. 2012), can give rise to microscale instabilities such as the firehose and mirror instability (Schekochihin & Cowley 2006; Schekochihin et al. 2010). These small-scale instabilities are only excited once the pressure anisotropy grows beyond $| {p}_{\parallel }-{p}_{\perp }| /P\gtrsim {\beta }^{-1}$. They do not appear in the tests of the linear regime of the MTCI, HPBI, and the diffusion modes presented in Section 4 because we terminate the simulations before the stability criterion is violated. They are not present in the simulations of the nonlinear regime in Section 5 because we take the pressure to be isotropic in these simulations (no Braginskii viscosity). The problem with microscale instabilities is that they are not correctly described by the framework of kinetic MHD (Schekochihin et al. 2005), an issue that will need to be addressed for simulations of the nonlinear evolution of the MTCI and the HPBI when Braginskii viscosity is included (see Kunz et al. 2012 for a discussion of these issues in the context of the MTI and the HBI).

All this being said, our study suggests that, at least in the idealized settings that we considered here, gradients in composition are able to drive turbulent mixing of the composition in weakly collisional, magnetized plasmas. This motivates future work on the generation and sustainment of both temperature and composition gradients in galaxy clusters and their potential influence on the global dynamics of the ICM. We envision that the modified version of Athena that we developed will be a useful asset in this context. In order to model more realistically the physics of the ICM, future improvements could include extending the simulations to three dimensions and adding optically thin cooling in order to study the cooling flow problem. Furthermore, the equations of kinetic MHD, as embodied in Equations (1)–(5), cannot account for the slow sedimentation of helium that is the core feature in the model of Peng & Nagai (2009). An extension of the framework of kinetic MHD to include this effect would allow us to self-consistently include sedimentation in the simulations (Bahcall & Loeb 1990; Berlok & Pessah 2015) and study the effects of the instabilities described in this paper in a dynamic, slowly varying background.

We acknowledge useful discussions with Daisuke Nagai, Matthew Kunz, Prateek Sharma, Ellen Zweibel, and Ian Parrish during the 3rd ICM Theory and Computation Workshop held at the Niels Bohr Institute in 2014. We thank Sagar Chakraborty, Colin McNally, Gareth C. Murphy, and Henrik Latter for valuable discussions and comments. We are grateful to Oliver Gressel for suggesting using the quasi-periodic boundary conditions for testing the linear theory and to Tobias Heinemann for aid in rendering magnetic field lines. We also thank the anonymous referee for a number of useful suggestions that helped improve the manuscript. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007–2013) under ERC grant agreement 306614. T.B. also acknowledges support provided by a Lørup Scholar Stipend, and M.E.P. also acknowledges support from the Young Investigator Programme of the Villum Foundation.

In this appendix, we describe the numerical methods used in this paper. We use the publicly available MHD code Athena, which solves the conservative form of the MHD equations. The algorithms used are described in Gardiner & Stone (2005) and Stone & Gardiner (2009), and the implementation of Athena along with tests is described in detail in Stone et al. (2008). Athena is a finite-volume code, which uses the Godunov method. We use the directionally unsplit corner transport upwind method along with constrained transport (CTU + CT), which is the recommended setting. We furthermore use the anisotropic heat conduction module that was implemented in Athena by Parrish & Stone (2005) using operator splitting.

Appendix A explains the implementation of a spatially varying mean molecular weight, μ, the anisotropic diffusion of helium and tests cases, in Athena. In Appendix B we discuss in detail the boundary conditions used in the simulations.

Let us consider the equation describing the evolution of the helium mass concentration, $c={\rho }_{{\rm{He}}}/\rho$, given by⁶

$$\begin{eqnarray}&&\displaystyle \frac{\partial c}{\partial t}+({\boldsymbol{v}}\;{\boldsymbol{\cdot }}\;{\rm{\nabla }})c=-{\boldsymbol{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{c}}}.\end{eqnarray} \tag{ 25 }$$

Athena has an option for adding passive scalars, which we use for adding the helium mass concentration. This option turns on an extra equation,

$$\begin{eqnarray}&&\displaystyle \frac{\partial (\rho {c}_{{\rm{n}}})}{\partial t}+{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;(\rho {c}_{{\rm{n}}}{\boldsymbol{v}})=0,\end{eqnarray} \tag{ 26 }$$

where ρ is the total density and ${c}_{{\rm{n}}}$ is the mass concentration of the nth scalar. We only add a single scalar, namely, the helium mass concentration, c. This built-in function takes care of the Lagrangian part of Equation (25). The diffusion term is then solved using a finite-difference scheme and operator splitting.

Anisotropic diffusion of helium is described by the right-hand side of Equation (25), which, when ${\boldsymbol{v}}=0$, reduces to

$$\begin{eqnarray}&&\displaystyle \frac{\partial c}{\partial t}=-{\boldsymbol{\nabla }}\;{\boldsymbol{\cdot }}\;{{\boldsymbol{Q}}}_{{\rm{c}}}=D{\rm{\nabla }}\;{\boldsymbol{\cdot }}\;(\hat{{\boldsymbol{b}}}\hat{{\boldsymbol{b}}}\;{\boldsymbol{\cdot }}\;{\rm{\nabla }}c).\end{eqnarray} \tag{ 27 }$$

The composition flux for anisotropic helium diffusion has the same form as the heat flux for anisotropic heat conduction, as seen by comparing Equations (12) and (11), respectively. We can therefore use the same method to calculate the two physically different anisotropic fluxes. The original implementation of anisotropic heat conduction was done by Parrish & Stone (2005) using an asymmetric finite-difference scheme (Sharma & Hammett 2007; van Es et al. 2014).

Nonideal effects are computationally expensive because they are generally described by parabolic operators that cannot be added to the hyperbolic fluxes used in the Godunov scheme. The parabolic operators can be shown to have a very prohibitive time-step constraint (Durran 2010) for heat conduction and concentration diffusion as given by, respectively,

$$\begin{eqnarray}&&{\rm{\Delta }}{t}_{{\chi }_{\parallel }}\lt \displaystyle \frac{b}{{\kappa }_{\parallel }}\displaystyle \frac{{({\rm{\Delta }}x)}^{2}}{\gamma -1},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\qquad {\rm{\Delta }}{t}_{D}\lt b\displaystyle \frac{{({\rm{\Delta }}x)}^{2}}{D}.\end{eqnarray} \tag{ 29 }$$

Here ${\kappa }_{\parallel }={\chi }_{\parallel }T/P$ is the heat diffusivity, the parameter b is $b=1/2,1/4,1/6$ in one, two, and three dimensions, respectively, and ${\rm{\Delta }}x$ is the grid size.

The Courant number, C, is defined to be the ratio of the applied time step to the allowed time step. We use C = 0.4 in all our simulations. Because ${\rm{\Delta }}{t}_{{\rm{MHD}}}\propto {\rm{\Delta }}x$, the very prohibitive constraints on the time step for the parabolic operators will generally lead to ${\rm{\Delta }}{t}_{{\rm{MHD}}}\gg {\rm{\Delta }}{t}_{\chi }\sim {\rm{\Delta }}{t}_{D}$. In order to partially circumvent this problem, we use subcycling, taking up to 10 diffusion steps for each MHD step, as suggested in Kunz et al. (2012).

Sharma & Hammett (2007) found that the finite-difference approximation can lead to unphysical behavior with diffusion in the wrong direction. In the context of heat diffusion this problem can lead to negative temperatures and therefore an imaginary sound speed. The same problem arises when considering helium diffusion, and we use van Leer limiters on the derivatives to circumvent it (Sharma & Hammett 2007).

The publicly available version of Athena works with constant viscosity, ${\nu }_{\parallel }$, and heat diffusivity, ${\kappa }_{\parallel }={\chi }_{\parallel }T/P$. However, these coefficients, as well as the diffusion coefficient D, do in general depend on temperature, density, and composition; see, for instance, the discussion in the appendix of Berlok & Pessah (2015). Accounting for this dependence is not crucial in local simulations, but it becomes essential in global simulations. We have modified Athena to use spatially varying coefficients by using a harmonic average of the coefficients (Sharma & Hammett 2007). This makes the time step computed from Equations (28) and (29) spatially dependent. We therefore calculate the time step at each cell and use the minimum value. This implementation will be useful in future global studies.

A.1. Tests of the Implementation of Anisotropic Diffusion

In order to verify the implementation of anisotropic diffusion of helium, we performed three different test problems with a known analytical solution. These tests were carried out with the MHD solver turned off.

A.1.1. One-dimensional Diffusion

We consider the diffusion of a step function as in Rasera & Chandran (2008) using a one-dimensional grid with 100 cells on the domain $x=[0,1]$ with D = 1 and run the simulation up to t = 0.0028. The analytical solution to the diffusion of a step function is

$$\begin{eqnarray}&&c(x,t)={c}_{0}+\displaystyle \frac{{\rm{\Delta }}c}{2}{\rm{erf}}\left(\pm \displaystyle \frac{x-{x}_{0}}{\sqrt{4{Dt}}}\right),\end{eqnarray} \tag{ 30 }$$

where ${c}_{0}=3/2$ and ${\rm{\Delta }}c=1$. The "+" sign is used with ${x}_{0}=0.25$ for $x\lt 0.5$ and the "−" sign is used with ${x}_{0}=0.75$ for $x\gt 0.5$. The numerical result matches the analytical solution, as seen in Figure 11, implying that the method works well in one dimension.

Zoom In Zoom Out Reset image size

A.1.2. Diffusion of a 2D Gaussian

A more challenging test can be posed by considering the magnetic field to be inclined at an angle, θ, with respect to the grid. We consider an initially isotropic, 2D Gaussian distribution of helium diffusing out along an inclined magnetic field. The analytical solution is⁷

$$\begin{eqnarray}c(x,y,t) & = & \displaystyle \frac{1}{2\pi a(t){a}_{0}}\mathrm{exp}\left\{-\displaystyle \frac{{(x\mathrm{cos}\theta +y\mathrm{sin}\theta )}^{2}}{2a{(t)}^{2}}\right\}\\ & & \times \mathrm{exp}\left\{-\displaystyle \frac{{(y\mathrm{cos}\theta -x\mathrm{sin}\theta )}^{2}}{2{a}_{0}^{2}}\right\},\end{eqnarray} \tag{ 31 }$$

where $a{(t)}^{2}={a}_{0}^{2}+2{Dt}$ and a₀ is the initial standard deviation of the Gaussian.

The computational domain is a $[-1,1]\times [-1,1]$ Cartesian box. We use ${a}_{0}=1/8$ and D = 0.001. The errors at t = 4 are compared to the analytical solution in Figure 12. In the left panel the L₂ errors are shown as a function of the magnetic field inclination and resolution. These errors are smallest when $\theta =0$ or $\theta =\pi /2$, corresponding to the grid and the magnetic field being aligned. In the right panel we show that the solution for $\theta =40^\circ$ converges as ${L}_{2}\propto {({\rm{\Delta }}x)}^{m},$ where ${\rm{\Delta }}x$ is the (uniform) grid spacing and m = 1.9 is the order of convergence.

Zoom In Zoom Out Reset image size

A.1.3. Diffusion of a High-concentration Patch in a Circular Magnetic Field

The final and most challenging test that we carry out for anisotropic transport was introduced in Parrish & Stone (2005). We consider a Cartesian box of size $[-1,1]\times [-1,1]$ with a patch with higher concentration c, specifically⁸

$$\begin{eqnarray} & & \\ & c\quad = & \left\{\begin{array}{cc}12 & {\rm{if}}\quad 0.5\lt r\lt 0.7\quad {\rm{and}}-\pi /12\lt \theta \lt \pi /12,\\ 10 & {\rm{otherwise}}.\end{array}\right.\end{eqnarray} \tag{ 32 }$$

The density is uniform with $\rho =1$, and the magnetic field is circular. In order to ensure ${\rm{\nabla }}\;{\boldsymbol{\cdot }}\;{\boldsymbol{B}}=0$, the magnetic field was initialized with a vector potential satisfying ${\rm{\nabla }}\times {\boldsymbol{A}}={\boldsymbol{B}}$.

We considered the value D = 0.01 and run the simulation until t = 200. The overconcentration diffuses out along the magnetic field lines, as observed in Figure 13. We have run this test problem with the same resolutions as Sharma & Hammett (2007) and obtain the exact same values quoted there for the error norms associated with the resolutions 200 × 200 and 400 × 400. For instance, for a resolution of 200 × 200, we obtain ${L}_{1}=0.0264,{L}_{2}=0.0407$, ${L}_{\infty }=0.0928,{c}_{{\rm{min}}}=10$, and ${c}_{{\rm{max}}}=10.1016$ at t = 200 as stated in Sharma & Hammett (2007).

Zoom In Zoom Out Reset image size

It is evident from Figure 13 that, even though only anisotropic diffusion is explicitly turned on, there is still a small amount of numerical, perpendicular diffusion. This is undesired in simulations of instabilities because isotropic diffusion will lower the growth rates or even quench the instabilities. This was investigated by Parrish & Stone (2005) for the MTI, who found that it is, however, insensitive to perpendicular diffusion provided that ${\chi }_{\perp }/{\chi }_{\parallel }\lt {10}^{-3}$. The fact that we find the correct analytical growth rates for all simulations discussed in Section 4 shows that the numerical perpendicular diffusion is not a problem for our present purposes.

In the horizontal direction we use the periodic boundary conditions that Athena provides as a standard option. In the vertical direction, boundary conditions that maintain hydrostatic equilibrium are required. In this section we describe the conventional, reflective boundary conditions, as well as the quasi-periodic boundary conditions alluded to in Section 3.2.

B.1. Reflective Boundary Conditions

Our implementation of reflective boundary conditions follows the description in Zingale et al. (2002). Hydrostatic equilibrium requires that Equation (13) is satisfied. This requirement can be approximated by

$$\begin{eqnarray}&&{P}_{i}-{P}_{i+1}=\displaystyle \frac{{ag}\;{\rm{\Delta }}z}{12}(5{\rho }_{i+1}+8{\rho }_{i}-2{\rho }_{i-1}),\end{eqnarray} \tag{ 33 }$$

where a = 1 ($a=-1$ ) at the top (bottom) of the domain. In this notation, i refers to cell i and $i+1$ refers to one cell further up (down) when a = 1 ($a=-1$). This equation is then solved for ${\rho }_{i+1}$ using that

$$\begin{eqnarray}&&{P}_{i+1}=\displaystyle \frac{{\rho }_{i+1}{T}_{i+1}}{{\mu }_{i+1}},\end{eqnarray} \tag{ 34 }$$

along with an assumption on ${\mu }_{i+1}$ and ${T}_{i+1}$. Either one can assume ${\mu }_{i}={\mu }_{i+1}$ and ${T}_{i}={T}_{i+1}$, or one can prescribe the values at the boundaries to be equal to their initial values, i.e., ${\mu }_{i+1}={\mu }_{0}$ and ${T}_{i+1}={T}_{0}$. In the case where the mean molecular weight μ is not included, Parrish & Stone (2005) refer to these boundary conditions as adiabatic and conductive, respectively. A combination of these two boundary conditions is also possible (i.e., fixing μ and varying T or vice versa). We have implemented all four combinations but will only discuss the conducting boundary conditions (${\mu }_{i+1}={\mu }_{0}$ and ${T}_{i+1}={T}_{0}$) in the following, since these are the boundary conditions used in Section 5.

Solving Equations (33) and (34) for ${\rho }_{i+1}$ and ${P}_{i+1}$, we find

$$\begin{eqnarray}&&{\rho }_{i+1}=\displaystyle \frac{{P}_{i}+\alpha (8{\rho }_{i}-{\rho }_{i-1})}{{T}_{0}/{\mu }_{0}-5\alpha }\end{eqnarray} \tag{ 35 }$$

and ${P}_{i+1}={\rho }_{i+1}{T}_{0}/{\mu }_{0}$, where

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{a{\rm{\Delta }}z\;g}{12}.\end{eqnarray} \tag{ 36 }$$

These relations are used to calculate the density and pressure of the four ghost zones at the top and bottom of the computational domain. At the same time, velocity is reflected symmetrically around z = 0 and $z={L}_{z}$. The magnetic field components are also mirrored. In the case of an initially vertical magnetic field, we let the B_x component change sign, whereas in the case of an initially horizontal magnetic field, we let the B_z component change sign. This forces the field to remain vertical (horizontal) at the boundary in the case of an initially vertical (horizontal) field.

Athena uses the Godunov scheme, which is known not to be optimal at maintaining hydrostatic equilibrium (Zingale et al. 2002). The reason is that the pressure term in the momentum equation is not solved simultaneously with the gravity term. There are ways to modify a Godunov scheme such that this problem is circumvented; see, for instance, Zingale et al. (2002) and Käppeli & Mishra (2014). We use a high numerical resolution and a low Courant number (C = 0.4) in order to maintain hydrostatic equilibrium as well as possible. The minimum amplitude we can use for perturbations in ${\boldsymbol{v}}$ is, however, limited by the numerical noise caused by the inability of the code to perfectly maintain hydrostatic equilibrium, in agreement with the findings of Parrish et al. (2012b).

B.2. Quasi-periodic Boundary Conditions

A key assumption in standard local linear mode analysis, such as presented in Pessah & Chakraborty (2013) and Berlok & Pessah (2015), is that the perturbations have the spatial dependence $\mathrm{exp}(i{\boldsymbol{k}}\;{\boldsymbol{\cdot }}\;{\boldsymbol{x}})$. This assumption is not fulfilled for the reflective boundary conditions, and it is thus impossible to cleanly excite a single eigenmode, the problem being that the boundary conditions excite other modes in an uncontrolled way. We originally realized this problem when we studied the HBI, but it persists in the case of the HPBI and its diffusive variant. The problem is not present for the MTI and the MTCI in the case k_z = 0 because the boundaries are periodic in x.

We have developed special boundary conditions that are consistent with the assumptions used in the local mode analysis. One of the key assumptions here is that the perturbed quantities $\delta {v}_{x},\delta {v}_{z},\delta {B}_{x}$, $\delta {B}_{z},\delta \rho /\rho ,\delta \mu /\mu$, and $\delta T/T$ are periodic. In the following, the values outside the computational domain (the ghost zones) are denoted by a subscript g, and the values on the inside are denoted by a subscript i. The subscript "eq" refers to the value of the equilibrium background (as given in Sections 3.1.1 or 3.1.2). The mapping from interior to ghost zones ($i\to g$) is the same as for periodic boundary conditions. Instead of directly mapping the interior values to the ghost zones, we let the ghost zones depend on the change in the interior values with respect to the equilibrium background. The quasi-periodic boundary conditions are then defined as

$$\begin{eqnarray}&&{\rho }_{g}={\rho }_{g,{\rm{eq}}}\left(1+\displaystyle \frac{{\rho }_{i}-{\rho }_{i,{\rm{eq}}}}{{\rho }_{i,{\rm{eq}}}}\right),\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{T}_{g}={T}_{g,{\rm{eq}}}\left(1+\displaystyle \frac{{T}_{i}-{T}_{i,{\rm{eq}}}}{{T}_{i,{\rm{eq}}}}\right),\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{\mu }_{g}={\mu }_{g,{\rm{eq}}}\left(1+\displaystyle \frac{{\mu }_{i}-{\mu }_{i,{\rm{eq}}}}{{\mu }_{i,{\rm{eq}}}}\right),\end{eqnarray} \tag{ 39 }$$

with the pressure given by ${P}_{g}={\rho }_{g}{T}_{g}/{\mu }_{g}$. The equilibrium magnetic field and velocity do not have a gradient, and so their boundary conditions are simply periodic, i.e., ${{\boldsymbol{v}}}_{g}={{\boldsymbol{v}}}_{i}$ and ${{\boldsymbol{B}}}_{g}={{\boldsymbol{B}}}_{i}$. These are the boundary conditions we used in Section 4.