PIC Simulations of Velocity-Space Instabilities in a Decreasing Magnetic Field: Viscosity and Thermal Conduction

We use particle-in-cell (PIC) simulations of a collisionless, electron-ion plasma with a decreasing background magnetic field, $B$, to study the effect of velocity-space instabilities on the viscous heating and thermal conduction of the plasma. If $B$ decreases, the adiabatic invariance of the magnetic moment gives rise to pressure anisotropies with $p_{||,j}>p_{\perp,j}$ ($p_{||,j}$ and $p_{\perp,j}$ represent the pressure of species $j$ ($=i$ or $e$) parallel and perpendicular to the magnetic field). Linear theory indicates that, for sufficiently large anisotropies, different velocity-space instabilities can be triggered. These instabilities, which grow on scales comparable to the electron and ion Larmor radii, in principle have the ability to pitch-angle 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 the magnetic field by an imposed shear in the plasma. Our results show that, in the regime $2 \lesssim \beta_j \lesssim 20$ ($\beta_j \equiv 8\pi p_j/B^2$), 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: $\Delta p_e/p_{||,e} \approx \Delta p_i/p_{||,i}$ (where $\Delta p_j=p_{\perp,j} - p_{||,j}$). We also quantify the thermal conduction of the plasma by directly calculating the mean free path of electrons along the mean magnetic field, which we find strongly depends on whether $B$ 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.


INTRODUCTION
In low-collisionality plasmas, the change in the magnitude of the local magnetic field (B ≡ |B|) generically drives a pressure anisotropy with p ,j = p ⊥,j (where p ⊥,j and p ,j correspond to the pressure of species j perpendicular and parallel to B). This is a consequence of the adiabatic invariance of the magnetic moment of particles, µ j ≡ v 2 ⊥,j /B, in the absence of collisions (where v ⊥,j is the velocity of species j perpendicular to B).
These pressure anisotropies can trigger various velocityspace instabilities, which are in principle expected to pitch-angle scatter the particles, to some extent mimicking the effect of collisions. The combined effect of pressure anisotropies and velocity-space 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 low-luminosity 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 B, which naturally drives p ⊥,j > p ||,j . In this paper we study the opposite case, where B decreases and p ⊥,j < p ||,j . In this case several velocity-space instabilities can be excited, ultimately regulating the extent to which the p ⊥,j < p ||,j 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: i) the oblique electron firehose (OEF) modes, which are purely growing modes, and ii) the Alfvén/ion-cyclotron (A/IC) modes, which are quasi-parallel, propagating waves, driven unstable by cyclotron-resonant electrons (Li et al. 2000;Camporeale et al. 2008). 5 Similarly, in presence of an ion pressure aniostropy p ⊥,i < p ||,i , there are also two types of modes that can grow unstable: i) the oblique ion firehose (OIF) modes, which are purely growing modes, and ii) the fast-magnetosonic/whistler (FM/W) modes, which are quasi-parallel, propagating waves, excited by cyclotron-resonant ions (Quest et al. 1996;Gary et al. 1998;Hellinger et al. 2000). 6 5 Although the Alfvén/ion-cyclotron (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 (k) Alfvén mode and then becomes the ioncyclotron mode at higher k. 6 Although gyrokinetic theory suggests that OIF and OEF modes correspond to the same Alfvén-mode branch (Kunz et al. 2015;Verscharen et al. 2017), we will consider them as separate instabilities, identifying them as two different growth-rate maxima in k-space.
In this work we studied the nonlinear, saturated properties of these instabilities making use of particle-in-cell (PIC) simulations. This is achieved by continuously decreasing the strength of the background magnetic field by externally imposing a shear motion in the plasma. This setup is interesting since in realistic astrophysical scenarios the pressure anisotropies are expected to be driven (via B decrease) for a time significantly longer than the initial regime where the instabilities grow exponentially.
Previous works have already studied this long term regime by simulating an expanding (instead of shearing) plasma.
These works have used both hybrid-PIC simulations, which focused on the evolution of the ion anisotropy-driven instabilities (Matteini et al. 2006;Hellinger et al. 2008), and PIC simulations that mainly captured the role of electron anisotropy-driven modes (Camporeale et al. 2010). Thus, our work is intended to study the combined effect of the electron and ion pressure anisotropies on the nonlinear, saturated regime of the different unstable modes. This aspect of our study is motivated in part by previous linear dispersion analyses showing that the electron pressure anisotropy can significantly influence the evolution of both the FM/W and OIF modes (Michno et al. 2014;Maneva et al. 2016).
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 low-collisionality plasmas (Sharma et al. 2006(Sharma et al. , 2007Squire et al. 2017). Also, the nonlinear evolution of the different velocity-space instabilities sets the pitch-angle 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 ∆p j 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 low-collisionality astrophysical plasmas.

SIMULATION SETUP
We use the electromagnetic, relativistic PIC code TRISTAN-MP (Buneman 1993;Spitkovsky 2005) in two dimensions. The simulation box consists of a square box in the x-y plane, containing an initially isotropic plasma with a homogeneous initial magnetic field B 0 . We simulate a decreasing magnetic field by imposing a velocity shear given by v = −sxŷ, where s is the shear rate of the plasma and x is the distance alongx (x andŷ are the unit vectors parallel to the x and y axes, respectively). From flux conservation, the x and y components of the mean field evolve as d B Note. -A summary of the physical and numerical parameters of our simulations. These are the mass ratio m i /me, the initial electron magnetization ωc,e/s, the ratio between electron temperature and rest mass energy, kTe/mec 2 , the number of particles per cell Nppc (including ions and electrons), and the box size in units of the typical initial electron Larmor radius L/R L,e (R L,e = v th,e /ωc,e, where v 2 th,e = k B Te/me, k B is Boltzmann's constant, and Te is the electron temperature). All of the runs initially have β i = βe = 2 and c = 0.225∆x/∆t, where ∆t is the simulation time step. The runs shown in this Table share the same electron skin depth c/ωp,e/∆x = 5 (where ∆x is the grid point separation), but we used several other simulations to confirm numerical convergence by varying c/ωc,e/∆x, Nppc, and L/R L,e .
positive, there will be a decrease of B y and, therefore, of | B |. Therefore, we initially choose B 0 ∝x + 3.3ŷ, which garantees a decrease of | B | and a p ⊥,j < p ||,j anisotropy during a simulation time ∼ 3s −1 .
By resolving the x-y plane, our simulations can capture the quasi-parallel 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 B . 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, ω c,j /s (j = i, e), and the ion to electron mass ratio, m i /m e . In typical astrophysical cases, ω c,j ≫ s and m i /m e = 1836. Due to computational constraints, however, we will use values of ω c,j /s and m i /m e 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 k B T e /m e c 2 = 0.28, which implies ω c,e /ω p,e = 0.53 (where k B , T e , and ω p,e are Boltzmann's constant, the electron temperature, and the electron plasma frequency, respectively). We will change our simulation conditions by varying: ω c,e /s and m i /m e (which uniquely fix ω c,i /s and k B T i /m i c 2 ). 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 electron-scale physics. These provide a useful contrast with our finite m i /m e runs and allow us to isolate the impact of ion physics on the electrons. The numerical parameters in our simulations will be: N ppc (number of particles per cell), c/ω p,e /∆ x (the electron skin depth in terms of grid size), L/R L,i (box size in terms of the initial ion Larmor radius for runs with finite -The initial evolution of the electron pressures perpendicular (p ⊥,e ; black-solid) and parallel (p ||,e ; red-solid) to B for runs I1 and I2 in Table 1 (with ωc,e/s = 3600 and ωc,e/s = 7200, respectively). The black-and red-dotted lines show the expected p ⊥,e and p ||,e evolutions from the CGL or double adiabatic limit (Chew et al. 1956). A significant deviation from adiabatic evolution can be seen at t · s 0.7.
, and L/R L,e (box size in terms of the initial electron Larmor radius for runs with infinite m i /m e ). Table  1 shows a summary of our key runs. We ran a series of simulations ensuring that the numerical parameters (e.g., different N ppc ) do not significantly affect our results. Note that most runs used just for numerical convergence are not in Table 1.

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 m i /m e = ∞ Figure 1 shows the early time evolution (until t · s ≈ 1.3) of the electron pressures perpendicular (p ⊥,e ; black-solid line) and parallel (p ||,e ; red-solid line) to B for runs I1 and I2 of Table 1. These runs have infinite mass ions so that the electrons can only be affected by the electron anisotropy-driven OEF and A/IC instabilities, and their magnetizations are ω c,e /s = 3600 and 7200, respectively. The black-and red-dotted lines show the expected evolutions of p ⊥,e and p ||,e from the CGL or double adiabatic limit (Chew et al. 1956). We see that this adiabatic limit is reasonably well satisfied in the early stage of the simulations (until t · s ∼ 0.7), regardless of the magnetization ω c,e /s. After that, the growth of the electron anisotropy-driven instabilities provides enough pitch-angle scattering to stop the adiabatic evolution of the electron pressure.
The presence of the OEF and A/IC instabilities can be seen from Figure 2, which shows the magnetic field fluctuations and plasma density in simulation I1. The upper row in Figure 2 corresponds to t · s = 1, i.e., after one shear time, while the lower row corresponds to t · s = 2. At all times the magnetic fluctuations are dominated by their δB z component, with wavenumbers k (≡ |k |, where k is the wave vector) satisfying kR L,e ∼ 0.2, and with k being mainly oblique to the mean direction of B. The presence of these oblique modes can also be seen from Figure 3b, which shows the Fourier transform of δB z at t · s = 1 as a function of the wavenumbers parallel and perpendicular to B (k || and k ⊥ , respectively). This suggests that the OEF modes contribute the most to the amplitude of the magnetic fluctuations. However, although smaller in amplitude, quasi-parallel A/IC modes can also be seen especially in the δB x component (see Figure 2a). This is also seen from Figure 3a, which shows that the Fourier transform of δB x (as a function of k || and k ⊥ ) is dominated by quasi-parallel modes. Figure 4 shows the time evolution of the magnetic energy in B , the volume-averaged pressure anisotropy, and the electron magnetic moment for two runs, one with ω c,e /s = 3600 and the other with ω c,e /s = 7200 (runs I1 and I2 in Table 1, respectively). Panels c and d show the volume-averaged pressure anisotropy − ∆p e / p ||,e for these two runs, where ∆p e = p ⊥,e − p ||,e . 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, s. These thresholds were calculated using the linear Vlasov solver developed by Verscharen et al. (2013). The calculations use mass ratio m i /m e = 1836 and assume very cold ions (β i = 10 −4 ), which seeks to resemble our simulated m i /m e = ∞ situation where the ions only provide a neutralizing charge. 7 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 ω c,e /s = 3600. This implies that both modes should play some role in regulating the electron anisotropy, with their relative importance depending weakly on the ratio ω c,e /s. Also, for both values of ω c,e /s 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 s.
The contribution of the different components of δB can be seen from Figures 4a and 4b, which show the magnetic energy along different axes as a function of time, normalized by the average magnetic energy in the simulation, B 2 /8π. δB is decomposed in terms of δB z (component perpendicular to the simulation plane), δB xy,⊥ (component parallel to the simulation plane but perpendicular to B ), and δB || (component parallel to B ). Clearly, δB is dominated by its z component (as already seen in Figure 2). 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 δB is quite different. Indeed, our linear calculations show that, for the plasma parameters of runs I1 and I2, the OEF modes should satisfy |δB z |/(|δB || | + |δB xy,⊥ |) ∼ 4. Thus the fact that in our runs δB 2 z ∼ 10δB 2 xy,⊥ implies that 7 In order to make sure that using mildly relativistic electrons in our runs (where k B Te/mec 2 = 0.28) does not invalidate our comparison with the calculated thresholds (which assume nonrelativistic electrons), in Figure 4d we added in solid-red line the pressure anisotropy for run I3, which uses k B Te/mec 2 = 0.1 (while keeping the same ωc,e/s and βe). The very small difference between the two cases suggests that the effect of having mildly relativistic electrons is fairly small. -The three components of δB and plasma density fluctuations δρ at two different times: t · s = 1 (upper row) and t · s = 2 (lower row), for a simulation with infinite ion mass (run I1). Fields and density are normalized by the initial magnetic field, B 0 , and the average density, ρ 0 . The arrows in panels d and h show the magnetic field direction on the x − y simulation plane. For this m i /me = ∞ case, the magnetic fluctuations are dominated by the oblique OEF modes. The presence of the quasi-parallel A/IC and OEF modes for run I1 is shown by the magnitude of the Fourier transform of δBx (panel a) and δBz (panel b) at t · s = 1. These quantities are plotted as a function of the wavenumbers parallel and perpendicular to B (k || and k ⊥ , respectively), are normalized by their maximum value, and are raised to the 1/5th power to provide better dynamical range. The contributions from (quasiparallel) A/IC and (oblique) OEF modes with wavevectors satisfying kR L,e ∼ 0.2 are most clearly seen from panels a and b, respectively. most of the magnetic energy is being contributed by the OEF modes. The quasi-parallel A/IC modes, which are most visible in the δB x component as can be seen from Figure 2a, make a significantly smaller contribution to δB. Another sign of the dominance of the OEF modes is the growing and damping phases of δB z observed in Figures 4a and 4b, which are likely related to the conversion of the saturated OEF modes into propagating waves that are rapidly damped through scattering with electrons, as it has been observed in previous initial value PIC simulations (e.g. Hellinger et al. 2014). The different contributions to δB 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 pitch-angle gradients of the distribution function. Therefore, the A/IC instability can relax the distribution faster toward a stable configuration through pitch-angle scattering than the OEF instability.
Finally, Figures 4e and 4f show the volume-average magnetic moment of the electrons, µ e (≡ p ⊥,e /B ; black-solid 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 (t · s ≈ 0.7), µ e is fairly constant, implying the lack of efficient pitch-angle scattering. After that, µ e tends to increase at a rate close to the shear rate s. This implies the appearance of an effective pitch-angle scattering rate for the electrons, ν ef f,e , due to their interaction with the OEF and A/IC modes.
In order to help us to understand the way ∆p j is regulated by the different velocity-space instabilities, we propose a second way to calculate the average magnetic moment of species j: (1) ωc,e/s = 7200 (run I2; right column), which use m i /me = ∞ (the ions simply provide a neutralizing charge). Upper row: the volumeaveraged magnetic energy components, δB 2 r pointing: along the z axis (δB 2 z ; green), parallel to B (δB 2 || ; black), and perpendicular to B but in the x-y plane (δB 2 xy,⊥ ; red), normalized by B 2 . Middle row: the evolution of the electron pressure anisotropy (green line), with the linear OEF and A/IC instabilities thresholds for growth rates γ = s (dotted-green and dotted-black lines, respectively). The pressure anisotropy saturates at a value roughly consistent with the linear OEF and A/IC instabilities growing at the rate γ = s. Additionally, panel d) shows in solid-red line the electron anisotropy for run I3, with the same parameters as run I2 but with kTe/mec 2 = 0.1 instead of 0.28. The small difference between the solid-green and solid-red lines shows that the electrons being mildly relativistic should not affect substantially the evolution of ∆pe. Lower row: the electron magnetic moment; see equation 1 and associated discussion for definitions of µe (solid) and µ ef f,e (dotted).
This definition is useful because there can be cases where µ j,ef f = µ j . This is expected when, besides pitchangle scattering, ∆p j is partly regulated by relatively large fluctuations in B, which may spatially correlate with p ⊥,j in a µ j conserving way. This occurs, for instance, in the presence of large amplitude mirror modes (Kunz et al. 2014;Riquelme et al. 2015Riquelme et al. , 2016. Figures  4e and 4f show that, after µ e conservation is broken, µ e,ef f (dotted-black) and µ e (solid-black) are almost indistinguishable. This confirms that ∆p e is regulated by an effective pitch-angle scattering provided by the OEF and A/IC modes, with no significant contribution from fluctuations in B.

3.2.
Simulations with finite m i /m e 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 m i /m e . Since using m i /m e ≃ 1836 is computationally infeasible with our current resources, we have instead tried to ensure that our simulation results do not depend significantly on m i /m e , which is reasonably well achieved for m i /m e = 25 and 64.
As an example, Figure 5 shows the three components of δB for run F1 of Table 1 (m i /m e = 64 and ω c,e /s = 7200). The upper and lower rows correspond to t · s = 1 and t · s = 2, respectively. At both times, quasi-parallel and oblique modes are present with similar amplitudes, and with wavenumbers satisfying kR L,i ∼ 0.4 (where R L,i is the ion Larmor radius). While the quasi-parallel modes are apparent in the three components of δB, the oblique modes appear mainly in δB z . This can be seen more clearly from Figure 6, which shows the Fourier transform of δB x (Fig. 6a) and δB z (Fig. 6b) at t · s = 1 as a function of k || and k ⊥ . The presence of quasi-parallel modes is clear in both panels, while the oblique modes mainly appear in δB z . 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 Figures 5 and 6. We will come back to this question below.
3.2.1. The Role of mi/me Figure 7 compares the evolution of the energy in δB, the ion and electron anisotropies, and µ i and µ e for simulations with different mass ratios and electron magnetization. The first and second columns compare simulations with the same electron conditions (ω c,e /s = 7200, k B T e = 0.28m e c 2 , and β e = β i = 2) but with different mass ratios: m i /m e = 25 (run F2) and m i /m e = 64 (run F1), respectively.
Figures 7d and 7e show the volume-averaged 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 s. In dotted lines we show thresholds that assume ∆p e = ∆p i , and correspond to the instabilities: OEF (dotted-green), OIF (dotted-red), 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 cyclotron-resonant ion damping when T i ∼ T e , becoming stable in our runs with finite m i /m e (this is true even for m i /m e = 1836). 8 We see that our simulation results are fairly independent of the mass ratio, and can be summarized as follows: 1. The ion and electron anisotropies evolve quite similarly in both cases. (This justifies using instability  for run F1 (m i /me = 64) at t · s = 1, as a function of k || and k ⊥ and normailized by their maximum value (these quantities are raised to the 2/5th power to provide better dynamical range). The presence of the (quasi-parallel) FM/W modes with kR L,e ∼ 0.4 is clearly seen in both panels. The presence of the (oblique) OIF modes is apparent mainly in δBz .
thresholds that assume ∆p i = ∆p e for comparison.) 2. The obtained electron anisotropy is a factor ∼ 1.5 smaller than the expected OEF threshold.
3. The ion and electron anisotropies are best described by the thresholds of the OIF and FM/W instabilities with ∆p i = ∆p e .
4. The OIF and FM/W thresholds with ∆p i = ∆p e are quite similar, which is consistent with the simultaneous presence of these modes in Figure 5.
The fact that the electron anisotropy is close to the OIF and FM/W thresholds, and a factor ∼ 1.5 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 ∆p i thresholds for the OIF and FM/W instabilities with ∆p e = 0 are a factor of ∼ 2 larger than in the ∆p i = ∆p e 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 m i /m e = 1836 with ω c,e /γ = 10 6 (where γ represents the growth rate of the different instabilities), which are shown in Figure 9a. We find that the thresholds of the OIF (dotted-red) and FM/W (dotted-blue) modes with ∆p e = ∆p i continue to be similar and smaller than the OEF (dotted-green) threshold by a factor ∼ 1.5 (the dotted-green line almost coincides with the solid-red line). Also, the OIF and FM/W thresholds with ∆p e = ∆p i are ∼ 1.5 times smaller than in the ∆p e = 0 case. This implies that, for realistic mass ratios and magnetizations, the dominant instability for . The upper row shows the volume-averaged magnetic energy in three components of δB : δBz (green), δB || (black), and δB xy,⊥ (red), normalized by B 2 /8π. The middle row shows the ion (black) and electron (green) pressure anisotropies, ∆p j /p ||,j , along with the anisotropy thresholds for the growth of different instabilities (using a growth rate of s). The dotted lines show thresholds assuming ∆pe = ∆p i , which correspond to the instabilities: OEF (dotted-green), OIF (dotted-red), and FM/W (dotted-blue). The solid-red and solid-blue lines also show the threshold for the growth of the OIF and FM/W instabilities, respectively, assuming ∆pe = 0. Our results are consistent with the OIF and FM/W modes dominating the pressure anisotropies of ions and electrons, which satisfy ∆pe ≈ ∆p i . The lower row shows the volume-averaged ion (solid-red) and electron (solid-black) magnetic moments, as well as the "effective" averages defined as in equation 1 for ions (dotted-red) and electrons (dotted-black), and normalized by the initial value of µ j .
the regulation of ion and electron anisotropies should continue to be the OIF and FM/W instabilities. 9 9 Although 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 (t · s 1.8, β i 7), while the opposite situation happens at late times (t · s 1.8, β i 7). This implies that in the more realistic cases there could be a clearer dominance of the FM/W (OIF) modes for β i 7 ( 7). This is consistent with the OIF dominance shown in the hybrid-PIC simulations with    Verscharen et al. (2013) and assuming the same β j evolution of our simulations (γ represents the growth rate for the different instabilities). The cases ∆pe = ∆p i are shown in dotted lines and ∆pe = 0 in solid lines. The OEF, OIF, and FM/W instabilities are represented by green, red, and blue colors. For these more realistic parameters, the OIF and FM/W instabilities are expected to continue to dominate in the regulation of both ion and electron anisotropies, rendering ∆pe ≈ ∆p i . Panel b: The growth rate γ as a function of kR L,e for OIF and FM/W modes (green and black lines, respectively), assuming β i = βe = 10, ∆pe = ∆p i . We consider two regimes: i) m i /me = 64 with ∆p j such that the maximum growth rate γmax = ωc,e/7200 (solid lines), and ii) m i /me = 1836 with ∆p j such that γmax = ωc,e/10 6 (dotted lines). We choose the maximum and minimum |k| = k for each γ.
64, respectively (for comparison the case m i /m e = 25 is replicated in Figure 7b using dotted lines). For the two mass ratios, the two components perpendicular to δB dominate, with δB 2 z being most of the time ∼ 3 times larger than δB 2 xy,⊥ . This result implies that the OIF and FM/W modes are contributing comparable energy to δB, as was also noticeable from Figure 5. Indeed, our linear calculations show that for the OIF modes δB 2 z ≫ δB 2 xy,⊥ , while for the FM/W modes δB 2 z ∼ δB 2 xy,⊥ , 10 which implies that most of the δB xy,⊥ 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 δB 2 z in the m i /m e = 64 case is on average ∼ 1.5 times larger than in the m i /m e = 25 case. Since the magnetizations ω c,i /s of the two runs differ by a factor ∼ 2.6 (= 64/25), this is roughly consistent with previous studies of the OIF instability that show that δB 2 at saturation should scale as δB 2 /B 2 ∝ (s/ω c,i ) 1/2 (Kunz et al. 2014). On the other hand, δB 2 xy,⊥ in the m i /m e = 64 case is about ∼ 2 times larger than for m i /m e = 25, which is roughly consistent with the expectation for the FM/W modes to have a saturation amplitude that satisfies δB 2 /B 2 ∝ s/ω c,i . This scaling can be obtained from the expected effective ion scattering frequency by resonant waves, ν ef f,i , which scales as where k || and v || are the wave vector component and the particle velocity component parallel to B (Marsch 2006). For the case of the quasi-parallel FM/W waves, we obtained from the simulations that k || R L,i ∼ 0.3, meaning that for most particles k || v || ∝ ω c,i and that ν ef f,i ∝ (δB 2 /B 2 )ω c,i . Thus, since at FM/W saturation one expects ν ef f,i ∼ s · p ||,i /∆p i (see Equation 4 below), this implies that δB 2 /B 2 ∝ s/ω c,i (considering that the change in ∆p i /p ||,i between the m i /m e = 25 and 64 cases is small).
Finally, in panels 7g and 7h we compare µ j = p ⊥,j /B and µ j,ef f = p ⊥,j / B for both ions and electrons. We see that the change in µ j,ef f tends to be somewhat larger than the one in µ j (by ∼ 20%) for the two mass ratios tested. This implies that the combined effect of the OIF and FM/W modes reduces p ⊥,j in a way that mainly breaks the adiabatic invariance of µ j , with the preservation of µ j due to changes in the field configuration playing a small role.

The Role of ωc,e/s
It is also important to understand the role of electron magnetization, ω c,e /s, while keeping the same mass ratio. This can be done by looking at the first and third columns in Figure 7, which compares simulations with m i /m e = 25 but with electron magnetization ω c,e /s = 7200 and 3600 (runs F2 and F3 in Table 1, 10 Indeed, using the Vlasov solver of Verscharen et al. (2013) one can obtain that, for the parameters of runs F1 (m i /me = 64) and F2 (m i /me = 25), the OIF modes satisfy |δBz |/(δB 2 xy,|| + δB 2 xy,⊥ ) 1/2 ∼ 5. respectively). We see that the two runs reproduce essentially the same results in terms of µ j and ∆p j evolution, with the only difference being in the amplitude of the magnetic fluctuations. Apart from some significant time variability, the amplitude of the δB z and δB xy,⊥ components in the m i /m e = 25, ω c,e /s = 3600 run are quite similar to the case m i /m e = 64, ω c,e /s = 7200. Since these two runs have a very similar ratio ω c,i /s (= 144 and 113, respectively), this result is consistent with the dependence of the OIF and FM/W saturated amplitude on ω c,i /s mentioned in §3.2.1.

Breaking of µe Adiabatic Invariance
An important question is whether the ion-scale instabilities alone are capable of explaining the break in the electron magnetic moment shown in Figure 7g, 7h, and 7i (which starts at t · s ∼ 0.7).
We explore this issue by comparing the power spectra of the fluctuations in our finite m i /m e runs with the power spectrum produced in the case with m i /m e = ∞. This is done in Figures 8a and 8b, where the magnetic energy per logarithmic unit of k || and k ⊥ is plotted at t · s = 2 for runs with m i /m e = 10, 25, 64, and ∞ (runs F4, F2, F1, and I2, respectively; k || and k ⊥ are the magnitude of the wave vector components parallel and perpendicular to B , respectively). The electrons in these simulations have the same conditions (k B T e /m e c 2 = 0.28, ω c,e /s = 7200, and initial β e = 2), so the different results are only due to the different m i /m e . We see that: 1. In the cases with finite mass ratio, as m i /m e increases the peaks of the spectra shift to longer wavelengths (in units of R L,e ) in a way consistent with the growth of the R L,i /R L,e ratio.
2. In the same way, as m i /m e 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 ω c,i /s decreases.
3. The energy of the magnetic fluctuations on scales of k ⊥ R L,e ∼ k || R L,e ∼ 0.2 is quite similar regardless of the used mass ratio.
4. For finite mass ratios, the power spectra develop, via power cascade, a tail that behaves as: 11 The similar amplitude of the magnetic fluctuations on electron scales suggests that the break in the adiabatic invariance of µ e can in principle be caused by the OIF and FM/W instabilities via a three steps scenario: i) ion and electron anisotropies create magnetic fluctuations through the OIF and FM/W instabilities, ii) part of the energy in the magnetic fluctuations is transferred to electron scales via power cascade, and iii) electrons are 11 This behavior differs from the one obtained from hybrid-PIC simulations (Kunz et al. 2014), where dδB 2 /dln(k || ) ∝ k −3.8 || and dδB 2 /dln(k ⊥ ) ∝ k −3.8 ⊥ . This possibly denotes the influence of the electrons in the cascading process. pitch-angle scattered by these fluctuations producing the break in µ e 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: δB 2 z , δB 2 || , and δB 2 xy,⊥ . This is done in Figures 8c and 8d for the cases m i /m e = 25 and 64, respectively. These figures show that δB 2 z and δB 2 xy,⊥ are comparable within the power-law tails. Since δB 2 z and δB 2 xy,⊥ 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 power-law tail.
A remaining question is whether the presented scenario is plausible in the more realistic case with m i /m e = 1836, where a larger scale separation is expected between R L,i and R L,e . We explore this question using Figure 9b, where we plot the growth rate γ as a function of kR L,e for OIF and FM/W modes (green and black lines, respectively). This is done assuming β i = β e = 10, ∆p e = ∆p i , and in two regimes: i) m i /m e = 64 with ∆p j such that the maximum growth rate γ max = ω c,e /7200 (solid lines), and ii) m i /m e = 1836 with ∆p j such that γ max = ω c,e /10 6 (dotted lines). Since for each value of γ there are multiple OIF and FM/W wavevectors k, in Figure 9b we chose the maximum and minimum |k| = k 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: 1. In the case with m i /m e = 64 and γ max = ω c,e /7200, the fastest growing OIF and FM/W modes have 0.03 kR L,e 0.06. This range roughly coincides with the peak of the spectra shown in Figures 8a and 8b. 2. In the more realistic case with m i /m e = 1836 and γ max = ω c,e /10 6 , the fastest growing OIF and FM/W modes appear at 0.006 kR L,e 0.01. This implies that both modes have a factor ∼ 6 larger wavelengths than in the m i /m e = 64, γ max = ω c,e /7200 case, which is expected because of the factor ∼ 6 increase in the ratio R L,i /R L,e .
Thus, for m i /m e = 64 the scale separation between R L,i and R L,e allows the generation of magnetic fluctuations at electron scales with enough energy to pitch-angle scatter the electrons. This result relies on the existence of a power cascade with dδB 2 /dln(k) ∝ k −2.8 , which is observed in Figures 8a and 8b. However, when m i /m e = 1836 this scenario seems less likely. Indeed, at electron scales (kR L,e ∼ 0.2) the cascade of OIF modes can produce an amount of energy ∼ (0.2/0.01) 2.8 ∼ 4400 times smaller than at the ion scales (kR L,i ∼ 0.01; see dotted-green line in Figure 9b). However, at saturation the OIF and OEF modes are expected to satisfy δB 2 /B 2 ∝ (s/ω c,j ) 1/2 (j = i and e, respectively). Thus, in order to provide enough ion and electron pitch-angle scattering, the energy at electron scales should be a factor ∼ (ω c,e /ω c,i ) 1/2 ∼ 1836 1/2 ∼ 43 smaller than at ion scales, which is ∼ 100 times larger than what can be for ions and electrons are shown with the dotted-black and the dotted-green lines, respectively. All quantities are normalized by p 0 s, where p 0 is the initial particle pressure in the simulation. For both species there is reasonably good agreement between the particle heating in the simulation and the contribution from the anisotropic stress.
produced through the cascade of OIF modes. 12 This difficulty may get ameliorated if the power cascade process were further modified when m i /m e = 1836, or if in 3D the spectral index of the cascade power-law 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 electron-scale 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 ∆p e = ∆p i ). Thus it seems likely that the electron anisotropy should continue to be determined by the OIF and FM/W marginal stability condition with ∆p e = ∆p i .

VISCOUS HEATING
The existence of electron and ion pressure anisotropies in general implies the presence of non-diagonal 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 B , so the relevant pressure tensor component is p xy,j ∝ (p ⊥,j − p ,j )B x B y /B 2 . It can be shown that this pressure component can tap into 12 We did not include in this analysis the possible cascade of FM/W waves since, for realistic values of s/ω c,i , their amplitude δB 2 /B 2 (∝ s/ω c,i ) should be much smaller than the one of the OIF modes (δB 2 /B 2 ∝ (s/ω c,j ) 1/2 ). Also, we did not include the possibility of electron scattering via cyclotron resonances (for which δB 2 /B 2 ∝ s/ωc,e), since we do not expect the cascade of OIF or FM/W modes to produce waves with the right polarization to resonate with electrons. the velocity shear of the plasma, producing an increase in the internal energy of the particles. In our case, assuming no heat flux, the internal energy density of species j, U j (= p ⊥,j + p ,j /2), evolves as (Kulsrud 1983;Snyder et al. 1997;Sharma et al. 2007): where q = −sB x B y /B 2 corresponds to the growth rate of B. In the present context, both q and ∆p j are negative, which implies an increase in U j . Before the onset of the instabilities, this process is adiabatic and therefore it is a reversible energy gain (in the sense that the increase in U j would be reverted by reversing the direction of the plasma shear velocity). Indeed, as shown in Figures 1a and 1b for the case of electrons, the early evolution of p ⊥,j and p ,j follows the CGL or double adiabatic behavior with p ⊥,j ∝ B and p ,j ∝ 1/B 2 (which gives rise to a net growth of U j = p ⊥,j + p ,j /2 since the p ,j growth occurs faster than the decrease in p ⊥,j ). Thus, only after the instabilities start keeping ∆p j /p ,j in a quasi-stationary regime by breaking µ j invariance (after t · s ≈ 0.7 in Figures 1a and 1b), the increase in U j can be considered as irreversible heating. Also it is important to point out that the role of the instabilities after t · s ≈ 0.7 is not the direct heating of the particles by wave-particle interactions. Instead, the role of the instabilities is to limit the pressure anisotropy and, therefore, to regulate the viscous heating provided by Equation 3. Figure 10 quantifies the importance of this heating mechanism by showing the volume-averaged ion (solidblack) and electron (solid-geen) heating rates for run F2. We also show the heating rate by anisotropic viscosity predicted by Equation 3 for ions (dotted-black) and electrons (dotted-green). 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 T e ∼ T i ), 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).

ELECTRON MEAN FREE PATH
Besides regulating the effective plasma viscosity, pitchangle scattering by velocity-space instabilities is also expected to limit the mean free path of the particles. To quantify this effect we used the distance D j (t) traveled along B by 2×10 4 ions and electrons 13 . We calculated their mean free paths assuming that D 2 j = tv th,j λ j (where λ j represents the average mean free path over species j, while v th,j = (k B T j /m j ) 1/2 is their thermal speed), which is valid if the particles move diffusively. This allows us to estimate the average mean free path of species j as λ j = d D 2 j /dt/v th,j .
where v j is the particle's velocity. Fig. 11.-Electron (black) and ion (red) mean free paths (normalized by v th,j /s), calculated via the time derivative of the mean squared distance traveled by particles along B ( λ j = d D 2 j /dtv th,j ). Results correspond to runs with m i /me = 64 (dotted lines; run F1) and m i /me = 25 (solid lines; run F2). At early times the particles stream freely, d D 2 j /dt ∝ t. After the velocity-space instabilities saturate, pitch-angle scattering leads to saturation of the mean free path. Figure 11 shows our estimates of the electron (black) and ion (red) mean free paths for two simulations with m i /m e = 25 and 64 (runs F2 and F1, respectively), normalized by v th,j /s. We see that in both cases there is an initial period when λ j /(v th,j /s) increases as ∼ 2t · s, which is consistent with the free streaming behavior D j ≈ tv th,j . This is followed by the saturation of λ j , expected to start after a time of the order of the collision time of particles. The behaviors of λ e and λ i in this stage are quite similar, with λ e /(v th,e /s) being somewhat smaller than λ i /(v th,i /s) (by a factor ∼ 1.5).
The evolutions of λ j seen in Figure 11 are expected to be influenced by the pitch-angle scattering of the particles caused by the velocity-space 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 j was incorporated using an effective scattering rate ν ef f,j (Kulsrud 1983;Snyder et al. 1997). In this model, valid for ∆p j /p ||,j ≪ 1, λ j behaves as: 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 λ j 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 ∆p j /p ||,j for runs F2 and F1 (see Figures 7d and 7e) one obtains that, at t · s = 1.5, λ j ≈ 0.15v th,j /s. However, Figure 11 shows that the average mean free paths of ions and electrons are ∼ 10 times larger than this simple estimate. On the other hand, considering the evolution of ∆p j /p ||,j and B 2 /B x B y (needed to determine q) from t · s = 1.5 to t · s = 3, λ i should decrease by a factor ∼ 2, which is essentially what Figure 11 shows. Thus, putting aside the factor ∼ 10 difference, the scaling of λ j on ∆p j /p ||,j and q presented in Equation 4 is reasonably well reproduced by the simulations with decreasing B.
The factor ∼ 10 discrepancy is interesting, partly because the behavior of ν ef f,j suggested by Equation 4 (ν ef f,j ≈ 3qp ||,j /∆p j ) is well reproduced in the case of ions in previous hybrid-PIC simulations that studied the saturated state of the firehose and mirror instabilities (Kunz et al. 2014). In that case, however, ν ef f,i is not measured using D i . Instead, they constructed a distribution of the times taken by each ion to change its µ i by a factor of e, and then approximated ν −1 ef f,i 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 µ e (we will refer to this estimate as λ e,µ ) and the other one using D e ( λ e ). These measurements of λ e and λ e,µ are not defined for different times (as in Figure 11), but they correspond to averages between t · s = 1 and 3. 14 The comparison was made for different groups of electrons, defined by the parameter A e ≡ v 2 ⊥,e / v 2 ,e , where represents the average between t · s = 1 and 3 for each electron. We use A e as a way to quantify the pitch-angle of electrons, which, as we will see below, affects the behaviors of λ e and λ e,µ in different ways.
Our results are shown in Figure 12a. We see that λ e and λ e,µ roughly coincide for A e 1 (pitch-angle 45 o ). However, for A e 1 (pitch-angle 45 o ), 14 For a given electron population, λe is estimated by first calculating for each electron the quantity λe = (D(t · s = 3) 2 − D(t · s = 1) 2 )/( |v ||,e | ∆t) (where ∆t is simply the time elapsed between t · s = 1 and 3 and |v ||,e | is the average magnitude of the electron velocity parallel to B in the same period), and then taking the average over the population. For λe,µ , we construct a distribution of the times taken by each electron to change µe by a factor e with respect to its value at t · s = 1, and then λe,µ is estimated multiplying the width of the distribution by v th,e . ⊥,e are respectively the electron velocity parallel to B and to the z axis, which are mutually perpendicular). Run F1 shows that, for ve ⊥,e ) 1/2 ), the distribution is dominated by electrons with small pitch-angle, while in run MW1 the velocity distribution appears more similar to a bi-Maxwellian distribution with p ⊥,e > p ||,e . λ e becomes significantly larger than λ e,µ . This is consistent with the fact that, for electrons with small pitch-angle, an order unity variation in v 2 ⊥,e due to scattering (which implies an order unity variation in µ e ) does not imply that they reverse their velocity along B. Also, when averaged over the entire A e distribution (shown by the black line in Figure 12a), λ e is a factor ∼ 10 larger than λ e,µ , showing that using λ e,µ would essentially eliminate the discrepancy between our estimated mean free path and Equation 4. 15 Given this difference between λ e,µ and λ e , it is important to understand why in the case of growing magnetic field studied by Riquelme et al. (2016) the behavior of λ e is well reproduced by Equation 4. Figure 12b shows the same quantities as Figure 12a but for run MW1 of Riquelme et al. (2016). We see that the trend for λ e to grow relative to λ e,µ as A e decreases is maintained in this case. However, for all values of A e , the growing B case tends to have a smaller ratio λ e / λ e,µ compared to the decreasing B case, making λ e ∼ λ e,µ if the average over the whole distribution of A e (black line) is considered.
This difference between the growing and decreasing B cases appears to be due to the specific effect of the relevant instabilities on the electron velocities. This is suggested by Figure 13a, which shows the electron velocity distribution f (v ||,e , v (z) ⊥,e ) for run F1 at t · s = 1.5 (corresponding to the saturated stage of the FM/W and OIF instabilities), where v ||,e and v (z) ⊥,e are respectively the electron velocity parallel to B and to the z axis (which are mutually perpendicular). We see that, ⊥,e ) is dominated by electrons with rather small pitch-angle ( 45 o ). This suggests that for v e 0.6c, the scattering process occurs in a way that disfavors the diffusion of electrons towards smaller values of |v ||,e |, which in turn precludes the reversal of v ||,e , contributing to increasing λ e . For comparison, in Figure 13b we show the analogous electron distribution for run MW1 of Riquelme et al. (2016) ⊥,e ) appears more similar to a bi-Maxwellian distribution with p ⊥,e > p ||,e . Notice that the dominance of small pitch-angle electrons for v e 0.6c seen in the decreasing field case is similar to the modification to the ion velocity distribution found by previous hybrid-PIC simulations of an expanding box, where the ion scattering is also provided by the FM/W and OIF instabilities (Matteini et al. 2006;Hellinger et al. 2008).
In summary, both for growing and decreasing B, the estimate of the electron mean free path provided by Equation 4 is fairly well reproduced by λ e,µ (as it was shown by Kunz et al. (2014) in the case of ions). However, since λ e 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. In the decreasing field case, λ e is a factor ∼ 10 larger than λ e,µ , most likely due to the specific electron scattering mechanism provided by the FW/W and OIF instabilities. Thus, in this case the electron mean free path is best described by the relation: which is valid for 2 β e 20, and only differs from Equation 4 by its prefactor ∼ 3 instead of ∼ 0.3.

DISCUSSION AND IMPLICATIONS
We used particle-in-cell (PIC) plasma simulations to study the nonlinear, saturated stage of various ion and electron velocity-space instabilities relevant for collisionless plasmas. We focused on instabilities driven by pressure anisotropy with p ⊥,j < p ||,j . To capture the nonlinear regime in a self-consistent way, we imposed a shear velocity in the plasma, which decreases the background magnetic field. This drives p ⊥,j < p ||,j 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 velocity-space instabilities, which inhibit the growth of pressure anisotropy. The relevant instabilities in this regime, as suggested by linear theory, are: i) the purely growing oblique ion-firehose (OIF) and the resonant fast magnetosonic/whistler (FM/W) modes, which are mainly driven by the ions, and ii) the purely growing oblique electron-firehose (OEF) and the resonant Alfvén/ion-cyclotron (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 m i /m e = 25 and 64. Our results, valid for the regime 2 β i ≈ β e 20, 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 ∆p e ≈ ∆p i . 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 ∆p e = ∆p i and with growth rate ∼ s. 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 ∼ 2 larger anisotropy than in the cases with finite m i /m e . We attribute this result to the rather strong destabilizing effect of the electron pressure anisotropy, ∆p e , on the OIF and FM/W modes (as already suggested by previous linear dispersion analyses Michno et al. (2014); Maneva et al. (2016)), which in turn maintains ∆p e at a value significantly lower than the one necessary to make the OEF and A/IC modes grow at a rate ∼ s.
Although the amplitude of the OIF and FM/W modes depend on the ratio ω c,i /s, the value of the parameters m i /m e and ω c,e /s 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 m i /m e = 1836, highly magnetized (ω c,i /s ≫ 1) 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 pitch-angle scattered in the saturated stage of the OIF and FM/W instabilities in the case of m i /m e = 1836. Answering this question 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 10 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 j, λ j , 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 j is determined by an effective scattering rate ν ef f,j that sets the rate at which the µ j invariance is broken (Kulsrud 1983;Snyder et al. 1997). However, the prefactor in Equation 5 is ∼ 10 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 µ e provides a good estimate of λ e only for electrons with relatively large pitch-angle ( 45 o ). For electrons with pitch-angles smaller than ∼ 45 o , λ e can become a factor ∼ 10 larger. This is in contrast with the case of growing B (Riquelme et al. 2016), where Equation 4 provides a reasonably good estimate for λ j , despite the fact that λ j also grows as the pitch-angle decreases. This difference is likely due to the specific electron scattering mechanism provided by the FW/W and OIF instabilities, which tends to preclude the reversal of v ||,e , contributing to increase λ e .
The results shown in this work, as well as those presented in Riquelme et al. (2016) for the growing B case, are relevant for quantifying the viscous heating and thermal conductivity in various low-collisionality astrophysical plasmas, including low-luminosity accretion flows around compact objects, the ICM, and the heliosphere. This work was supported by NSF grants AST 13-33612 and 1715054, a Simons Investigator Award to EQ from the Simons Foundation and the David and Lucile Packard Foundation.
DV also acknowledges support from NSF/SHINE grant AGS-1460190, NSF/SHINE grant AGS-1622498, NASA grant NNX16AG81G, and a UK STFC Ernest Rutherford Fellowship (ST/P003826/1). We are also grateful to the UC Berkeley-Chile Fund for support for collaborative trips that enabled this work. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575.