Turbulence in collisionless plasmas: statistical analysis from numerical simulations with pressure anisotropy

In the past years we have experienced an increasing interest in understanding of the physical properties of collisionless plasmas, mostly because of the large number of astrophysical environments, e.g. the intracluster medium (ICM), containing magnetic fields which are strong enough to be coupled with the ionized gas and characterized by densities sufficiently low to prevent the pressure isotropization with respect to the magnetic line direction. Under these conditions a new class of kinetic instabilities arises, such as firehose and mirror ones, which were extensively studied in the literature. Their role in the turbulence evolution and cascade process in the presence of pressure anisotropy, however, is still unclear. In this work we present the first statistical analysis of turbulence in collisionless plasmas using three dimensional double isothermal magnetohydrodynamical with the Chew-Goldberger-Low closure (CGL-MHD) numerical simulations. We study models with different initial conditions to account for the firehose and mirror instabilities and to obtain different turbulent regimes. We study the probability distribution functions, spectra, structure functions and anisotropy of density and velocity fluctuations. The results indicate that in some cases the instabilities significantly modifies the statistical properties of turbulence and even though preliminary and restricted to very specific conditions, show that the physical properties of turbulence in collisionless plasmas, as those found in the ICM, may be very different from what has been largely believed. Implications can range from interchange of energies to cosmic rays acceleration.


Introduction
Magnetized and low density (weakly collisional) plasmas are known to present anisotropic pressures with respect to the magnetic field orientation, which can survive considerably long times compared to the dynamical timescales of certain systems. In astrophysical environments such pressure anisotropy may be generated by several different processes, such as kinetic pressure of cosmic rays, supernovae explosions, stellar winds, or anisotropic turbulent motions (see Quest and Shapiro 1996).
Under certain conditions gyrotropic plasmas give rise to new wave modes and instabilities, which cannot be studied by standard isotropic magnetohydrodynamic (MHD) model (Hasegawa 1969, Wang and Hau 2003, Passot and Sulem 2006. For instance, Hau and Wang (2007) showed that gyrotropic magnetohydrodynamic equations closed by the Chew-Goldberger-Low laws (CGL-MHD) lead to a positive density n versus magnetic field strength B correlations for the slow magnetosonic mode under certain conditions in contradiction to the standard MHD model. It was commonly believed that the detected absence of slow modes in several astrophysical sites was related to the strong damping of these waves. However, in the collisionless plasmas it could be explained by wrong identification of the positive n − B correlations for fast magnetosonic modes.
The pressure anisotropies give rise to plasma instabilities depending on the anisotropy ratio, e.g., p > p ⊥ and p ⊥ > p for the firehose and mirror instabilities, respectively. They are responsible for the growth of the magnetic energy and acceleration of particles. The predictions of the CGL-MHD model, including new plasma instabilities, are also important in weakly magnetized environments, since even a weak magnetic field is enough to change the motion of the charged particles and therefore increase the pressure anisotropy.
As an example of the possible applications Sharma et al. (2003) and Sharma et al. (2006) showed the importance of the collisionless plasma approach in protostellar disks. They studied the role of the kinetic instabilities in the magneto-rotational instability (MRI) showing that the transport of angular momentum in the disk may be efficiently increased under certain circumstances.
The intracluster medium (ICM) is possibly the most suitable environment for studies of gyrotropic plasma effects (Schekochihin et al. 2005). Considering typical parameters of n ∼ 10 −3 cm −3 , T ∼ 10 7 K and B ∼ 1 µG (Ensslin and Vogt 2006), it is possible to show that the cyclotron frequency (Ω) is much larger than the collision frequency (ν ii ). Under such conditions, the plasma fluctuations with wave numbers k ≥1 kpc will be subject to different processes related to the pressure anisotropy. The turbulent cascade, for example, may be modified by the new wave modes and instabilities resulting in a different picture of the energy budget in these environments. It may be particularly important in understanding of the cooling flow and the cosmic rays acceleration processes. Schekochihin et al. (2008) showed that in the ICM the anisotropy-induced firehose and mirror instabilities may grow non-linearly up to the saturation at δB/B ∼ 1. This growth is very fast. They estimated the increase of power at small scales. As a consequence, the excess of small scale fluctuations of B and the energy transport at these environments may drastically change the picture. In certain cases the small scales may also be considered as collisionless, what may be important, for instance, in the development of turbulent cascade. For wavelengths smaller than the Larmor radius the kinetic treatment of plasma is necessary, as shown by Howes et al. (2006). However, in this work we focus on the large scale so the CGL-MHD approximation may be used instead.
The aim of this work is to provide an extensive statistical analysis of the MHD turbulence in collisionless plasmas and to study the role of the different instabilities in the evolution of the system. To accomplish that we performed the first 3dimensional simulations focusing on the evolution of turbulence in the presence of pressure anisotropy. The description of the model as well as the presentation of the governing equations is given in Section 2 with additional discussion of the instabilities and the double-isothermal approximation. In Section 3 we describe the numerical simulations. The results and the extensive statistical analysis including the derivation of probability distribution functions, spectra and structure functions of the fluctuations of density and velocity are presented in Section 4. We discuss the most important results in Section 5 followed by Section 6, where we draw main conclusions.

Governing Equations
In the fluid approximation, a gyrotropic plasma can be described by the Chew et al. (1956) magnetohydrodynamic (CGL-MHD henceforth) equations expressed in the conservative form as follows: ∂ρv ∂t where ρ and v are the plasma density and velocity, respectively, B is the magnetic field, P = p ⊥Î + (p − p ⊥ )bb is the pressure tensor,b = B/|B| is the unit vector along the magnetic field, p and p ⊥ are the pressure components parallel and perpendicular tob, respectively, and f represents the forcing term. The above set of equations is completed by the description of the parallel and perpendicular pressure components. To avoid the complexities related to explicitly calculated processes that may be important for the description of the energies, e.g., the anisotropic heat conduction and the emission and absorption of the radiation, we can make use of the double-polytropic equations instead, as suggested by Chew et al. (1956). Tests of this approximation using the solar magnetosheath data confirmed its validity under the form (see Hau and Sonnerup 1993): where γ ⊥ and γ are the polytropic exponents for the perpendicular and parallel pressures, respectively. These equations, according to Hau (2002), can be expressed in the conservative form, as follows: where S ⊥ = p ⊥ B 1−γ ⊥ , S = p (B/ρ) γ −1 , and B = |B| is the strength of magnetic field. In this form, the double-polytropic equations can be solved numerically similarly as the continuity equation.

Double-Isothermal Closure
Under the double-isothermal closure we have γ = γ ⊥ = 1 and the equation of state is described by two relations, p ⊥ = a 2 ⊥ ρ and p = a 2 ρ, where a ⊥ and a are constants and represent speeds of sound along the perpendicular and parallel directions to the magnetic field, respectively. In such situation the conservative form of the momentum equation for the double-isothermal CGL-MHD model can be rewritten as follows, ∂ρv ∂t where α ≡ 1 2 β − β ⊥ ≡ 1 2 β ⊥ (ξ − 1), ξ ≡ p /p ⊥ is the pressure anisotropy ratio, and β and β ⊥ are the plasma betas in the parallel and perpendicular directions to the local field, respectively.
The main consequence of the double-isothermal closure is that the pressure anisotropy is kept constant, i.e. ξ = const, independently of the evolution of the kinetic instabilities that may arise. Therefore, the stability condition is fulfilled by the local decrease of density and the increase of magnetic pressure and not due to the local changes of the pressure tensor.
Furthermore, this work is focused on studying the differences of the turbulent regimes in the collisionless plasmas and the standard MHD turbulence using as the reference the studies done by Kowal et al. (2007), Kowal and Lazarian (2010), where the authors presented an extensive statistical analysis of density, velocity and magnetic field distributions in the isothermal MHD simulations. For this reason we chose the double-isothermal CGL-MHD model as the natural extension of the isothermal MHD in order to understand the importance of pressure anisotropy and its consequences in the turbulent plasmas. On the horizontal axis we show the plasma beta parameter related to the parallel pressure and defined as β ≡ p mag /p = 1 2 c 2 A /a 2 , where p mag ≡ 1 2 |B| 2 is the magnetic pressure. On the vertical axis we show the squared sinus of the angle between the directions of the perturbation and mean field.

Firehose and Mirror Instabilities
While in the MHD model the dispersion relations do not exhibit any instabilities, in the CGL-MHD equations the term of the pressure anisotropy introduces significant changes to the dynamical system. The detailed linearization and wave analysis of the doubleisothermal CGL-MHD equations is provided in the Appendix. Here we describe the final dispersion formulas resulting from the presence of pressure anisotropy. For example, in the case of the incompressible Alfvén mode, the dispersion relation can be written as where α ≡ 1 2 β − β ⊥ and θ is the angle between the mean field and the wave vector of the perturbation. Equation (5) becomes negative when β − β ⊥ > 2, what results in the occurrence of the firehose instability and the growth of the magnetic field fluctuations. As a consequence, the field lines bend and the magnetic pressure increases reducing the parameter α to its saturation value α ≈ 2. This is known as the kinetic Alfvén firehose instability. Two other instabilities related to the slow mode are called the compressible firehose and mirror instabilities and occur when p > p ⊥ or p ⊥ > p , respectively, and the dispersion relation for magnetosonic waves becomes negative (see Appendix).
In Figure 1 we show the stability regimes for two cases of the pressure anisotropy studied in this paper, a /a ⊥ = 0.5 in the left panel and a /a ⊥ = 2.0 in the right one (corresponding to p /p ⊥ = 0.25 and p /p ⊥ = 4.0, respectively), plotted as functions of the plasma beta β related to the parallel pressure and the squared sinus of the angle between the directions of the perturbation propagation and magnetic field. When the perpendicular pressure dominates (left panel) only the mirror instability corresponding to the slow mode can occur. When perturbations propagate in directions almost perpendicular to the magnetic field the range of unstable β grows. On the contrary, when the direction of propagation is close to the direction of magnetic field the instability does not occur. In the case of dominating parallel pressure (right panel) the plasma becomes unstable due to two firehose instabilities related to the incompressible Alfvén and compressible slow modes, but the Alfvén mode firehose instability extends over larger region of the parameter space up to β = 0.375 independently of the angle between the perturbation propagation and magnetic field directions.
Both instabilities have the growth rate γ ≡ Im(ω) larger at smaller scales, i.e., increasing with the wave number k. This dependence introduces a stiffness in the numerical integration since the micro instabilities, which may be due to the numerical imprecision, tend to grow very fast and destroy the configuration of the studied problem. As a consequence, the pressure anisotropy tends to disappear quickly and the problem of interest cannot be studied. Sharma et al. (2006) studying the magneto-rotational instability (MRI) in protostellar disks pointed out this important problem. They bypassed it by implementing a quasi-stability condition for the computational cells where the pressure anisotropy was larger than the threshold for instability. Physically, it can be understood as an almost "instantaneous" evolution of the system to the quasi-stable condition. Under such condition all anisotropic effects are kept at the large scales and the problem of interest could be analyzed. The downside of this method is that it artificially removes the free energy of pressure anisotropy without, as a counterpart, increasing the kinetic and/or magnetic energies. That is because, in the physical sense, the instability increases the magnetic energy and/or accelerate the gas.
In the case of turbulence this kind of artificial removal of the energy at small scales may not be the best approach since the small scale structures in turbulence are generated by the cascade process operating in the turbulent models and therefore by artificially influencing the dissipation of energy we can distort its physical meaning and obtain incorrect conclusions. Nevertheless, under the double-isothermal approximation the pressure anisotropy is kept even after the growth of the small scale perturbations saturates, so the stability condition at large scales may lay in the unstable region even after the saturation at small scales. Therefore, under certain conditions, the doubleisothermal approximation is valid and might be an interesting area of studies of the evolution of turbulence in a collisionless plasma.

Numerical Simulations
The simulations of turbulence in collisionless plasma were performed solving the set of double-isothermal CGL-MHD equations in a conservative form given by equations (1a)-(1c) with an addition term f in the motion equation representing the turbulence driving. The numerical integration of the system evolution governed by the CGL-MHD equations were performed using the second order shock-capturing Godunov-scheme code (Kowal et al. 2007, 2009, Kowal and Lazarian 2010. We incorporated the field interpolated Table 1. Description of the performed simulations of the double-isothermal CGL-MHD turbulence with the resolution 512 3 . Initial density for all models was set to 1.0. constrained transport (CT) scheme (see, e.g., Tóth, 2000) into the integration of the induction equation to maintain the ∇·B = 0 constraint numerically, the general Harten-Lax-van Leer Riemann solver (Einfeldt et al. 1991) to obtain the numerical fluxes. The time integration was done with the second order Runge-Kutta method (see e.g. Press et al. 1992). On the right-hand side, the source term f represents a random solenoidal large-scale driving force. The rms velocity δv is maintained to be approximately unity, so that v can be viewed as the velocity measured in units of the rms velocity  (Kowal et al. 2007, Kowal andLazarian 2010), the sound speeds and the strength of the external field B 0 are the controlling parameters defining the sonic Mach number M s = δv/a and the Alfvénic Mach number M A = δv/c A , respectively. The angle brackets signify the averaging over the volume. M s < 1 and M s > 1 define subsonic and supersonic regimes, respectively, and M A < 1 and M A > 1 define another two regimes, sub-Alfvénic and super-Alfvénic, respectively. Since these two parameters are independent we can analyze, e.g., supersonic sub-Alfvénic turbulence, which signifies that M s > 1 and M A < 1. In the case of the CGL-MHD simulations the regimes are defined by two sonic Mach numbers corresponding to the parallel and perpendicular sounds speeds.
We drove turbulence solenoidally at wave scale k equal to about 2.5 (2.5 times smaller than the size of the box). This scale defines the injection scale in our models. We did not set the viscosity and diffusion explicitly in our models. The scale at which the dissipation starts to act is defined by the numerical diffusivity of the scheme.
We performed six three-dimensional CGL-MHD simulations using the resolution 512 3 for different initial conditions, as shown in Table 1. We simulated the clouds up to t max ∼ 6, i.e. 6 times longer than the dynamical timescale, to ensure a full development of the turbulent cascade. The computational time required for each CGL-MHD simulation with intermediate resolution was equivalent to a MHD 512 3 , as performed by Kowal et al. (2007) and Kowal and Lazarian (2010).
Models 1 and 2 are examples of weak turbulence and belong to subsonic and subAlfvénic regimes. The difference between these models is the pressure anisotropy accounting for the mirror and firehose instabilities. the comparison of both models gives us an insight into the different evolution of turbulence in these two cases. Similarly we calculated Models 3 and 4 for strong turbulence resulting in the supersonic and superAlfvénic regimes. Again, we tested these regimes for different instabilities. Model 5 belongs to the subsonic and superAlfvénic regime representing the physical conditions of the ICM, and can be used as reference for further studies on this subject. Finally, Model 6 belongs to the supersonic and subAlfvénic turbulent regime and is the only simulation initiated with magnetic pressure larger than the thermal one. In this model all cells are initially stable but as the turbulence develops unstable conditions can occur. In the following sections we present the results obtained for each model, as well as a direct comparison with the standard MHD turbulent simulations with similar initial conditions.

Distribution of Density and Velocity
The results obtained for the distribution of density and velocity show strong differences between the CGL-MHD and standard MHD models. However, the kinetic instabilities play a role in the evolution of turbulence under specific conditions. In Figure 2 we present the column density obtained for each CGL-MHD model, as as well as the MHD models for comparison. Each case is presented as labeled in Table 1, with the MHD case shown in the left, and the CGL-MHD in the right.
From the column density maps it is possible to distinguish models 2 and 5. Both models are subsonic, where initially we set p /p ⊥ = 2, which resulted in strong firehose instabilities. Here, the firehose instability is responsible for a deformation of the magnetic field lines. The curved magnetic lines tend to slow down and trap the flowing gas. Since the growth rate is larger at small scales we expect this effect to create more granulated maps. This is not seen in model 4 because ρδv 2 > p δB 2 and, therefore, the turbulence is able to destroy the configuration of the growing instability. The same occurs for model 3. Even though it is not visible in the column density maps, the kinetic instabilities are responsible for changes in the statistics of the turbulence, as addressed below. For models 1 and 6, in which p /p ⊥ = 0.5, the column density maps show mild differences. In these cases we have a mirror instability operating, which is responsible for changes in the velocity distribution. For model 1 with a weak turbulence, the mirror instability is responsible for the acceleration of the gas resulting in the increase of the effective sonic Mach number. In model 6, where initially all cells were stable, the evolution of turbulence causes the instability in most of the computational domain. Column density for the different MHD and CGL-MHD models corresponding to the models presented in Table 1.
The instability is responsible here for slowing the gas and reducing the effective sonic Mach number.
In Figure 3 we present the probability distribution functions (PDFs) of density obtained from each model. In general, the distribution of gas exhibits an increasing contrast with increasing sonic Mach number. This result is in agreement with the Figure 3. PDFs of density for the studied models. The left column shows models 1 and 2, the middle one shows models 3 and 4, and the right column shows models 4 and 5, according to Table 1. For each case both, the CGL-MHD (solid lines) and MHD (dashed lines) models are shown for comparison. Figure 4. PDFs of velocity for the studied models. The left column shows models 1 and 2, the middle one shows models 3 and 4, and the right column shows models 4 and 5, according to Table 1. For each case both, the CGL-MHD (solid lines) and MHD (dashed lines) models are shown for comparison.
MHD models. However, the kinetic instabilities cause even larger density contrast in the models with weak turbulence, compared to the MHD models. As stated above, in models 1 and 2 the instabilities could freely grow without being suppressed by the turbulent motions of the gas. The PDFs of density for the remaining models show no difference when compared to the MHD models. Surprisingly, in model 5 the firehose instability is responsible for the granulated map of column density, but the PDF of density remains very similar to the MHD case.
The PDFs of velocities are shown in Figure 4. Again, the supersonic and super-Alfvenic models (3 and 4) show small differences compared to the standard MHD case. On the other hand, the subsonic and sub-Alfvenic models (1 and 2) present distorted velocity PDFs, with an increase in the high velocity tail. Both instabilities are responsible for more effective acceleration of the plasma, broadening the distribution of velocities. In the case of the firehose instability, the weakly magnetized model presents more evident changes, as noted in model 5. Here, since the magnetic field is weak the perpendicular flows are able to destroy the magnetic field configuration. Therefore, the magnetic breaking is not important, what results in flows with higher velocities. The same process is responsible for the changes in the PDF of model 6. Interestingly, for model 6 we obtained narrower PDF with lower velocities when compared to the MHD case. As the turbulence develops and unstable cells arise the local magnetic field grows. The strong magnetic breaking takes place resulting in a low velocity distribution.

Power Spectra of Density and Velocity
In order to characterize and study the changes in the energy cascade, as well as the correlations between different scales in CGL-MHD models, we analyze the power spectra. In Figures 5 and 6 we present the spectra of fluctuations of density and velocity for all models. Within each plot, we also present the anisotropy of the spectra for parallel and perpendicular directions with respect to the global magnetic field.
At the center, we show the plots for Models 3 and 4. As also shown previously, there is no substantial difference between the CGL-MHD and the MHD simulations. The density spectra present similar slope (α ∼ −0.5) at the inertial range, while a slope of ∼ −2 is obtained for the velocity spectra. Regarding the anisotropy (l ∝ l ζ ⊥ ), we found the typical relations ζ ∼ 1 at small scales, but the Goldreich and Sidhar (1995) slope ζ ∼ 2/3 at larger scales, for both MHD and CGL-MHD simulations. A similar result in spectra and anisotropy of density was obtained for Model 6, with exception of a slight increase in the spectrum at small scales.
Clearly, differences appear between the spectra of the MHD and CGL-MHD simulations in Models 1, 2 and 5. In these cases, it is noticeable the power excess at large values of k due to the fast growth of the instabilities in these scales. The slopes of velocity and density spectra for the MHD models range between -1.7 to -2.0 at the inertial range, while CGL-MHD simulations show positive slopes (∼ +1 for Models 1 and 5) in velocity spectra, and flat density spectra (α ∼ 0). Regarding the spectral anisotropy, as an interesting result, we see that the firehose instability reduces the anisotropies regarding the global magnetic field lines. A more detailed study of the anisotropy of density and velocity perturbations regarding the local magnetic field is given below. Figure 5. Spectra of density for the studied models, models 1 to 3 in the upper row and models 4 to 6 in the lower row, according to Table 1. For each case both, the MHD (solid lines) and CGL-MHD (dashed lines) models are shown for comparison.

Structure functions
In the previous section we showed, from power spectra of density and velocity fluctuation, that the instabilities in CGL-MHD models efficiently transport power from large to small scales, where their concentration is increased. Interesting modifications from MHD turbulence was also shown from the spectral anisotropy regarding the global magnetic field. However, since the instabilities are dominant at small scales, mapping the structure functions regarding the local magnetic field seems to be a better approach Figure 6. Spectra of velocity for the studied models, models 1 to 3 in the upper row and models 4 to 6 in the lower row, according to Table 1. For each case both, the MHD (solid lines) and CGL-MHD (dashed lines) models are shown for comparison.
if we want to determine the role of the instabilities on the isotropization of fluctuations.
The second order structure function of a given parameter f is defined as: where r is the referenced position and l is the distance calculated along the magnetic field line. The SF is calculated by randomly choosing a large number of referenced positions for each studied correlation length l. Here, to account for the importance of magnetic field fluctuation at small scales, we calculate the correlation length l along the field lines, i.e. in the magnetic field reference frame. Figure 7. Structure functions of the density in the local reference frame for the studied models. In the left column we show models 1 and 2, inte middle column models 2 and 3, and in the right column, models 5 and 6, according to Table 1. For each case both, the CGL-MHD (solid lines) and MHD (dashed lines) models are shown for comparison.
The obtained SFs for density and velocity fluctuations are shown in Figures 7 and  8, respectively. In all plots, the solid lines represent the results obtained from the MHD simulations, while the dashed lines were obtained from the CGL-MHD simulations. The numbers refer to each model, as described in Table 1.
Similar to the previous calculations, Model 3 shows no difference between the two theoretical approaches, independent on the scale. The same result is obtained for Model 6, which also showed similar spectra for the CGL-MHD and MHD models. Surprisingly, the structure functions of Model 4, which showed no differences in PDF and spectra, present more isotropic maps in the CGL-MHD model. The same behavior is obtained for Models 2 and 5. These are the models presenting p /p ⊥ = 2, i.e. the firehose instability. Here, the firehose instability is responsible for changes in the magnetic field topology, tangling the field lines, resulting in an increase of the perpendicular pressure in the local reference frame. Obviously, this effect is over estimated in these calculations because of the double-isothermal approximation. Otherwise, the magnetic field topology would not change drastically, but the interchange of energy would cause a reduction of parallel pressure anyway. In this sense, the result would be the same, i.e. the firehose instability is responsible for a isotropization of the fluctuations with respect to the Figure 8. Structure functions of velocity in the local reference frame for the studied models. In the left column we show models 1 and 2, inte middle column models 2 and 3, and in the right column, models 5 and 6, according to Table 1. For each case both, the CGL-MHD (solid lines) and MHD (dashed lines) models are shown for comparison. magnetic field lines. Model 1, on the other hand, presented a larger anisotropy for the CGL-MHD model. Here, as p /p ⊥ = 0.5, the free energy is mostly converted to kinetic pressure. The mirror instability is responsible for an increase in the acceleration of the plasma along the field lines (as already seen in the PDF of velocity), increasing the anisotropy of the fluctuations. The result is more elongated structures, mostly at small and intermediate scales, as also noticed in the power spectra.

Time Evolution of Unstable Systems
In order to understand the evolution of the system during the growth and saturation of the instabilities, we calculated the instability condition for the Alfvén and the compressible modes for each cell of the computational domain. The dispersion relation of these modes reveal the cells where the instability condition is fulfilled.
As a general behavior, we found that the turbulence increases the range of unstable cells, even for the case with initially stable cells (Model 6). This process is illustrated in Figure 9 where we show the time evolution of the dispersion relation of Model 4 (similar distributions were obtained for all models). Here, we plot the histograms of the stability condition for different times. The Alfvén and magnetosonic modes are independently shown in the left and right plots, respectively. Unstable cells populate the negative range of the dispersion relation. The time is shown in units of the dynamical timescale τ d ∼ δV /L.
As the turbulence is injected in the simulation domain, fluctuations of density and magnetic field change the local characteristic speeds and, as a consequence, the stability conditions. This process is fast compared to the dynamical timescale of the system (t < 0.5τ d ). As a result, the dispersion relation spread over a larger range including the negative values of (ω/k) 2 . Then, the instabilities start to grow and saturate, at different timescales for different scales. The saturation of the instability brings the cells towards positive values of (ω/k) 2 . At t ∼ 2.5, most of cells have already reached the saturation condition for both modes. We believe that the few still unstable conditions are related to the large scale fluctuations, which evolve slow and has higher saturation values. It also seems that the firehose instability of Alfvén modes saturates faster than the magnetosonic mode, though the differences between the modes could be also be related to the driving mechanism. We plan to address this possibility in a future work.
In the present calculations we did not perform any wave mode decomposition and, therefore we were unable to analyze the time evolution of the individual modes. Also, for this reason we are unable to characterize the time evolution for each wavelength (scale). Needless to say that this is an interesting subject for future studies on CGL-MHD turbulence.

Conclusions
In this work we presented the first extensive statistical analysis of three-dimensional simulations of turbulence in collisionless plasma. We studied the case of double-isothermal closure of the CGL-MHD equations in order to compare our results with previous isothermal simulations of MHD turbulence. We performed simulations with different characteristic sonic and Alfvenic Mach numbers, as well as different pressure anisotropies to account for both firehose and mirror instabilities. As main results we showed that: • we obtained firehose/mirror unstable conditions in all simulations. The unstable conditions may be created, even for stable initial conditions, as a consequence of the driving and evolutions of turbulent cascade; • the supersonic and super-Alfvenic models showed no significant differences when compared to standard MHD models. Basically, strong turbulence is able to destroy the changes in the topology resulted from the instabilities; • the PDF's of density showed broadened profiles for the subsonic and sub-Alfvenic cases. The PDF's of velocity showed changes for sub-Alfvenic models. Specifically, we obtained an increase on the number high velocity flows in subsonic models, while its decrease for supersonic turbulence; • the spectra of density and velocity showed an increase of power in small scales for subsonic models; • the structure functions of velocity and density fluctuations revealed that the firehose instability tends to isotropize the fluctuations regarding the local reference frame, i.e. along the magnetic field lines. On the other hand, the mirror instability increases the elongation of the fluctuations along the magnetic field lines; • the dynamical timescale (τ d ∼ δV /L), may also be a good estimate for the saturation timescale of the instability growth of most of the wavelengths. This fast evolution of the system implies in interesting physical processes in, e.g. interchange of energy and acceleration of cosmic rays in the collisionless plasma at intracluster medium of galaxies. The growth rate of long wavelengths may be much larger than τ d .
The ideal double-isothermal CGL-MHD equations with the pressure isotropy assumption p ≡ a 2 ρ = p ⊥ ≡ a 2 ⊥ ρ results in P = a 2 ⊥ ρI and can be written in the non-conservative form as follows ∂ρ ∂t where ρ is the density, v is the velocity field, B is the magnetic field, and a ⊥ and a are the sound speeds in the perpendicular and parallel directions with respect to B, respectively. We assume that all variables can be separated into the uniform and fluctuating components, i.e., ρ → ρ 0 + δρ, v → v 0 + δv, B → B 0 + δB. We also assume the lack of uniform flow, i.e., v 0 = 0. Substituting the variable separation in equations (A.1)-(A.3) and removing all non-linear terms leads to the following set of equations ∂δρ ∂t − k (k · B 0 ) (B 0 · δv) + (k · B 0 ) (B 0 · k) δv} = 0, (A.10) Without loosing the generality we can assume that the mean magnetic field is parallel to the X direction, i.e. B 0 = B 0x , and k lays in the XY-plane under an angle θ with the respect to B 0 , i.e. k = k (cos θx + sin θŷ). We also introduce the Alfvén speed c A ≡ B 0 / √ 4πρ 0 . Substituting these assumptions and diving the dispersion relation by k 2 we can express it in a matrix form − ω 2 k 2 + a 2 ⊥ cos 2 θ a 2 ⊥ sin θ cos θ 0 a 2 ⊥ sin θ cos θ − ω 2 k 2 + a 2 ⊥ sin 2 θ + c 2 Determinant of the matrix A gives the dispersion relation det A = − ω 2 k 2 + c 2 A cos 2 θ ω 4 k 4 − ω 2 k 2 a 2 ⊥ + c 2 A + c 2 A a 2 ⊥ cos 2 θ = 0, (A.12) with the eigenvalues (A.13) corresponding to the Alfvén wave and corresponding the fast and slow magnetosonic waves, respectively.
In the next step we include the pressure anisotropy term in the motion equation (A.2) which for the double-isothermal approximation can be written as Substituting the variables separation in this term we obtain whereb 0 ≡ B 0 /B 0 , and the first term corresponds to the perturbation of the density, and two remaining terms correspond to the perturbation of magnetic field. Introducing the perturbed variables of the form as before, i.e. δf (x, t) ∝ exp [i (k · x − ωt)], the term (A.16) can be written in Fourier space as follows, Next, we multiply it by ω and substitute equations (A.4) and (A.6) to obtain simpler relation ωk · δ p − p ⊥ bb = ρ 0 a 2 − a 2 ⊥ k ·b 0 2 2b 0 b 0 · δv − δv . (A.18) Finally, assuming B 0 = B 0x and k = k (cos θx + sin θŷ) the term reduces to a very simple expression ω k 2 1 ρ 0 k · δ p − p ⊥ bb = 2δv xb0 − δv a 2 − a 2 ⊥ cos 2 θ.