PIC Simulations of VelocitySpace Instabilities in a Decreasing Magnetic Field: Viscosity and Thermal Conduction
Abstract
We use particleincell (PIC) simulations of a collisionless, electronion plasma with a decreasing background magnetic field, B, to study the effect of velocityspace instabilities on the viscous heating and thermal conduction of the plasma. This work is complementary to Riquelme et al. (2016), where the case of an amplifying magnetic field was studied. If decreases, the adiabatic invariance of the magnetic moment gives rise to pressure anisotropies with ( and represent the pressure of species [electron or ion] parallel and perpendicular to B). Linear theory indicates that, for sufficiently large anisotropies, different velocityspace instabilities can be triggered. These instabilities, which grow on scales comparable to the electron and ion Larmor radii, in principle have the ability to pitchangle scatter the particles, limiting the growth of the anisotropies. Our PIC simulations focus on the nonlinear, saturated regime of the instabilities. This is done through the permanent decrease of by an imposed shear in the plasma. Our results show that, in the regime (), the saturated ion and electron pressure anisotropies are controlled by the combined effect of the oblique ion firehose (OIF) and the fast magnetosonic/whistler (FM/W) instabilities. These instabilities grow preferentially on the ion Larmor radius scale, and make the ion and electron pressure anisotropies nearly equal: (where ). We also quantify the thermal conduction of the plasma by directly calculating the mean free path of electrons, , along the mean magnetic field. A simple physical prescription for is obtained, which strongly depends on whether decreases or increases. Our results can be applied in studies of low collisionality plasmas such as the solar wind, the intracluster medium, and some accretion disks around black holes.
Subject headings:
plasmas – instabilities – accretion disks – solar wind1. Introduction
In lowcollisionality plasmas, the change in the magnitude of the local magnetic field () generically drives a pressure anisotropy with (where and correspond to the pressure of species perpendicular and parallel to B). This is a consequence of the adiabatic invariance of the magnetic moment of particles, , in the absence of collisions (where is the velocity of species perpendicular to B).
These pressure anisotropies can trigger various velocityspace instabilities, which are in principle expected to pitchangle scatter the particles, to some extent mimicking the effect of collisions. The combined effect of pressure anisotropies and velocityspace instabilities can affect various large scale properties of the plasma, including its effective viscosity (Sharma et al., 2006; Squire et al., 2017) and thermal conductivity (see, e.g., Riquelme et al., 2016; Komarov et al., 2016). This weakly collisional behavior is expected to be important in several astrophysical systems, including lowluminosity accretion flows around compact objects (Sharma et al., 2007), the intracluster medium (ICM) (Schekochihin et al., 2005; Lyutikov, 2007), and the heliosphere (Maruca et al., 2011; Remya et al., 2013).
In a previous work (Riquelme et al., 2016) we studied how the plasma viscosity and thermal conductivity are affected by an increase of , which naturally drives . In this paper we study the opposite case, where decreases and . In this case several velocityspace instabilities can be excited, ultimately regulating the extent to which the anisotropy can grow. When only the electron dynamics is considered, two types of plasma waves are expected to be driven unstable by the electron pressure anisotropy: the oblique electron firehose (OEF) modes, which are purely growing modes, and the Alfvén/ioncyclotron (A/IC) modes, which are quasiparallel, propagating waves, driven unstable by cyclotronresonant electrons (Li et al., 2000).^{1}^{1}1Although the Alfvén/ioncyclotron (A/IC) modes grow at wavelengths comparable to the electron Larmor radius, their name indicates that they correspond to the Alfvén branch that starts as low wavenumber () Alfvén mode and then becomes the ioncyclotron mode at higher . Similarly, in presence of an ion pressure aniostropy , there are also two types of modes that can grow unstable: the oblique ion firehose (OIF) modes, which are purely growing modes, and the fastmagnetosonic/whistler (FM/W) modes, which are quasiparallel, propagating waves, excited by cyclotronresonant ions (Quest et al., 1996; Gary et al., 1998; Hellinger et al., 2000).^{2}^{2}2Although the OIF and OEF modes are technically the same mode (they correspond to the same Alfvénmode branch), we will consider them as separate instabilities since they are identified as two different growthrate maxima in kspace.
Although the linear properties of these instabilities have been relatively well studied, their long term, nonlinear behavior is less clear. This is interesting since in realistic astrophysical settings the pressure anisotropies are expected to be driven (via decrease) for a time significantly longer than the initial regime where the instabilities grow exponentially. Also, since the decreasing magnetic field should simultaneously drive electron and ion anisotropies, we would like to understand the combined effect of these two species on the different unstable modes. In this paper, we achieve this goal by continuously decreasing the strength of the backgroud magnetic field in our simulations through externally imposing a shear motion in the plasma. This setup allows us to drive and maintain the relevant instabilities in their nonlinear, saturated regime during most of the simulation time.
There are two important applications of our work. One is to quantify the so called “anisotropic viscosity” of the plasma, which is controlled by the pressure anisotropies of the particles. This viscosity is believed to contribute significantly to the heating of electrons and ions in accretion disks and other lowcollisionality plasmas (Sharma et al., 2006, 2007; Squire et al., 2017). Also, the nonlinear evolution of the different velocityspace instabilities sets the pitchangle scattering rate of electrons, which is key to determine their mean free path and, therefore, the thermal conductivity of the plasma.
The paper is organized as follows. In §2 we describe our simulation setup and strategy. In §3 we determine the saturated pressure anisotropy of ions and electrons, describing the responsible physical mechanisms. In §4 we quantify the ion and electron heating. In §5 we measure the mean free path of electrons and ions, and determine their dependence on the physical parameters of the plasma. In §6 we summarize our results and discuss their implications for various lowcollisionality astrophysical plasmas.
2. Simulation Setup
We use the electromagnetic, relativistic PIC code TRISTANMP (Buneman, 1993; Spitkovsky, 2005) in two dimensions. The simulation box consists of a square box in the  plane, containing an initially isotropic plasma with a homogeneous initial magnetic field . We simulate a decreasing magnetic field by imposing a velocity shear given by , where is the shear rate of the plasma and is the distance along ( and are the unit vectors parallel to the and axes, respectively). From flux conservation, the and components of the mean field evolve as and . Thus, if and are positive, there will be a decrease of and, therefore, of . Therefore, we initially choose , which garantees a decrease of and a anisotropy during a simulation time .
By resolving the  plane, our simulations can capture the quasiparallel A/IC and FM/W modes, as well as the oblique OEF and OIF modes with their wave vectors k forming any angle with the mean magnetic field . The key parameters in our simulations are: the particle magnetization, quantified by the ratio between the initial cyclotron frequency of each species and the shear rate of the plasma, (), and the ion to electron mass ratio, . In typical astrophysical cases, and . Due to computational constraints, however, we will use values of and much larger than unity, but still much smaller than expected in real environments. This limitation will be taken into account when applying our simulation results to relevant astrophysical cases.
Our simulations initially have (). In almost all of our runs , which implies (where , , and are Boltzmann’s constant, the electron temperature, and the electron plasma frequency, respectively). We will change our simulation conditions by varying: and (which uniquely fix and ). Some of our simulations use “infinite mass ions” (the ions are technically immobile, so they just provide a neutralizing charge), with the goal of focusing on the electronscale physics. These provide a useful contrast with our finite runs and allow us to isolate the impact of ion physics on the electrons. The numerical parameters in our simulations will be: N (number of particles per cell), (the electron skin depth in terms of grid size), (box size in terms of the initial ion Larmor radius for runs with finite ; , where ), and (box size in terms of the initial electron Larmor radius for runs with infinite ). Table 1 shows a summary of our key runs. We ran a series of simulations ensuring that the numerical parameters (e.g., different N) do not significantly affect our results. Note that most runs used just for numerical convergence are not in Table 1.
Runs  N  L/R  

I1  3600  0.28  40  210  
I2  7200  0.28  40  210  
I3  7200  0.1  40  210  
F1  64  7200  0.28  40  640 
F2  25  7200  0.28  40  400 
F3  25  3600  0.28  40  400 
F4  10  7200  0.28  40  250 
Note. – A summary of the physical and numerical parameters of our simulations. These are the mass ratio , the initial electron magnetization , the ratio between electron temperature and rest mass energy, , the number of particles per cell N (including ions and electrons), and the box size in units of the typical initial electron Larmor radius (, where , is Boltzmann’s constant, and is the electron temperature). All of the runs initially have and , where is the simulation time step. The runs shown in this Table share the same electron skin depth (where is the grid point separation), but we used several other simulations to confirm numerical convergence by varying , N, and L/R.
3. Pressure Anisotropies
In this section we focus on the nonlinear evolution of the ion and electron pressure anisotropies. As stated above, we will begin by showing simulations where ions have infinite mass.
3.1. Simulations with
Figure 1 shows the magnetic field fluctuations and plasma density in a simulation in which the ions have infinite mass (run I1 in Table 1). Thus the electrons can only be affected by the electron anisotropydriven OEF and A/IC instabilities. The upper row in Figure 1 corresponds to , i.e., after one shear time, while the lower row corresponds to . Figure 1 shows that at all times the magnetic fluctuations are dominated by their component, with wavenumbers (, where k is the wave vector) satisfying , and with k being mainly oblique to the mean direction of B. This suggests that the OEF modes contribute the most to the amplitude of the magnetic fluctuations. However, although smaller in amplitude, the presence of quasiparallel A/IC modes can also be seen especially in the component.
Figure 2 shows the time evolution of the magnetic energy in , the volumeaveraged pressure anisotropy, and the electron magnetic moment for two runs, one with and the other with (runs I1 and I2 in Table 1, respectively).
Panels and show the volumeaveraged pressure anisotropy for these two runs, where . For comparison, in both cases we plot the anisotropy threshold that would make the OEF and A/IC modes grow at a rate equal to the shearing rate, . These thresholds were calculated using the linear Vlasov solver developed by Verscharen et al. (2013). The calculations use mass ratio and assume very cold ions (), which seeks to resemble our simulated situation where the ions only provide a neutralizing charge.^{3}^{3}3In order to make sure that using mildly relativistic electrons in our runs (where ) does not invalidate our comparison with the calculated thresholds (which assume nonrelativistic electrons), in Figure 2 we added in solidred line the pressure anisotropy for run I3, which uses (while keeping the same and ). The very small difference between the two cases suggests that the effect of having mildly relativistic electrons is fairly small. We see that the OEF and A/IC thresholds are quite similar, although with the OEF mode having a slightly lower threshold, especially in the case . This implies that both modes should play some role in regulating the electron anisotropy, with their relative importance depending weakly on the ratio . Also, for both values of there is a reasonably good agreement between the electron anisotropy obtained from the simulations and the linear OFE and A/IC instability thresholds. This is thus consistent with the electron anisotropy being maintained at the level for the OEF and A/IC modes to grow at a rate close to .
The contribution of the different components of can be seen from Figures 2 and 2, which show the magnetic energy along different axes, normalized by the average magnetic energy in the simulation, . is decomposed in terms of (component perpendicular to the simulation plane), (component parallel to the simulation plane but perpendicular to ), and (component parallel to ). Clearly, is dominated by its component (as already seen in Figure 1). This shows that, although the OEF and A/IC modes are expected to contribute to limiting the electron anisotropy, their contribution to the magnetic energy in is quite different. Indeed, our linear calculations show that, for the plasma parameters of runs I1 and I2, the OEF modes should satisfy . Thus the fact that in our runs implies that most of the magnetic energy is being contributed by the OEF modes. The quasiparallel A/IC modes, which are most visible in the component as can be seen from Figure 1, make a significantly smaller contribution to . These different contributions to are likely due to the slightly different anisotropy thresholds of the OEF and A/IC instabilities, as well as to the different amplitude at saturation expected for these two modes. Another possible factor is that the A/IC instability growth rate is very sensitive to the orientation of the pitchangle gradients of the distribution function. Therefore, the A/IC instability can relax the distribution faster toward a stable configuration through pitchangle scattering than the OEF instability.
Finally, Figures 2 and 2 show the volumeaverage magnetic moment of the electrons, (; blacksolid line), for the same runs I1 and I2, respectively. It can be seen that until the onset of the OEF and A/IC exponential growth (), is fairly constant, implying the lack of efficient pitchangle scattering. After that, tends to increase at a rate close to the shear rate . This implies the appearance of an effective pitchangle scattering rate for the electrons, , due to their interaction with the OEF and A/IC modes.
In order to help us to understand the way is regulated by the different velocityspace instabilities, we propose a second way to calculate the average magnetic moment of species :
(1) 
This definition is useful because there can be cases where . This is expected when, besides pitchangle scattering, is partly regulated by relatively large fluctuations in , which may spatially correlate with in a conserving way. This occurs, for instance, in the presence of large amplitude mirror modes (Kunz et al., 2014; Riquelme et al., 2015, 2016). Figures 2 and 2 show that, after conservation is broken, (dottedblack) and (solidblack) are almost indistinguishable. This confirms that is regulated by an effective pitchangle scattering provided by the OEF and A/IC modes, with no significant contribution from fluctuations in .
3.2. Simulations with finite
In order to study the effect of ions in regulating both the ion and electron pressure anisotropies, we now focus on simulations with finite ion to electron mass ratios . Since using is computationally infeasible with our current resources, we have instead tried to ensure that our simulation results do not depend significantly on , which is reasonably well achieved for and 64.
As an example, Figure 3 shows the three components of for run F1 of Table 1 ( and ). The upper and lower rows correspond to and , respectively. At both times, quasiparallel and oblique modes are present with similar amplitudes, and with wavenumbers satisfying (where is the ion Larmor radius). While the quasiparallel modes are apparent in the three components of , the oblique modes appear mainly in . These features are consistent with the simultaneous presence of both OIF and FM/W modes. The presence or absence of fluctuations on electron Larmor radius scale is less clear from the figure. We will come back to this question below.
3.2.1 The Role of
Figure 4 compares the evolution of the energy in , the ion and electron anisotropies, and and for simulations with different mass ratios and electron magnetization. The first and second columns compare simulations with the same electron conditions (, , and ) but with different mass ratios: (run F2) and (run F1), respectively.
Figures 4 and 4 show the volumeaveraged electron and ion pressure anisotropies as a function of time (green and black lines, respectively) for runs F2 and F1, respectively. These figures also show the anisotropy threshold for the growth of different instabilities, using a growth rate of . In dotted lines we show thresholds that assume , and correspond to the instabilities: OEF (dottedgreen), OIF (dottedred), and FM/W (dottedblue). We do not include A/IC thresholds in this case. This is because our linear calculations show that the A/IC modes are subject to cyclotronresonant ion damping when , becoming stable in our runs with finite (this is true even for ).^{4}^{4}4The A/IC modes are relevant for the simulations with fixed ions because these particles behave like very cold ions and can not damp the A/IC modes. We see that our simulation results are fairly independent of the mass ratio, and can be summarized as follows:

The ion and electron anisotropies evolve quite similarly in both cases. (This justifies using instability thresholds that assume for comparison.)

The obtained electron anisotropy is a factor smaller than the expected OEF threshold.

The ion and electron anisotropies are best described by the thresholds of the OIF and FM/W instabilities with .

The OIF and FM/W thresholds with are quite similar, which is consistent with the simultaneous presence of these modes in Figure 3.
The fact that the electron anisotropy is close to the OIF and FM/W thresholds, and a factor smaller than the expected OEF threshold, shows that the OIF and FM/W modes are the ones with the largest effect on the electron anisotropy. This can be understood as due to the significant contribution of the electron pressure anisotropy to the growth of OIF and FM/W modes. Indeed, our linear calculations show that the thresholds for the OIF and FM/W instabilities with are a factor of larger than in the case. This conclusion is also supported by the fact that, for the obtained electron anisotropy, the OEF modes are stable, indicating that the contribution from the OEF modes to the scattering of electrons is not expected to be important.
Finally, we performed similar linear threshold calculations for the case with (where represents the growth rate of the different instabilities), which are shown in Figure 6. We find that the thresholds of the OIF (dottedred) and FM/W (dottedblue) modes with continue to be similar and smaller than the OEF (dottedgreen) threshold by a factor (the dottedgreen line almost coincides with the solidred line). Also, the OIF and FM/W thresholds with are times smaller than in the case. This implies that, for realistic mass ratios and magnetizations, the dominant instability for the regulation of ion and electron anisotropies should continue to be the OIF and FM/W instabilities.^{5}^{5}5Although the OIF and FM/W thresholds are similar, there is the trend for the FM/W threshold to be smaller than the OIF threshold at early times (, ), while the opposite situation happens at late times (, ). This implies that in the more realistic cases there could be a clearer dominance of the FM/W (OIF) modes for (). This is consistent with the OIF dominance shown in the hybridPIC simulations with fluid electrons presented by Kunz et al. (2014), which use and where the ions are significantly more magnetized than in our simulations.
Figures 4 and 4, show the magnitude of the volumeaveraged magnetic energy in the same three components of plotted in Figure 2: (solidgreen), (solidblack), and (solidred), for and 64, respectively (for comparison the case is replicated in Figure 4). For the two mass ratios, the two components perpendicular to dominate, with being most of the time times larger than . This result implies that the OIF and FM/W modes are contributing comparable energy to , as was also noticeable from Figure 3. Indeed, our linear calculations show that for the OIF modes , while for the FM/W modes ,^{6}^{6}6Indeed, using the Vlasov solver of Verscharen et al. (2013) one can obtain that, for the parameters of runs F1 () and F2 (), the OIF modes satisfy . which implies that most of the component is being produced by FM/W modes.
The amplitudes of the OIF and FM/W modes appear to depend on the mass ratio. Although time dependent, the magnitude of in the case is on average times larger than in the case. Since the magnetizations of the two runs differ by a factor (), this is roughly consistent with previous studies of the OIF instability that show that at saturation should scale as (Kunz et al., 2014). On the other hand, in the case is about times larger than for , which is roughly consistent with the expectation for the FM/W modes to have a saturation amplitude that satisfies . This scaling can be obtained from the expected effective ion scattering frequency by resonant waves, , which scales as
(2) 
where and are the wave vector component and the particle velocity component parallel to (Marsch, 2006). For the case of the quasiparallel FM/W waves, we obtained from the simulations that , meaning that for most particles and that . Thus, since at FM/W saturation one expects (see Equation 4 below), this implies that (considering that the change in between the and 64 cases is small).
Finally, in panels 4g and 4h we compare and for both ions and electrons. We see that the change in tends to be somewhat larger than the one in (by ) for the two mass ratios tested. This implies that the combined effect of the OIF and FM/W modes reduces in a way that mainly breaks the adiabatic invariance of , with the preservation of due to changes in the field configuration playing a small role.
3.2.2 The Role of
It is also important to understand the role of electron magnetization, , while keeping the same mass ratio. This can be done by looking at the first and third columns in Figure 4, which compares simulations with but with electron magnetization and 3600 (runs F2 and F3 in Table 1, respectively). We see that the two runs reproduce essentially the same results in terms of and evolution, with the only difference being in the amplitude of the magnetic fluctuations. Apart from some significant time variability, the amplitude of the and components in the , run are quite similar to the case , . Since these two runs have a very similar ratio ( and , respectively), this result is consistent with the dependence of the OIF and FM/W saturated amplitude on mentioned in §3.2.1.
3.2.3 Breaking of Adiabatic Invariance
An important question is whether the ionscale instabilities alone are capable of explaining the break in the electron magnetic moment shown in Figure 4, 4, and 4 (which starts at ).
We explore this issue by comparing the power spectra of the fluctuations in our finite runs with the power spectrum produced in the case with . This is done in Figures 5 and 5, where the magnetic energy per logarithmic unit of and is plotted at for runs with and (runs F4, F2, F1, and I2, respectively; and are the magnitude of the wave vector components parallel and perpendicular to , respectively). The electrons in these simulations have the same conditions (, , and initial ), so the different results are only due to the different . We see that:

In the cases with finite mass ratio, as increases the peaks of the spectra shift to longer wavelengths (in units of ) in a way consistent with the growth of the ratio.

In the same way, as increases there is a growth in the amplitude of the peak of the spectra, which accounts for the expected increase in the amplitude of the OIF and FM/W modes as decreases.

The energy of the magnetic fluctuations on scales of is quite similar regardless of the used mass ratio.

For finite mass ratios, the power spectra develop, via power cascade, a tail that behaves as:^{7}^{7}7This behavior differs from the one obtained from hybridPIC simulations (Kunz et al., 2014), where and . This possibly denotes the influence of the electrons in the cascading process. and .
The similar amplitude of the magnetic fluctuations on electron scales suggests that the break in the adiabatic invariance of can in principle be caused by the OIF and FM/W instabilities via a three steps scenario: ion and electron anisotropies create magnetic fluctuations through the OIF and FM/W instabilities, part of the energy in the magnetic fluctuations is transferred to electron scales via power cascade, and electrons are pitchangle scattered by these fluctuations producing the break in invariance. We can compare the contributions from the OIF and FM/W modes to the cascading process by looking at the power spectra by components: , , and . This is done in Figures 5 and 5 for the cases and 64, respectively. These figures show that and are comparable within the powerlaw tails. Since and are expected to be dominated by the OIF and FM/W modes, respectively, this result suggests that these two modes contribute similar amount of energy to the powerlaw tail.
A remaining question is whether the presented scenario is plausible in the more realistic case with , where a larger scale separation is expected between and . We explore this question using Figure 6, where we plot the growth rate as a function of for OIF and FM/W modes (green and black lines, respectively). This is done assuming , , and in two regimes: with such that the maximum growth rate (solid lines), and with such that (dotted lines). Since for each value of there are multiple OIF and FM/W wavevectors , in Figure 6 we chose the maximum and minimum for each . This way we explore the possibility that the modes with a given could have wavelengths close to both the electron and ion Larmor radii. We obtained that:

In the more realistic case with and , the fastest growing OIF and FM/W modes appear at . This implies that both modes have a factor larger wavelengths than in the , case, which is expected because of the factor increase in the ratio .
Thus, for the scale separation between and allows the generation of magnetic fluctuations at electron scales with enough energy to pitchangle scatter the electrons. This result relies on the existence of a power cascade with , which is observed in Figures 5 and 5. However, when this scenario seems less likely. Indeed, at electron scales () the cascade of OIF modes can produce an amount of energy times smaller than at the ion scales (; see dottedgreen line in Figure 6). However, at saturation the OIF and OEF modes are expected to satisfy ( and , respectively). Thus, in order to provide enough ion and electron pitchangle scattering, the energy at electron scales should be a factor smaller than at ion scales, which times larger than what can be produced through the cascade of OIF modes.^{8}^{8}8We did not include in this analysis the possible cascade of FM/W waves since, for realistic values of , their amplitude () should be much smaller than the one of the OIF modes (). Also, we did not include the possibility of electron scattering via cyclotron resonances (for which ), since we do not expect the cascade of OIF or FM/W modes to produce waves with the right polarization to resonate with electrons.
This difficulty may get ameliorated if the power cascade process were further modified when , or if in 3D the spectral index of the cascade powerlaw tail were different from the one obtained in our 2D simulations. Unfortunately, our current simulations can not clarify this aspect of the interplay between the electrons and the OIF and FM/W instabilities. It is important to point out, however, that in realistic settings we do not expect the OEF modes to produce the necessary electronscale fluctuations either, since our linear calculations show that these modes are stable for the electron anisotropy set by the OIF and FM/W instabilities (with ). Thus it seems likely that the electron anisotropy should continue to be determined by the OIF and FM/W marginal stability condition with .
4. Viscous Heating
The existence of electron and ion pressure anisotropies in general implies the presence of nondiagonal terms in the pressure tensor, which give rise to an effective viscosity for both species. In our case, particle velocities are nearly gyrotropic with respect to , so the relevant pressure tensor component is . It can be shown that this pressure component can tap into the velocity shear of the plasma, producing particle heating. In our case, assuming no heat flux, the internal energy density of species , (), evolves as (Kulsrud, 1983; Snyder et al., 1997; Sharma et al., 2007):
(3) 
where corresponds to the growth rate of . In the present context, both and are negative, which implies a positive particle heating.
Figure 7 quantifies the importance of this heating mechanism by showing the volumeaveraged ion (solidblack) and electron (solidgeen) heating rates for run F2. We also show the heating rate by anisotropic viscosity predicted by Equation 3 for ions (dottedblack) and electrons (dottedgreen). For both species there is reasonably good agreement between the particle heating in the simulation and the contribution from the anisotropic stress. This shows that anisotropic viscosity contributes most of the ion and electron heating in collisionless plasmas (with ), both in the case of decreasing magnetic field presented in this work, as well as in the growing field regime shown in Riquelme et al. (2016).
5. Electron Mean Free Path
Besides regulating the effective plasma viscosity, pitchangle scattering by velocityspace instabilities is also expected to limit the mean free path of the particles. To quantify this effect we used the distance traveled along by ions and electrons ^{9}^{9}9, where is the particle’s velocity.. We calculated their mean free paths assuming that (where represents the average mean free path over species , while is their thermal speed), which is valid if the particles move diffusively. This allows us to estimate the average mean free path of species as .
Figure 8 shows our estimates of the electron (black) and ion (red) mean free paths for two simulations with and (runs F2 and F1, respectively), normalized by . We see that in both cases there is an initial period when increases as , which is consistent with the free streaming behavior . This is followed by the saturation of , expected to start after a time of the order of the collision time of particles. The behaviors of and in this stage are quite similar, with being somewhat smaller than (by a factor ).
The evolutions of seen in Figure 8 are expected to be influenced by the pitchangle scattering of the particles caused by the velocityspace instabilities. In Riquelme et al. (2016) we estimated this effect assuming an incompressible fluid with no heat flux, where the scattering produced by the instabilities on species was incorporated using an effective scattering rate (Kulsrud, 1983; Snyder et al., 1997). In this model, valid for , behaves as:
(4) 
Equation 4 provides a good approximation to the mean free path of particles in the case of growing magnetic fields (Riquelme et al., 2016), where is regulated by the whistler and mirror instabilities. In the cases of decreasing magnetic field, however, this is not the case. Using our measurements of for runs F2 and F1 (see Figures 4 and 4) one obtains that, at , . However, Figure 8 shows that the average mean free paths of ions and electrons are times larger than this simple estimate. On the other hand, considering the evolution of and (needed to determine ) from to , should decrease by a factor , which is essentially what Figure 8 shows. Thus, putting aside the factor difference, the scaling of on and presented in Equation 4 is reasonably well reproduced by the simulations with decreasing .
The factor discrepancy is interesting, partly because the behavior of suggested by Equation 4 () is well reproduced in the case of ions in previous hybridPIC simulations that studied the saturated state of the firehose and mirror instabilities (Kunz et al., 2014). In that case, however, is not measured using . Instead, they constructed a distribution of the times taken by each ion to change its by a factor of e, and then approximated by the width of the distribution. Thus, in order to clarify this discrepancy, we compared two different measurements of the electron mean free path from run F1, one using the variations of (we will refer to this estimate as ) and the other one using (). These measurements of and are not defined for different times (as in Figure 8), but they correspond to averages between and 3.^{10}^{10}10For a given electron population, is estimated by first calculating for each electron the quantity (where is simply the time elapsed between and 3 and is the average magnitude of the electron velocity parallel to B in the same period), and then taking the average over the population. For , we construct a distribution of the times taken by each electron to change by a factor e with respect to its value at , and then is estimated multiplying the width of the distribution by . The comparison was made for different groups of electrons, defined by the parameter , where represents the average between and 3 for each electron. We use as a way to quantify the pitchangle of electrons, which, as we will see below, affects the behaviors of and in different ways.
Our results are shown in Figure 9. We see that and roughly coincide for (pitchangle ). However, for (pitchangle ), becomes significantly larger than . This is consistent with the fact that, for electrons with small pitchangle, an order unity variation in due to scattering (which implies an order unity variation in ) does not imply that they reverse their velocity along B, since . Also, when averaged over the entire distribution (shown by the black line in Figure 9), is a factor larger than , showing that using would essentially eliminate the discrepancy between our estimated mean free path and Equation 4.^{11}^{11}11The factor diference between our measured and the estimate given by Equation 4 was also obtained in simulations with infinite mass ions, where the electron scattering is dominated by the OEF modes. Also, the same difference is obtained in simulations with infinite mass ions and initial , suggesting that this discrepancy is fairly insensitive to , at least in the moderate regime that we studied.
Finally, Figure 9 shows the same quantities as Figure 9 but for run MW1 of Riquelme et al. (2016) (with growing mean magnetic field). We see that the trend for to grow relative to as decreases is maintained in this case. However, when averaged over the distribution of (black line), . The different behaviors of the growing and decreasing runs can be explained partly because particles with low pitchangle are more abundant in the decreasing field case (as can be seen comparing in Figures 9 and 9). However, regardless of , for all values of the decreasing case tends to have a larger ratio compared with the growing case. This trend is likely because, besides the effect of scattering, the pitchangle of particles is also affected by the permanent driving of the pressure anisotropy. Thus, in the case of decreasing , this driving tends to reduce the pitchangle of particles, possibly reducing their ability to reverse their velocities along B (compared with a population with similar and pitchangle, but in a growing ).
We notice that, both for growing and decreasing , the estimate of provided by Equation 4 is fairly well reproduced by using the change rate of to measure the scattering rate (as it was shown by Kunz et al. (2014) in the case of ions). However, since is based on the direct calculation of the distance travelled by the electrons along B, this quantity provides a more meaningful measurement of the electron mean free path for the purpose of quantifying the thermal conductivity of the plasma. Thus, for the decreasing field case, the scaling for the electron mean free path is best described by the relation:
(5) 
which is valid for , and only differs from Equation 4 by its prefactor instead of .
6. Discussion and Implications
We used particleincell (PIC) plasma simulations to study the nonlinear, saturated stage of various ion and electron velocityspace instabilities relevant for collisionless plasmas. We focused on instabilities driven by pressure anisotropy with . To capture the nonlinear regime in a selfconsistent way, we imposed a shear velocity in the plasma, which decreases the background magnetic field. This drives due to the adiabatic invariance of the magnetic moment (the driving timescale is much longer than the gyroperiod of the particles). This, in turn, drives velocityspace instabilities, which inhibit the growth of pressure anisotropy. The relevant instabilities in this regime, as suggested by linear theory, are: the purely growing oblique ionfirehose (OIF) and the resonant fast magnetosonic/whistler (FM/W) modes, which are mainly driven by the ions, and the purely growing oblique electronfirehose (OEF) and the resonant Alfvén/ioncyclotron (A/IC) modes, which are driven by the electrons. The nonlinear state of these instabilities is expected to be influenced by the simultaneous presence of ion and electron anisotropies on the different modes. In order to achieve reasonable scale separation between these modes, we mainly used and 64. Our results, valid for the regime , showed no significant difference between these two mass ratios.
We found that the mechanism for regulating the ion and electron anisotropies consists in the growth of OIF and FM/W modes, which affect equally the ions and electrons, rendering . The numerically obtained ion and electron anisotropies are well approximated by the linear threshold for the growth of the OIF and FM/W modes with and with growth rate . The electron pressure anisotropy in simulations with infinite mass ions (where the ions only provide a neutralizing charge) is dominated by the OEF and A/IC modes, giving a factor larger anisotropy than in the cases with finite . We attribute this result to the rather strong destabilizing effect of the electron pressure anisotropy, , on the OIF and FM/W modes, which in turn maintains at a value significantly lower than the one necessary to make the OEF and A/IC modes grow at a rate .
Although the amplitude of the OIF and FM/W modes depend on the ratio , the value of the parameters and used in the simulations do not affect our conclusions. Also, based on our linear Vlasov calculations (Verscharen et al., 2013), we infer that the presented scenario should hold in the , highly magnetized () case relevant for real astrophysical plasmas. However, an important point that our simulations could not completely clarify (due to the lack of sufficient ion and electron scale separation) is the mechanism by which the electrons would be pitchangle scattered in the saturated stage of the OIF and FM/W instabilities in the case of . Answering this question completely requires using significantly larger mass ratios and magnetizations (and possibly 3D runs), so we differ this aspect of the study to a future work.
We have also used our simulations to verify the expected viscous heating of particles, which is described in Equation 3, and arises due to pressure anisotropies tapping into the free energy in the shear motion of the plasma. Figure 7 shows a good agreement between the heating of the particles in our simulations and the expectation from Equation 3. This result is valid for decreasing magnetic fields, as shown here, and for growing fields as shown by Riquelme et al. (2016).
With the intention to quantify the thermal conductivity in these plasmas, we have also computed the mean free path of species , , during the nonlinear stage of the OIF and FM/W instabilities. The average mean free path of both ions and electrons is reasonably well described by Equation 5. The scaling factors in this equation are the same as in Equation 4, which is based on a model where the mean free path of species is determined by an effective scattering rate that sets the rate at which the invariance is broken (Kulsrud, 1983; Snyder et al., 1997). However, the prefactor in Equation 5 is times larger than in the case of Equation 4. We explained this discrepancy for the case of the electrons by noticing that the variation of provides a good estimate of only for electrons with relatively large pitchangle (). For electrons with pitchangles smaller than (which are abundant in the case of decreasing mean magnetic field), can become a factor larger. This is in contrast with the case of growing (Riquelme et al., 2016), where Equation 4 provides a reasonably good estimate for , despite the fact that also grows as the pitchangle decreases. This difference is likely because in our present setup there are significantly more particles with low pitchangle, and also because the permanent driving of pressure anisotropy may decrease the ability of particles to reverse velocity along B (compared with the case of growing ).
The results shown in this work, as well as those presented in Riquelme et al. (2016) for the growing case, are relevant for quantifying the viscous heating and thermal conductivity in various lowcollisionality astrophysical plasmas, including lowluminosity accretion flows around compact objects, the ICM, and the heliosphere.
References
 Buneman (1993) Buneman, O. 1993, “Computer Space Plasma Physics”, Terra Scientific, Tokyo, 67
 Gary et al. (1998) Gary, S. P., Li, H., O?Rourke, S. & Winske, D. 1998, J. Geophys. Res., 103, 14,567
 Kulsrud (1983) Kulsrud, R. M. 1983, in Handbook of Plasma Physics, ed. M. N. Rosenbluth & R. Z. Sagdeev (Amsterdam: North Holland), 115
 Hellinger et al. (2000) Hellinger, P. & Matsumoto, H. 2000, J. Geophys. Res., 105, 10,519
 Komarov et al. (2016) Komarov, S. V., Churazov, E. M., Kunz, M. W., & Schekochihin, A. A. 2016, MNRAS, 460, 467
 Kunz et al. (2014) Kunz, M. W., Schekochihin, A. A., & Stone, J. M. 2014, Physical Review Letters, 112, 205003
 Li et al. (2000) Li, X. & Habbal, S. R. 2000, J. Geophys. Res., 105, 27377
 Lyutikov (2007) Lyutikov, M. 2007, ApJ, 668, L1
 Marsch (2006) Marsch, E. 2006, Living Rev. Solar Phys., 3, 1
 Maruca et al. (2011) Maruca, B. A., Kasper, J. C. & Bale, S. D. 2011, Phys. Rev. Lett. 107, 201101
 Quest et al. (1996) Quest, K. B. & Shapiro, V. D. 1996, J. Geophys. Res., 101,24,457
 Remya et al. (2013) Remya, B., Reddy, R. V., Tsurutani, B. T., Lakhina, G. S., & Echer, E. 2013, JGRA, 118, 785
 Riquelme et al. (2015) Riquelme, M. A., Quataert, & Verscharen, D. 2015, ApJ, 800, 27
 Riquelme et al. (2016) Riquelme, M. A., Quataert, & Verscharen, D. 2016, ApJ, 824, 123
 Schekochihin et al. (2005) Schekochihin, A. A., Cowley, S. C., Kulsrud, R. M., Hammett, G. W., & Sharma, P. 2005, ApJ, 629, 139
 Sharma et al. (2006) Sharma, P., Hammett, G. W., Quataert, E., & Stone, J. 2006, ApJ, 637, 952
 Sharma et al. (2007) Sharma, P., Quataert, E., Hammett, G. W., & Stone, J. 2007, ApJ, 667, 714
 Snyder et al. (1997) Snyder, P. B., Hammett, G. W., & Dorland, W. 1997, Phys. Plasmas, 4, 3974
 Spitkovsky (2005) Spitkovsky, A. 2005, AIP Conf. Proc, 801, 345, astroph/0603211
 Squire et al. (2017) Squire, J., Kunz, M., Quataert, E., & Schekochihin, A. 2017, arXiv:1705.01956
 Verscharen et al. (2013) Verscharen, D., Bourouaine, S., Chandran, B. D. G., & Maruca, B. A. 2013, ApJ, 773, 8