MAGNETIC FIELD AMPLIFICATION AND EVOLUTION IN TURBULENT COLLISIONLESS MAGNETOHYDRODYNAMICS: AN APPLICATION TO THE INTRACLUSTER MEDIUM

In the absence of the source terms in Equation (1), we recover the standard CGL-MHD model. In this model, first obtained by Chew et al. (1956), the evolution of the pressure components is expressed by

$$\begin{equation} \frac{d }{d t} \left(\frac{p_{\perp }}{\rho B} \right) = 0, \;\;\;\;\; \frac{d }{d t} \left(\frac{p_{\parallel } B^{2}}{\rho ^{3}} \right) = 0, \end{equation} \tag{ 4 }$$

which is also called the double-adiabatic law.

Linear perturbation analysis of the CGL-MHD equations reveals three waves, analogous to the Alfvén, slow, and fast magnetosonic MHD waves. These waves, however, can have imaginary frequencies for sufficiently high degrees of the pressure anisotropy. The corresponding dispersion relations can be found in Kulsrud (1983). For convenience, we reproduce here these relations (using the same notation as in Hau & Wang 2007):

$$\begin{equation} \left(\frac{\omega }{k} \right)^{2}_{a} = \left(\frac{B^{2}}{4 \pi \rho } + \frac{p_{\perp }}{\rho } - \frac{p_{\parallel }}{\rho } \right) \cos ^{2}\theta , \end{equation} \tag{ 5 }$$

$$\begin{equation} \left(\frac{\omega }{k} \right)^{2}_{f,s} = \frac{b \pm \sqrt{b^{2} - 4 c}}{2}, \end{equation} \tag{ 6 }$$

where cos θ = k · B/B (k being the wavevector of the perturbation) and

$$\begin{equation*} b = \frac{B^{2}}{4 \pi \rho } + \frac{2 p_{\perp }}{\rho } + \left(\frac{2 p_{\parallel }}{\rho } - \frac{p_{\perp }}{\rho } \right) \cos ^{2}\theta , \end{equation*}$$

$$\begin{equation*} c = - \cos ^{2}\theta \left[ \left(\frac{3 p_{\parallel }}{\rho } \right)^{2} \cos ^{2}\theta - \frac{3 p_{\parallel }}{\rho } b + \left(\frac{p_{\perp }}{\rho } \right)^{2} \sin ^{2}\theta \right]. \end{equation*}$$

The dispersion relation for the transverse (Alfvén) mode $(\omega / k)^{2}_{a}$ coincides with that obtained from the kinetic theory (in the limit when the Larmor radius goes to zero) and does not change when heat conduction is added to the system (see Kulsrud 1983); the criterion for the instability (named firehose instability), in terms of A = p_⊥/p_∥ and β_∥ = p_∥/(B²/8π), is in this case

$$\begin{equation} A < 1 - 2\beta _{\parallel }^{-1}. \end{equation} \tag{ 7 }$$

However, for the compressible modes $(\omega /k)^{2}_{f,s}$ (which include the mirror unstable modes), the linear dispersion relation of the CGL-MHD equations is known to deviate from the kinetic theory. The mirror instability criterion is

$$\begin{equation} A/6 > 1 + \beta _{\perp }^{-1}, \end{equation} \tag{ 8 }$$

while the one derived from the kinetic theory is

$$\begin{equation} A > 1 + \beta _{\perp }^{-1}, \end{equation} \tag{ 9 }$$

where β_⊥ = p_⊥/(B²/8π) in the last two expressions.

Taking into account the finite Larmor radius effects, Meng et al. (2012a) (see also Hall 1979, 1980, 1981) give the following expressions for the maximum growth rate γ_max (normalized by the ion gyrofrequency Ω_i) of the firehose and kinetic mirror instabilities:

$$\begin{equation} \frac{\gamma _{{\rm max}}}{\Omega _{i}} = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\displaystyle { \frac{1}{2} \frac{(1 - A - 2 \beta _{\parallel }^{-1})}{\sqrt{A-1/4}} } \;\;\; {\rm (firehose),} \\ \displaystyle { \frac{4 \sqrt{2}}{3 \sqrt{5}} \sqrt{A - 1 - \beta _{\perp }^{-1}} } \;\;\; {\rm (mirror),} \end{array}\right. \end{equation} \tag{ 10 }$$

which are achieved for k⁻¹ ∼ l_ci, the ion Larmor radius. These expressions are valid for the case of |A − 1| ≪ 1 and β_{∥, ⊥} ≫ 1.

Following previous works (Denton et al. 1994; see also Pudovkin et al. 1999; Samsonov & Pudovkin 2000; Samsonov et al. 2001; Meng et al. 2012a), whenever the plasma satisfies the firehose (Equation (7)) or kinetic mirror instability criteria (Equation (9)), we impose the following pressure anisotropy relaxation condition:

$$\begin{equation} \left(\frac{\partial p_{\perp }}{\partial t} \right)_{S} = - \frac{1}{2} \left(\frac{\partial p_{\parallel }}{\partial t} \right)_{S} = - \nu _{S} (p_{\perp } - p_{\perp }^{*} ), \end{equation} \tag{ 11 }$$

where $p_{\perp }^{*}$ is the value of p_⊥ for the marginally stable state (which is obtained from the equality in Equations (7) and (9) for each instability and with the conservation of the thermal energy w). The way by which Equation (11) is related to the rate of change of the anisotropy $\dot{A}_{S}$ of Equation (1) is discussed in Section 3.1.

It is not clear yet how the saturation and isotropization timescales are related to the local physical parameters. Some authors claim that the values of ν_S are of the order of the maximum growth rate of each instability γ_max, which in turn is a fraction of the ion Larmor frequency γ_max/Ω_i ∼ 10⁻²–10⁻¹ (see Gary et al. 1997, 1998, 2000). In the ICM, the frequency Ω_i is very large compared to the frequencies that we resolve numerically. This implies that ν_S → ∞ would be a good approximation, or in other words, the relaxation to the marginal values would be instantaneous (which is similar to the "hardwalls" employed in Sharma et al. 2006). However, it is not clear yet whether the extreme low density and weak magnetic fields of the ICM would result in isotropization timescales as fast as these. Therefore, we have also tested, for comparison, finite values for ν_S that are ≪Ω_i (see Section 3.3).

The CGL-MHD model neglects any heat conduction or radiative cooling mechanisms, which is not a realistic assumption for the ICM. Combining Equations (4), we find that

$$\begin{equation} w^{*} = \left[ \left(\frac{B}{B_{0}} \right) A_{0} + \frac{1}{2} \left(\frac{B}{B_{0}} \right)^{-2} \left(\frac{\rho }{\rho _{0}} \right)^{2} \right] \frac{w_{0}^{*}}{\left(A_{0} + 1/2 \right) }, \end{equation} \tag{ 12 }$$

where w* = (p_⊥ + p_∥/2)/ρ, and the subscripts 0 refer to the initial values in the Lagrangian fluid volume.

In the statistically steady state of the turbulence, the constant turbulent dissipation power leads to a secular increasing of the temperature of the gas, which can lead to heat conduction and radiative losses. In order to deal with these effects in a simplified way, we employ a term $\dot{w}$ that relaxes the specific internal energy w* to the initial value $w_{0}^{*}$ at a constant rate ν_th (see a similar approach in Brandenburg et al. 1995):

$$\begin{equation} \dot{w} = - \nu _{{\rm th}}(w^{*} - w_{0}^{*}) \rho. \end{equation} \tag{ 13 }$$

Although simplistic, this approximation is useful for two reasons: (1) it allows the system to dissipate the turbulent power excess; and (2) it helps to relax the local values of w*, which may become artificially high or low in the CGL-MHD formulation without constraints on the anisotropy growth (see discussion in Section 5.2).

Equations (1) were evolved in a 3D Cartesian box employing a modified version of the shock-capturing, second-order Godunov code (Kowal et al. 2011b). The numerical fluxes were calculated using the HLL Riemann solver, with the maximum characteristic speed evaluated from Equations (5) and (6). For the time integration we used the second-order Runge–Kutta method (RK2).

The induction equation was integrated in its "uncurled" form, in an equivalent way to the Constrained Transport method. In the simulations presented in this work, no explicit diffusion term was used, except for model A2 (see Table 1), where an ohmic dissipation term with a diffusivity η ∼ 10⁻⁴ (i.e., of the same order of the numerical diffusivity) was employed in order to prevent eventual negative values of the internal energy w. These can arise because the eventual diffusion of the magnetic energy, especially in the presence of very high spatial frequency instabilities, is not being explicitly taken into account in the energy equation. Nonetheless, numerical tests showed that the introduction of this diffusion term does not cause significant differences in the results for this model.

To prevent negative values of the anisotropy A due to precision errors during the numerical integration, we used an equivalent logarithmic formulation of the last equation in (1), which in the absence of source terms is given by

$$\begin{equation} \frac{\partial }{\partial t} [ \rho \ln (A \rho ^{2}/B^{3} ) ] + \nabla \cdot [ \rho \ln (A \rho ^{2}/B^{3} ) \mathbf {u} ] = 0. \end{equation} \tag{ 14 }$$

The pressure anisotropy relaxation was applied after each sub-time step of the RK2 method, by transforming the conservative variables e and A into the primitive ones p_⊥ and p_∥, calculating their relaxed values through Equation (11) (using the same implicit method as in Meng et al. 2012b)⁹ and then reconstructing back the conservative variables.

The turbulence (represented by the source term f in Equation (1)) is driven by adding a solenoidal velocity field to the domain at the end of each time step. This velocity field is calculated in the Fourier space with a random (but solenoidal) distribution in directions and sharply centered in a chosen value k_turb (being the injection scale l_turb = L/k_turb, where L is the size of the cubic domain). The forcing is approximately delta correlated in time.

The time-step constraint based on the Courant stability criterion δt_C is calculated taking into account both the real and the imaginary characteristic speeds of the linear modes (Equations (5) and (6)).

Another time-step constraint is considered due to the thermal relaxation (Equation (13)). At the end of each time step, we estimate the minimum characteristic time of the thermal relaxation δt_th for the next time step as given by

$$\begin{equation} \delta t_{{\rm th}} = \min \left(\frac{ w }{ \dot{w} } \right), \end{equation} \tag{ 15 }$$

where $\dot{w}$ is the value calculated during the time step and the minimum value is computed over the whole domain.

The next time step is then taken as the minimum between _Cδt_C and _thδt_th where, after performing several tests, we have chosen the factors _C = 0.3 and _th = 0.1.

In the next sections, all the physical quantities are given in code units and can be easily converted in physical units using the reference physical quantities described below.

We arbitrarily choose three representative quantities from which all the other ones can be derived: a length scale l_* (which is given by the computational box side), a density ρ_* (given by the initial ambient density of the system), and a velocity v_* (given by the initial sound speed in most of the models, but Model C3 for which the velocity unit is 0.3v_*; see Table 1). For instance, with such representative quantities the physical timescale is given by the time in code units multiplied by l_*/v_*; the physical energy density is obtained from the energy value in code units times $\rho _{*} v_{*}^{2}$, and so on. The magnetic field in code units is already divided by $\sqrt{4 \pi }$; thus, to obtain the magnetic field in physical units, one has to multiply the value in code units by $v_{*} \sqrt{4 \pi \rho _{*}}$.

Table 1 lists the simulated models and their initial parameters.

In Table 1, V_A0 is the Alfvén speed given by the initial intensity of the magnetic field directed along the x-axis. Initially, the gas pressure is isotropic for all the models with an isothermal sound speed V_S0. The parameter β₀ is the initial ratio between the thermal pressure and the magnetic pressure.

Turbulence was driven considering the same setup in all the models of Table 1. The injection scale is l_turb = 0.4. The power of injected turbulence _turb is kept constant and equal to unity. After t = 1 a fully turbulent flow develops in the system with an rms velocity v_turb close to unity. This implies a turbulent turnover (or cascading) time t_turb ≈ 0.4.

Models A, B, and C in Table 1 are collisionless MHD models with initial moderate, strong, and very small (seed) magnetic fields, respectively. For models A and C the injected turbulence is initially super-Alfvénic, while for models B it is sub-Alfvénic.

Amhd, Bmhd, and Cmhd correspond to collisional MHD models, i.e., have no anisotropy in pressure. The set of equations describing these models is identical to those in Equation (1), but dropping the equation for the evolution of the anisotropy A and replacing the thermal energy by w = 3p/2 (which corresponds to a polytropic gas index 5/3). Their corresponding dispersion relations are those from the usual collisional MHD approach (rather than Equations (5) and (6)).

To better explore the parameter space, we have also considered different values of the pressure anisotropy relaxation rate ν_S. Models with ν_S = ∞ (i.e., with instantaneous relaxation rate) represent conditions for which the relaxation time $\sim \nu _S^{-1}$ is much shorter than the minimum time step δt_min that our numerical simulations are able to solve (δt_min ∼ 10⁻⁶). Previous studies (Gary et al. 1997, 1998, 2000) suggest that the rate ν_S should be of the order of a few percent of Ω_p, the proton gyrofrequency (see Section 5). If we consider typical physical conditions for the ICM, in order to convert the code units into physical units (see Section 3.2), we may take l_* = 100 kpc, v_* = 10⁸ cm s⁻¹, and ρ_* = 10⁻²⁷ g cm⁻³ as characteristic values for the length scale, dynamical velocity, and density of the ICM, respectively. This implies a characteristic timescale t_* ∼ 10¹⁵ s, while for models A in Table 1, the proton Larmor period is τ_cp ∼ 10³ s. Using $\nu _S \sim 10^{-3} \tau _{cp}^{-1}$, we find $\nu _S^{-1} \sim 10^{-9} t_*$. Therefore, the models of Table 1 for which we assumed ν_S = ∞ are very good approximations to the description of the direct effect of plasma instabilities at the large-scale turbulent motions within the ICM. For comparison, we have also run models with no anisotropy relaxation, or ν_S = 0, which thus behave like standard CGL-models.

We notice that the turbulence in the ICM is expected to be trans- or even subsonic, and the plasma β is expected to be high (β ∼ 200). Therefore, models A are possibly more representative of the typical conditions in the ICM.

In the following section we will start by describing the results for models A and B, which have initial finite magnetic fields and therefore reach a nearly steady state turbulent regime relatively rapidly after the injection of turbulence. Then, we will describe models C, which start with seed magnetic fields and therefore undergo a dynamo amplification of field due to the turbulence and take much longer to reach a nearly steady state.

The injected turbulence produces shear and compression in the gas and in the magnetic field. Under the collisionless approximation, according to Equations (4) A∝B³/ρ²; therefore, one should expect that compressions along the magnetic field lines, which keep B constant but make ρ increase, cause a decrease of A, while compressions or shear perpendicular to the magnetic field lines, which make B increase but keep either B/ρ or ρ constant, cause an increase of A. Therefore, even starting with A = 1, parcels of the gas with A ≠ 1 will naturally develop. Inside these parcels, kinetic instabilities can be triggered, which in turn will inhibit the growth of the anisotropy.

Figure 2 presents the distribution of the anisotropy A as a function of β_∥ for the models with moderate initial magnetic field A1, A2, A3, A4, and A5 of Table 1. Model A2 (ν_S = 0) has an A distribution that nearly follows a line with negative inclination in the log–log diagram. This is consistent with the derived A dependence in the CGL models given by $A \propto (\rho /B) \beta _{\parallel }^{-1}$ (when the initial conditions are homogeneous; see Equations (4)). This model attains values of A spanning several orders of magnitude (from 10⁻² to 10³).

Zoom In Zoom Out Reset image size

Model A1 (ν_S = ∞), on the other hand, keeps A close to unity, varying by less than one order of magnitude.

Figure 2 also shows the distribution of A for the A3, A4, and A5 models, which have bounded anisotropy with finite anisotropy relaxation rates ν_S (see Table 1). We see that in these cases, a fraction of the gas has A values out of the stable zone. The model with smaller anisotropy relaxation rate (model A3) obviously presents a larger fraction of gas inside the unstable zones. We also note that the higher the value of β_∥, the larger the linear growth rate of the instabilities and the more gas is inside the unstable regions with A < 1. This is consistent with the CGL trend for which $A \propto \beta _{\parallel }^{-1}$.

Bottom right panel of Figure 2 shows the distribution of A versus β_∥ for the model B1 with strong initial magnetic field (small β₀ = 0.2). We see that, in this regime, model B1 has an A distribution inside the stable zone.

The spatial anisotropy distribution is illustrated in Figure 3 in two-dimensional (2D) maps that depict central slices of A in the XY-plane at the final time step for models A1, A2, A3, A4, A5, and B1. For the CGL model with moderate magnetic field (β₀ = 200), model A2, the A structures are thin and elongated. These small-scale structures probably arise from the fast fluctuations driven by the kinetic instabilities (see Figure 2). For the model with strong magnetic field (small β₀), model B1, the A structures are smoother. They are originated by small-amplitude magnetic fluctuations (Alfvén waves) and also compression modes at the large scales. The map of model A1 also shows thin and elongated structures, but with lengths of the order of the turbulence scale.

Zoom In Zoom Out Reset image size

As an illustration of the spatial distribution of the unstable gas, Figure 4 depicts maps of the maximum growth rate of both the firehose (left column) and the mirror (right column) instabilities given by Equations (10) for the models with moderate initial magnetic field (β₀ = 200) and different anisotropy relaxation rates ν_S.¹⁰ These maximum growth rates are normalized by the initial ion gyrofrequency Ω_i0 and occur for modes with wavelengths of the order of the ion Larmor radius. The first thing to note is that the mirror unstable regions have a larger volume filling factor than the firehose unstable regions for all the models in Figure 4. This is because the regions where the magnetic field is amplified have a large perpendicular pressure, and this happens on most of the turbulent volume. Regions with an excess of parallel pressure arise when the magnetic intensity decays, like in regions with magnetic field reversals (and also in local enhancements of density provided by parallel shocks). The correspondence of the low-intensity magnetic field with firehose unstable regions can be checked directly in model A2 by comparing the maps of Figures 4 and 1. The firehose unstable regions in models A2 and A3 in Figure 4 are small and fragmented, while in models A4 and A5 they are elongated (at lengths of the turbulent injection scale) and very thin (with thickness of the order of the dissipative scales) and are regions with magnetic field reversals and reconnection.

Zoom In Zoom Out Reset image size

Also, from Figure 4 we see that most of the volume of models A2 and A3 is mirror unstable; for models A4 and A5, the mirror unstable regions are elongated but with much larger thickness than in the firehose unstable regions. We must remember that the criterion for the mirror unstable regions in Figure 4 is the kinetic one (Equation (9)) rather than the CGL-MHD criterion (Equation (8); see also Figure 2).

The spatial dimensions of the unstable regions in Figure 4 also reveal the maximum wavelength of the unstable modes that should develop inside the turbulent domain. In the models with finite anisotropy relaxation rate ν_S, the larger the value of ν_S, the smaller the wavelength of the unstable modes. For realistic values of ν_S of the order of γ_max (the maximum frequency of the instabilities), there would have only been unstable modes with wavelengths below the spatial dimensions we can resolve.

The gyrotropic tensor gives the gas a larger (smaller) strength to resist against bending or stretching of the field lines if A > 1 (A < 1). This higher or smaller strength comes from the parallel anisotropic force

$$\begin{equation} f_A = (p_{\parallel } - p_{\perp }) \nabla _{\parallel } \ln {B}, \end{equation} \tag{ 16 }$$

where ∇_∥ ≡ (B/B) · ∇. The relative strength between this anisotropic force and the usual Lorentz curvature force can be estimated from α ≡ (p_∥ − p_⊥)/(B²/4π).

As a measure of the dynamical importance of the anisotropy, we calculated the average value of |α| for all the models of Table 1, and the values are listed in Tables 2, 3, and 4 for models A, B, and C, respectively.

First, let us consider the models with initial moderate magnetic field (β₀ = 200). For models A1, A6, A7, and A8 (ν_S = ∞), the anisotropic force is nondominant: 〈|α|〉 ≈ 0.4. For model A2 (ν_S = 0), on the other hand, the anisotropic force is dominant, with 〈|α|〉 ≈ 5. For models A3, A4, and A5, with finite isotropization rate, the anisotropic force is comparable to the curvature force, being smaller for the higher isotropization rate: 〈|α|〉 ≈ 4 for model A3 (ν_S = 10¹) and 〈|α|〉 ≈ 1.5 for model A5 (ν_S = 10³).

For the model with strong magnetic field (β₀ = 0.2), model B1, the anisotropic force is negligible compared to the Lorentz curvature force: 〈|α|〉 ≈ 0.04.

Figure 5 shows the normalized histograms of log ρ for models A and B of Table 1 having different rates of anisotropy relaxation ν_S. The upper panel shows models with initial moderate magnetic field intensity (β₀ = 200), and the lower panel the model with initial strong magnetic field intensity (β₀ = 0.2). The corresponding collisional MHD models are also shown for comparison.

Zoom In Zoom Out Reset image size

Examining the high β models in the top diagram, we note that all the models with anisotropy relaxation have similar distribution to the collisional model. Model A2, for which the anisotropy relaxation is null, has a much broader distribution, especially in the low-density domain. This difference is due to the presence of strong mirror forces in the A2 model, which expels the gas to outside of high magnetic field intensity regions, causing the formation of low-density zones. Consistently, we can check this effect in the bi-histograms of density versus magnetic field intensity in Figure 10 of the Appendix for model A2 (the bi-histogram for model Amhd is also shown for comparison). The lowest density points are correlated with high-intensity magnetic fields for the model A2.

The bottom panel of Figure 5 indicates that the low-β, strong magnetic field model B1 has density distribution only slightly narrower than the collisional MHD model Bmhd, especially at the high-density region. The slight difference with respect to the collisional model is possibly due to: (1) the sound speed parallel to the field lines is higher in the collisionless models, $c_{\parallel s} = \sqrt{3 p_{\parallel }/\rho }$ for the collisionless model, while for the collisional model $c_s = \sqrt{5 p/3 \rho }$; and (2) in the direction perpendicular to the magnetic field, the fast modes have characteristic speeds higher in the collisionless model: $c_f = \sqrt{B^{2}/4\pi \rho + 2p_{\perp }/\rho }$, while for the MHD model $c_f = \sqrt{B^{2}/4\pi \rho + 5p/3\rho }$. These larger speeds in the anisotropic model imply a larger resistance to compression and therefore smaller density enhancements (at least for our transonic models).

Power spectrum is an important characteristic of turbulence. For MHD turbulence substantial progress has been achieved recently as the Goldreich–Sridhar model has become acceptable. Recent numerical work has tried to resolve the controversies and confirmed the Kolmogorov −5/3 spectrum of Alfvénic turbulence predicted in the model (e.g., Beresnyak & Lazarian 2009, 2010; Beresnyak 2011, 2012b). This spectrum corresponds to the Alfvénic mode of the compressible MHD turbulence (Cho & Lazarian 2002, 2003; Kowal & Lazarian 2010; Beresnyak & Lazarian 2013).

Our goal here is to determine the power spectrum of the turbulence in collisionless plasma in the presence of the feedback of plasma instabilities on scattering.

Figure 6 compares, for different models of Table 1, the power spectra of the velocity (top row), magnetic field (middle row), and density (bottom row). The models starting with moderate magnetic field and β₀ = 200 (A1, A2, A3, A4, A5, Amhd), for which the turbulence is super-Alfvénic, are in the left column, and the models starting with strong magnetic field and β₀ = 0.2 (B1, Bmhd), for which the turbulence is sub-Alfvénic, are in the right column. Each power spectrum is multiplied by the factor k^5/3.

Zoom In Zoom Out Reset image size

The velocity power spectrum P_u(k) for the super-Alfvénic high-β collisional model Amhd (in the left top panel of Figure 6) is consistent with the Kolmogorov slope approximately in the interval 4 < k < 20 and decays quickly for k > 30. The power spectra of the collisionless models A1, A3, A4, and A5 are similar, but show slightly less power in the interval 4 < k < 30. In fact, in Table 2 in the Appendix, we find that the average values of u² for these models are smaller than the model Amhd. Models A3 and A4 evidence more power at the smallest scales, already at the dissipation range. This is due to the acceleration of gas produced by the firehose instability (see Figure 2). Model A2 (ν_S = 0) has a flatter velocity power spectrum than the collisional MHD model Amhd and much more power at the smallest scales. This excess of power comes from the firehose and mirror instabilities and is consistent with the trend reported in the previous sections and also in Kowal et al. (2011b).

The sub-Alfénic velocity power spectrum P_u(k) of the collisional MHD model Bmhd (top right panel in Figure 6) has a narrower interval of wavenumbers consistent with the Kolmogorov slope. The power spectra P_u(k) of the collisionless model B1 are almost identical, which is in agreement with the small dynamical importance of the anisotropy forces compared to the magnetic forces (see Section 4.2).

The power spectrum related to the compressible component of the velocity field P_C(k) is shown in Figure 11 in the Appendix, where it is divided at each wavenumber by the total power of the velocity field. For the high-β models, the ratio P_C(k)/P_u(k) for the collisionless models is similar to that of the collisional MHD model Amhd for almost every wavenumber k and is ≈0.15. For the low-β model, however, the collisionless model has a ratio P_C(k)/P_u(k) slightly higher than that of the collisional MHD model Bmhd for wavenumbers above k ≈ 10. The fractional power in the compressible modes in the interval 2 < k < 20 is smaller compared to the super-Alfvénic (high-β) models, but at larger wavenumbers it becomes higher.

The anisotropy in the structure function of the velocity is shown in Figure 7. The structure function of the velocity $S^{u}_{2}$ is defined by

$$\begin{equation} S^{u}_{2}(l_{\parallel }, l_{\perp }) \equiv \langle | \mathbf {u(r+l)} - \mathbf {u(r)} |^{2} \rangle , \end{equation} \tag{ 17 }$$

where the displacement vector l has the parallel and perpendicular components (relative to the local mean magnetic field) l_∥ and l_⊥, respectively. The local mean magnetic field is defined by (B(r + l) + B(r))/2 (as in Cho et al. 2002 and Zrake & MacFadyen 2012). The GS95 theory predicts an anisotropy scale dependence of the velocity structures (eddies) of the form $l_{\parallel } \propto l_{\perp }^{2/3}$. The axis in Figure 7 is in cell units. The collisional MHD model Amhd is consistent with the GS95 scaling for the interval 10Δ < l_⊥ < 40Δ, where Δ is one cell unit in the computational grid. For the sub-Alfvénic model Bmhd, however, this scaling is less clear, although the anisotropy is clearly seen.

Zoom In Zoom Out Reset image size

The collisionless models A1, A4, and A5 in Figure 7 evidence anisotropy in the velocity structures, which is identical to that of the collisional MHD model Amhd. Models A2 and A3, on the other hand, have more isotropized structures at small values of l. This effect is due to the action of the instabilities and is also observed in the high-β models in Kowal et al. (2011b) for both the firehose and mirror instability regimes.

Using the two-point statistics of the structure function, we can also estimate the velocity power spectrum in the directions parallel or perpendicular to the local mean magnetic field (see Cho et al. 2002 for a detailed discussion on this subject). The scaling laws $S^{u}_{2} (l_{\parallel }, 0) \propto l_{\parallel }^{\zeta }$ and $S^{u}_{2} (0, l_{\perp }) \propto l_{\perp }^{\xi }$ are related to the one-dimensional power spectra in the directions parallel and perpendicular to the local mean magnetic field: $E(k_{\parallel }) \propto k_{\parallel }^{-(1 + \zeta)}$ and $E(k_{\perp }) \propto k_{\perp }^{-(1 + \xi)}$, respectively. For the range of scales corresponding to the wavenumbers 4 < k < 20, we fitted ζ and ξ for the different models. The results are shown in Table 5 (we should observe that this fitting is sensitive to the range of scales chosen). For the model Amhd, our results are close to the expectations of the GS95 theory: $E(k_{\perp }) \propto k_{\perp }^{-5/3}$ and $E(k_{\parallel }) \propto k_{\parallel }^{-2}$ (also in agreement with the values obtained by Cho et al. 2002). Also, at least for the velocity field, the power-law indices of the collisionless models A1–A5 in Table 5 are similar to model Amhd, except for model A2, for which the values of ζ and ξ result in flatter power-law spectra.

The magnetic field power spectra P_B(k) of the collisional MHD models Amhd and Bmhd (middle row of Figure 6) show a power law consistent with the Kolmogorov slope at the same intervals of the velocity power spectra. As in the velocity power spectrum, in the high-β, super-Alfvénic cases, the collisionless models A1, A3, A4, and A5 have similar P_B(k) to the collisional model Amhd (although with slightly less power). Model A2 has a P_B(k) much flatter than that of Amhd and has less power (by a factor of two) at the inertial range interval. In the smallest scales (k > 50), however, its power is above that of the Amhd model (this is also observed for model A3). As in the velocity power spectrum, these small-scale structures are due to the instabilities that are present in this model.

For the sub-Alfvénic, low-β models (B), the magnetic field power spectrum P_B(k) of the collisionless model is again similar to the collisional MHD model Bmhd.

Figure 12 of the Appendix compares P_B(k) and P_u(k) for our models. For the super-Alfvénic, high-β models (A) that are in steady state, the magnetic field power spectrum is in super equipartition with the velocity power spectrum for k > 3 for all models, except the A2 model, which has P_B(k) < P_u(k) for all wavenumbers. Models A3, A4, and A5 show P_B(k)/P_u(k) decreasing values for larger wavenumbers, this effect being more pronounced in model A3, which has a smaller anisotropy relaxation rate. The sub-Alfvénic, low-β collisionless model B1 has the ratio P_B(k)/P_u(k) slightly smaller than unity for all wavenumbers and slightly smaller than the collisional Bmhd model at large k values.

The anisotropy in the structure function for the magnetic field shows a similar trend to the velocity field in all models and is not presented here. Likewise, the density power spectra P_ρ(k) for the super-Alfvénic, high-β models (bottom row in Figure 6) reveal the same trend of the velocity power spectra. For the sub-Alfvénic model, however, the smaller power in the larger scales compared to the collisional MHD model Bmhd is clearly evident, specially in the inertial range. This is consistent with the discussion following the presentation of the density distribution (Section 4.3), which evidenced that the collisionless models resist more to compression than the collisional model (see also Figure 11).

Figure 8 shows the magnetic energy evolution of the models having initially very weak magnetic (seed) field, models C1, C2, C3, C4, and Cmhd of Table 1. The kinetic energy of the models is not shown, but their values are approximately constant in time (after t ≈ 1), and their average values (taken during the last Δt = 10 for each model) 〈E_K〉 are shown in Table 4. For each of these models, Figure 9 shows the power spectrum of the magnetic field, from t = 2 until the final time, for every Δt = 2 (dashed lines). The final magnetic field power spectrum is the continuous line. Also for comparison, the final velocity power spectrum is plotted (dash-dotted line).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The collisional MHD model Cmhd shows an initial exponential growth of the magnetic energy until t ≈ 5. In this interval, the average magnetic energy grows from E_M = 5 × 10⁻⁷ to E_M ∼ 10⁻². After this, a slower (linear) growth rate of the magnetic energy takes place until t ≈ 10, as can be seen in the bottom panel of Figure 8. This is consistent with studies of turbulent dynamo amplification of magnetic fields in collisional plasmas (see, e.g., Cho et al. 2009). At the final times, the magnetic energy achieves the value E_M ≈ 9.0 × 10⁻² which is approximately four times smaller than the average kinetic energy E_K ≈ 0.38 (see Table 4). The bottom panel of Figure 9 shows that the final magnetic field power spectrum is peaked at k ≈ 20 above, which it is in super-equipartition with the velocity power spectrum.

The collisionless model C1 with instantaneous relaxation of the pressure anisotropy (ν_S = ∞) has a turbulent amplification of the magnetic energy very similar to that of the collisional MHD model Cmhd (Figure 8). The initial exponential growth rates are nearly indistinguishable between the two models, but in the linear stage the growth rate is slightly smaller in model C1 (see bottom panel of Figure 8) and also the final value of saturation of the magnetic energy: E_M ≈ 6.2 × 10⁻² (see Table 4). During the initial exponential growth of the magnetic energy, when the plasma still has high values of β, the pressure anisotropy relaxation due to the kinetic instabilities keeps the plasma mostly isotropic, explaining the similar behavior to the collisional MHD model. When β starts to decrease, the anisotropy A can increase (or decrease), spanning a range of A values in the stable zone (as in Figure 2 for model A1). Then, the anisotropic forces can start to have dynamical importance. At the final times, the value of 〈|p_∥ − p_⊥|/(B²/4π)〉, which measures the dynamical importance of the anisotropic forces compared to the Lorentz curvature force (see Section 4.2), is ≈0.5 (Table 4). The magnetic field power spectrum has an identical shape to the model Cmhd, specially in the final time step.

The turbulent dynamo is also tested for a model with a finite anisotropy relaxation rate, model C4, which has ν_S = 10². The growth rate of the magnetic energy (in the exponential and linear phases) is smaller compared to models Cmhd and C1 (top and bottom panels in Figure 8). In this case, the anisotropy A > 1 develops moderately during the magnetic energy amplification and gives the mirror forces some dynamical importance to change the usual collisional MHD dynamics. The value of the magnetic energy at the final time of the simulation is approximately one-third of the value for the collisional MHD model. The final magnetic power spectrum has a shape similar to the collisional MHD model Cmhd, but below the equipartition with the velocity field power spectrum, which has more power at the smallest scales due to the presence of the instabilities (Figure 9).

Model C2, a standard CGL model with no constraints on the growth of pressure anisotropy (ν_S = 0), shows no evidence of a turbulent dynamo amplification of its magnetic energy that saturates at very low values already at t ≈ 5 (Figure 8), when E_M ≈ 6.2 × 10⁻⁶, while the kinetic energy is E_K ≈ 0.32 (see Table 4). The reason is that the anisotropy A increases at the same time that the magnetic field is increased (A∝B³/ρ² in the CGL closure), giving rise to strong mirror forces along the field lines, which increase their resistance against bending or stretching. For this model, 〈|p_∥ − p_⊥|/(B²/4π)〉 ∼ 10⁵, that is, the anisotropic forces dominate over the Lorentz force. The magnetic field power spectrum (top right panel in Figure 9) is similar in shape (but not in intensity) to the Cmhd model, being peaked at k ≈ 40.

The saturated value of the magnetic energy for models without anisotropy relaxation is, nevertheless, sensitive to the initial plasma β. Model C3 is similar to model C2, but starts with a lower sound speed (V_S0 = 0.3), which makes β 10 times smaller (see Table 1). Turbulence is supersonic in this case, rather than transonic. The magnetic energy evolution is similar to that of model C2, but the magnetic energy saturates with a value about two orders of magnitude larger, although the anisotropic forces are still dominant, with 〈|p_∥ − p_⊥|/(B²/4π)〉 ∼ 10⁴ (see Table 4).

Although the electrons have a larger collisional rate than the protons in the ICM (${\sim} \sqrt{m_p/m_e}$), we have assumed in this work, for simplicity, that both species have the same anisotropy in pressure. Also, we assumed them to be in "thermal equilibrium." A more precise approximation would be to consider the electrons only with an isotropic pressure. This would require another equation to evolve the electronic pressure and additional physical ingredients in our model, such as a prescription on how to share the turbulent energy converted into heat at the end of the turbulent cascade or how to quantify the thermalization of the free-energy released by the kinetic instabilities, as well as a description of the cooling for each of the species. The assumption of same temperature and pressure anisotropy for both species has resulted in a force on the collisionless plasma due to the latter that is maximized. Nevertheless, since our results have shown that the dynamics of the turbulence when considering the relaxation of the anisotropy due to the instabilities feedback is similar to that of collisional MHD, we can conclude that if we had considered the electronic pressure to be already isotropic, then this similarity would be even greater.

Another relevant aspect that should be considered in future work regards the fact that the electron thermal speed achieves relativistic values for temperatures ∼10 keV, which are typical in the ICM. Thus, a more consistent calculation would require a relativistic treatment (see, for example, Hazeltine & Mahajan 2002).

Our model considers a thermal relaxation (Equation (13)), which ensures that the average temperature of the domain is maintained nearly constant, despite the continuous dissipation of turbulent power. This simplification allowed us to avoid a detailed description of the radiative cooling and its influence on the temperature anisotropy, even though the rate ν_th = 5 employed in most of our simulations is low enough to not perturb significantly the dispersion relation arising from the CGL-MHD equations. Timescales $\delta t \simeq \nu _{{\rm th}}^{-1} = 0.2$ are much larger than the typical time step of our simulations (∼10⁻⁵). This means that the maximum characteristic speeds calculated via relations (5) and (6) were more than appropriate for the calculation of the fluxes in our numerical scheme (see Section 3.1).

In order to evaluate the effects of the rate of the thermal relaxation on the turbulence statistics, we also performed numerical simulations of three models with different rates ν_th (namely, models A6, A7, and A8 of Table 1). Model A6 has no thermal relaxation (ν_th = 0), while model A7 has a slower rate than model A1, ν_th = 0.5, and both A6 and A7 systems undergo a continuous increase of the temperature as time evolves, which increases β and reduces the sonic Mach number of the turbulence. Model A8, on the other hand, has a faster rate ν_th = 50 and quickly converges to the isothermal limit. Despite different averages and standard deviations in their internal energy, models A6, A7, and A8 presented overall behavior similar to model A1 (see Table 2).

We have also tested models without anisotropy relaxation (not shown here) that employed the CGL-MHD equations of state for calculating the pressure components parallel and perpendicular to the magnetic field (Equations (4) accompanied by homogeneous initial conditions, rather than evolving the last two equations in Equation (1) for the anisotropy A and the internal energy, respectively). Although there are some intrinsic differences due to larger local values of the sound speeds, the overall behavior of these models was qualitatively similar to the models with ν_S = 0 presented here.

In spite of the results above, a more accurate treatment of the energy evolution will be desirable in future work. For instance, as discussed earlier, the lack of a proper treatment for the heat conduction makes the linear behavior of the mirror instability in a fluid description different from the kinetic theory, leading to an overstability of the system. A higher order fluid model reproducing the kinetic linear behavior of the mirror instability (see Equation (9)) would enhance the effects of this instability in the models with finite ν_S, probably producing more small-scale fluctuations compared to the present results (see Figure 4). The effects of the mirror instability on the turbulence statistics have been extensively discussed in Kowal et al. (2011b; see also next section), where a double-isothermal closure was used. This closure is able to reproduce the threshold of the mirror instability given by kinetic derivation.

Kowal et al. (2011b) studied the statistics of the turbulence in collisionless MHD flows assuming fixed parallel and perpendicular temperatures in the so-called double-isothermal approximation, but without taking into account the effects of anisotropy saturation due to the instability feedback. They explored different regimes of turbulence (considering different combinations of sonic and Alfvénic Mach numbers) and initially different (firehose or mirror) unstable regimes. They analyzed the power spectra of the density and velocity, and also the anisotropy of the structure function of these quantities, and found that super-Alfvénic, supersonic turbulence in these double-isothermal collisionless models does not evidence significant differences compared to the collisional-MHD counterpart.

In the case of subsonic models, they have also detected an increase in the density and velocity power spectra at the smallest scales due to the growth of the instabilities at these scales, when compared to the collisional-MHD counterparts. They found elongation of the density and velocity structures along the magnetic field in mirror unstable simulations and isotropization of these structures in the firehose unstable models. In the present study, the closest to their models is the high β, super-Alfvénic model A2, which is without anisotropy relaxation. As in their subsonic sub-Alfvénic mirror unstable case, the instabilities accumulate power in the smallest scales of the density and velocity spectra. However, we should note that the density and velocity structures in our Model A2 become more isotropic at these scales probably because it is in a super-Alfvénic regime.

Our simulations starting with initial seed magnetic field have revealed the crucial role of the pressure anisotropy saturation (due to the mirror instability) for the dynamo turbulent amplification of the magnetic field, which in turn increases the anisotropy A. In our seed field simulations without anisotropy constraints (models C2 and C3), where the mirror forces dominate the dynamics, the turbulent flow is not able to stretch the field lines; therefore, there is no magnetic field amplification. This is consistent with earlier results presented in Santos-Lima et al. (2011) and de Gouveia Dal Pino et al. (2013), and also with the findings of de Lima et al. (2009), where the failure of the turbulent dynamo using a double-isothermal closure for p_⊥ > p_∥ was reported. On the other hand, in model C1 where the pressure anisotropy growth is constrained by the instabilities, there is a dynamo amplification of the magnetic energy until nearly equipartition with the kinetic energy. This result is in agreement with 3D numerical simulations of MRI turbulence performed by Sharma et al. (2006), where a collisionless fluid model taking into account the effects of heat conduction was employed in a shearing box. They have found that the anisotropic stress stabilizes the MRI when no bounds on the anisotropy are considered, making the magnetic lines stiff and avoiding its amplification. When using bounds on the anisotropy, however, they found that the generated MRI is similar (but with corrections of order of unity) to the collisional-MHD case. Sharma et al. (2006), however, did not consider any cooling mechanism, so that the temperature increased continuously in their simulations. Besides, the simulations presented here have substantially larger resolution. Further, they have found that the system overall evolution is nearly insensitive to the adopted threshold values for the anisotropy. We have also found little difference in the turbulence statistics between models with different non-null values of the anisotropy relaxation rate.

Meng et al. (2012a) also employed a collisionless MHD model to investigate the Earth's magnetosphere by means of 3D global simulations. They employed the CGL closure, adding terms to constrain the anisotropy in the ion pressure only (the electronic pressure considered isotropic was neglected in their study). Using real data from the solar wind at the inflow boundary, they compared the outcome of the model in trajectories where data from spacecraft (correlated to the inflow data) were available. Then, they repeated the same calculation, but employing a collisional MHD model. They found better agreement with the collisionless MHD model in the trajectory passing by the bowshock region, where gas is compressed in the direction parallel to the radial magnetic field lines, producing a firehose (A < 1) unstable zone. However, in the trajectory passing by the magnetotail, the simulated data in the collisionless model were not found to be more precise than in the collisional case. In summary, the collisionless MHD description of the magnetosphere seems to differ little from the standard MHD model when the anisotropy is constrained. Even though, they have found quantitative differences in, for example, the thickness of the magnetosheath, which is augmented in the collisionless case, in better agreement with the observations. In the more homogeneous problem discussed here, in a domain with periodic boundaries and isotropic turbulence driving, we have found that both the evolution of the turbulence and the turbulent dynamo growth in the ICM under a collisionless-MHD description accounting for the anisotropy saturation due to the kinetic instability feedback behave similarly (both qualitatively and quantitatively) to the collisional-MHD description.¹¹

Brunetti & Lazarian (2011) appealed to theoretical arguments about the decrease of the effective mean free path and related isotropization of the particle distribution to argue that the collisionless damping of compressible modes will be reduced in the ICM compared to the calculations in earlier papers (Brunetti & Lazarian 2007).¹² Our present calculations do not account for the collisionless damping of compressible motions, but similar to Brunetti & Lazarian (2011), we may argue that this type of damping is not important, at least for the large-scale compressions.

The dynamics of the ICM plasmas is important for understanding most of the ICM physics, including the formation of galaxy clusters and their evolution. The relaxation that we discussed in this paper explains how clusters can have magnetic field generation, as well as turbulent cascade present there. We showed that for sufficiently high rates of isotropization arising from the interaction of particles with magnetic fluctuations induced by plasma instabilities, the collisionless plasma becomes effectively collisional and can be described by an ordinary MHD approach. This can serve as a justification for earlier MHD studies of the ICM dynamics and can motivate new ones.

In general, ICM studies face one major problem. The estimated Reynolds number for the ICM using the Coulomb cross sections is small (∼100 or less), so that one may even question the existence of turbulence in galaxy clusters. This is the problem that we deal with in the present paper and argue that the Reynolds numbers in the ICM may be much larger than the naive estimates above. The difference comes from the dramatic decrease of the mean free path of the particles due to the interaction of ions with fluctuations induced by plasma instabilities. In other words, our study shows that the collisional MHD approach may correctly represent properties of turbulence in the intracluster plasma. In particular, it indicates that MHD turbulence theory may be applicable to a variety of collisionless media. This is a big extension of the domain of applicability of the Goldreich & Sridhar (1995) theory of Alfvénic turbulence.