<figure><img src="image_1.png" /></figure>

HAL
open science

# Interplay between the non-resonant streaming instability and self-generated pressure anisotropies

Alexis Marret, Andrea Ciardi, Roch Smets

► To cite this version:

Alexis Marret, Andrea Ciardi, Roch Smets. Interplay between the non-resonant streaming instability and self-generated pressure anisotropies. Monthly Notices of the Royal Astronomical Society, 2024, 532 (4), pp.4082-4088. 10.1093/mnras/stae1773. hal-04679927

HAL Id: hal-04679927

https://hal.science/hal-04679927v1

Submitted on 30 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Interplay between the non-resonant streaming instability and self-generated pressure anisotropies

A. Marret<sup>1,2,3*</sup> A. Ciardi<sup>2</sup> and R. Smets<sup>3</sup>

Accepted 2024 July 19. Received 2024 July 3; in original form 2024 January 25

## ABSTRACT

The non-thermal particles escaping from collision-less shocks into the surrounding medium can trigger a non-resonant streaming instability that converts parts of their drift kinetic energy into large amplitude magnetic field perturbations, and promote the confinement and acceleration of high energy cosmic rays. We present simulations of the instability using an hybrid-Particle-in-Cell approach including Monte Carlo collisions, and demonstrate that the development of the non-resonant mode is associated with important ion pressure anisotropies in the background plasma. Depending on the initial conditions, the anisotropies may act on the instability by lowering its growth and trigger secondary micro-instabilities. Introducing collisions with neutrals yield a strong reduction of the magnetic field amplification as predicted by linear fluid theory. In contrast, Coulomb collisions in fully ionized plasmas are found to mitigate the self-generated pressure anisotropies and promote the growth of the magnetic field.

**Key words:** instabilities – magnetic fields – plasmas.

## 1 INTRODUCTION

A population of ions drifting at super-Alfvénic speeds with respect to a magnetized ambient plasma can trigger various streaming instabilities, leading to the exponential growth of electromagnetic perturbations across a wide range of astrophysical environments (Kulsrud & Pearce 1969; Gary et al. 1984; Winske & Leroy 1984; Bell 2004; Amato & Blasi 2009; Amato 2014; Cui, Pühlhofer & Santangelo 2016). Depending on the plasma conditions, a mode not involving particle-wave resonances can grow (Winske & Leroy 1984) and drive turbulent magnetic field amplification that can greatly exceed the initial seed field (Bell 2004). The significance of this non-resonant streaming instability (NR) was initially recognized for its role in the generation of diffuse ion distributions in the Earth's foreshock region, where reflected ions streaming back from Earth's bow shock encounter the incoming solar wind ions (Sentman, Edmiston & Frank 1981; Onsager, Winske & Thomsen 1991; Akimoto et al. 1993). The discovery of its ability to strongly amplify the magnetic field has now made it a key ingredient in the diffusive shock acceleration of cosmic rays in supernovae remnants of particles up to PeV energies (Bell 2004; Amato & Blasi 2009; Bell 2013; Marcowith et al. 2016). The NR instability has been studied numerically using a variety of computational techniques, including modified magnetohydrodynamics (MHD) (Bell 2004; Zirakashvili, Ptuskin & Völk 2008), hybrid-particle-in-cell (PIC ions and massless fluid electrons) (Winske & Leroy 1984; Akimoto et al. 1993; Haggerty & Caprioli 2019; Marret et al. 2021), full-PIC (Ohira et al. 2009; Riquelme & Spitkovsky 2009; Crumley et al. 2019),

and MHD-PIC (Bai et al. 2015; Mignone et al. 2018; van Marle, Casse & Marcowith 2018) simulations.

The MHD-PIC method has received growing attention as it combines the kinetic treatment of the cosmic rays while retaining the advantage of modelling the background plasma as a magnetofluid, over large spatial and temporal scales. Neglecting kinetic effects in the background plasma however is not always justified. For example, in the hot plasmas of superbubbles or in the intergalactic medium, the background's ions thermal Larmor gyroradius can become comparable to or larger than the unstable wavelengths. Under these conditions a kinetic treatment of the background ions population is necessary to account for resonance effects that can significantly reduce the growth of the NR instability (Reville et al. 2008; Zweibel & Everett 2010; Marret et al. 2021). In addition, the development of significant ion pressure anisotropies in the background plasma have been observed in collision-less hybrid-PIC simulations (Marret et al. 2021), and is found to be well described within the CGL adiabatic theory (Marret et al. 2022). This suggests that the assumption of an isotropic scalar pressure, often employed in fluid models, may not be well-suited. Pressure anisotropies can be suppressed by particle collisions or by micro-instabilities such as the mirror and ion-cyclotron mode (Gary et al. 1976), among other isotropization mechanisms. In this work we expand on previous studies (Marret et al. 2022) and investigate the interplay between magnetic field amplification and self-generated pressure anisotropies in plasmas with varying levels of collisionality.

The NR instability can be described by considering a modified MHD model (Bell 2004) with two components: a background plasma population, and a less dense cosmic ray population of protons, drifting with a velocity $u_{cr}$ parallel to an ambient magnetic field $B_0$. The perpendicular case is not considered here because of its slower growth rate (Matthews et al. 2017). The background plasma is electrically

<footnote n="1"><sup>1</sup>High Energy Density Science Division, SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="2"><sup>2</sup>Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, LERMA, F-75005 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Sorbonne Université, Ecole Polytechnique, CNRS, Observatoire de Paris, LPP, F-75005 Paris, France</footnote>

<footnote>*E-mail: alexis.marret@slac.stanford.edu</footnote>
charged to satisfy the quasi-neutrality condition $n_m - n_e = -n_{cr}$, where the subscripts $m$, $e$, and $cr$ denote the main protons, electrons, and proton cosmic rays, respectively. The background plasma carries a return current, produced by the electron drift relative to the main protons as $u_e = u_{cr}n_{cr}/n_m$, and which initially compensates the cosmic rays current. The electric field is obtained from an ideal Ohm's law with an additional term due to the presence of the low density drifting cosmic rays (Bai et al. 2015; van Marle et al. 2018): $E = -u \times B - j_{cr} \times B/\epsilon m$, where $j_{cr} = \epsilon m u_{cr}$ is the current carried by the cosmic rays, initially opposed to the background return current, $e$ is the elementary charge, and $u$ is the background fluid velocity. The essence of the instability can then be captured (Marret et al. 2021) by considering the simplified linearized momentum conservation for the background MHD fluid and Maxwell-Faraday's equations

$$ \rho_0 \frac{\partial u_1}{\partial t} = -j_{cr} \times B_1 \quad (1) $$

$$ \frac{\partial B_1}{\partial t} = (B_0 \cdot \nabla) u_1, \quad (2) $$

where $\rho$ is the main plasma mass density. The numerical subscripts refer to initial and perturbed quantities. We are interested in electromagnetic modes and neglect density fluctuations. In addition, the cosmic rays current is taken to be constant and the magnetic tension is neglected. These approximations remain valid (Bell 2004) for unstable wave numbers in the range $k_{min} < k < k_{max}$, where

$$ k_{min} = \Omega_0 / u_{cr} \quad (3) $$

$$ k_{max} = \frac{n_{cr}}{n_m} \frac{u_{cr}}{v_{A0}} \Omega_0 \quad (4) $$

$\Omega_0 = eB_0/m_p$ is the proton gyro-frequency, $m_p$ the proton mass, $v_{A0} = B_0/\sqrt{\mu_0\rho_0}$ the Alfvén velocity, and $\mu_0$ the vacuum permeability. Larger scale perturbations with $k \le k_{min}$ are unstable to resonant modes, which require a fully kinetic treatment of the cosmic rays population (Amato & Blasi 2009). The coupled equations (1) and (2) show that the NR mode is driven by the cosmic ray current through the action of the magnetic force, $-j_{cr} \times B_1$, which produces fluid velocity fluctuations, $u_1$, in the background plasma. The induced electric field, $-u_1 \times B_0$, feeds back and enhances the initial magnetic field perturbation via Faraday's law, promoting its exponential growth. Considering a right-hand circularly polarized parallel propagating magnetic field perturbations of the form $B_1 = \hat{B}_1 e^{i(kx-\omega t)}$, where the angular frequency $\omega = \omega_r + i\gamma$, $\omega_r \ll \gamma$ is taken to be positive, the growth rate for negative $k$ is given by retaining the magnetic tension force in equation (1) that dominates at the smallest scales yields a cutoff of the instability at $k = k_{max}$, and a fastest growing mode with $\gamma_{nr} = k_{max}v_{A0}/2$ for $k = k_{max}/2$. The magnetic field amplification can lead to pressure anisotropies due to adiabatic invariants conservation, which can in turn produce pressure gradients in the background plasma (Marret et al. 2022). The effect of pressure gradients on the instability cannot be captured by linear theory, since no transverse contribution remains in the linearized momentum conservation equation in the gyrotropic limit (Bell 2004). We will show that non-linear pressure gradients effects can none the less modify the growth of the NR mode by opposing the magnetic force driving the magnetic field amplification.

## 2 NUMERICAL SETUP

We present the results of 1D and 2D simulations performed with the hybrid-PIC-MCC code HECKLE (Smets & Aunai 2011), which solves the Vlasov-Maxwell system including Monte Carlo collisions,

using a predictor-corrector scheme for the electromagnetic field, a non-relativistic Boris pusher (Boris 1970), and a first order interpolation scheme for current deposition on the grid. The main protons and cosmic rays protons are described with macroparticles following the PIC method, and the electrons are treated as a mass-less neutralizing fluid. This hybrid approach is well suited to study the kinetic, non-linear evolution of systems at the ions temporal and spatial scales while avoiding prohibitive computational time. The density and magnetic field are normalized to their initial uniform values $n_0 = n_m(t_0)$ and $B_0 = B(t_0)$. Times and lengths are normalized to the initial inverse proton cyclotron angular frequency $\Omega_0^{-1}$ and proton inertial length $l_0 = v_{A0}\Omega_0^{-1}$.

We consider a population of cold main protons and electrons, traversed by a population of super-Alfvénic cosmic rays (protons) with a density $n_{cr}/n_m = 0.01$. The initial drift velocity, $u_{cr}/v_{A0} = 50$, is oriented parallel to the initial magnetic field $B_0 = B_0e_x$. The cosmic rays current is not maintained externally, and thus decreases strongly upon the NR mode growth and saturation. The plasma and field quantities are initially homogeneous, and periodic boundary conditions were used in all directions. The simulation domain has a length $L_x = 10^3l_0$ where $l_0$ is the proton inertial length, discretized with $10^3$ cells for 1D simulations, and extended to $L_y = 2 \times 10^2l_0$ in the 2D case with the same resolution as in the $x$ direction. We used 200 macroparticles per cells in 1D simulations, 100 for each proton species. We used 600 macroparticles per cell for 2D simulations, 500 for the background protons, and 100 for the cosmic rays. These values have been chosen to mitigate large numerical electric fields that may develop in the hybrid-PIC model in regions with few macroparticles per cell, while achieving numerical convergence and minimizing the numerical cost of the simulations. The electrons pressure is assumed to follow an isothermal behaviour $P_e(t) = n_e(t)T_e(t_0)$ with $T_e(t_0) = T_m(t_0) = m_p v_{A0}^2$. Unless stated otherwise, the initial plasma-$\beta$ is $\beta_0 = P_0/W_{B0} = 2$ where $P_0$ is the initial, isotropic main protons pressure and $W_{B0} = B_0^2/2\mu_0$ is the initial magnetic field energy density. We also investigated the case $\beta_0 = 10$ by increasing the main protons initial temperature. The initial main protons Larmor radius is then resolved with two cells for $\beta_0 = 10$ and 1 cell for $\beta_0 = 2$.

## 3 PRESSURE ANISOTROPIES GENERATION AND INTERPLAY WITH MICROINSTABILITIES

To emphasize the relationship between the development of the NR mode and the emergence of pressure anisotropies, we display in Fig. 1 the evolution of the magnetic field energy density ($W_B = B^2/2\mu_0$), together with the perpendicular and parallel components of the main protons' pressure tensor $P_m^{\perp}/P_m^{\parallel}$ with respect to the magnetic field, computed as the second-order moment of the velocity distribution. Between $t = 15$ and $t = 34 \Omega_0^{-1}$, the magnetic perturbations exhibit an exponential increase, accompanied by significant anisotropies, with a maximum spatially averaged value of $P_m^{\perp}/P_m^{\parallel} = 36.7$. Upon reaching saturation, the magnetic field energy density can be predicted by considering the energy exchange rates obtained within quasi-linear theory (Winske & Leroy 1984), which indicates that the rate of energy gained by the magnetic field is half that of the cosmic rays' drift kinetic energy loss. Extrapolating this result to the non-linear evolution and supposing that the initial cosmic rays drift kinetic energy density is entirely depleted at saturation, one obtains $W_{B,sat} = W_{cr}/2$ with $W_{cr} = n_{cr}m_p u_{cr}^2/2$. This estimate is close to the simulations results $W_{B,sat} = 0.55W_{cr}$ obtained by averaging the magnetic field energy density over space, consistent with the fact that
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Evolution of the spatially averaged magnetic field energy density $W_B$ (solid black line) and of the spatial average of the $P_m^{\perp}$ and $P_m^{\parallel}$ components of the main protons pressure tensor calculated from the local macroparticles distribution (solid lines), in a 1D collisionless simulation. The dotted lines show the CGL predicted pressures $P_m^{\perp} = P_0 \rho B / \rho_0 B_0$ and $P_m^{\parallel} = P_0 \rho^3 B_0^2 / \rho_0^2 B^2$, calculated using the density and magnetic field extracted from the simulations.</figcaption></figure>

the cosmic rays current is not externally maintained and decreases strongly over the course of the simulation.

The development of pressure anisotropies can be described within the adiabatic CGL theory (Chew et al. 1956), which may be interpreted as the conservation of the first and second adiabatic invariants in a slowly varying magnetic field (Kulsrud 1983) and predicts a decrease in parallel pressure and an increase in perpendicular pressure in regions of amplified magnetic field for a constant density. Taking the electrons to be magnetized, considering an ideal Ohm's law (without the cosmic rays contribution), and neglecting heat fluxes and non-gyrotropic components of the pressure tensor, the main protons pressure parallel and perpendicular to the magnetic field may be expressed as

$$ \frac{d}{dt} \left( \frac{P_m^{\parallel} B^2}{\rho^3} \right) = 0 \qquad (5) $$

$$ \frac{d}{dt} \left( \frac{P_m^{\perp}}{\rho B} \right) = 0, \qquad (6) $$

where $d/dt = \partial_t + u_m \cdot \nabla$ denotes the material derivative. The advective part of the derivative may be neglected by integrating over the periodic simulation domain, and by neglecting the remaining velocity divergence term, which is a well verified assumption in the simulations. In this case the pressure components may be evaluated directly as $P_m^{\perp}/P_0 = \rho B/\rho_0 B_0$ and $P_m^{\parallel}/P_0 = \rho^3 B_0^2/\rho_0^2 B^2$. These equations yield the well-known effect that pressure anisotropies $P_m^{\perp}/P_m^{\parallel} > 1$ can be produced in an amplified magnetic field at constant density, and are a good model for the simulations. In Fig. 1, we display the evolution of the pressure components $P_m^{\perp}$ and $P_m^{\parallel}$ predicted from the CGL model and calculated using the density and magnetic field from the simulation. The pressure anisotropies driven by the NR mode evolve according to the CGL equations up to saturation ($t \approx 35 \Omega_0^{-1}$). This result confirms that the observed pressure anisotropies are not a byproduct of the restricted 1D and 2D geometry of the simulations. We note that the CGL equations do not account for the cosmic rays contribution in Ohm's law, but still describe well the anisotropies.

Pressure anisotropies $P_m^{\perp}/P_m^{\parallel} > 1$ can play an important role in the plasma dynamics by triggering the growth of the ion-cyclotron (ic) and mirror (mi) microinstabilities (Gary et al. 1976) upon reaching the threshold $P_m^{\perp}/P_m^{\parallel} - 1 \ge S_p/(\beta_m^{\parallel} - \beta_p)^{\alpha_p}$ (Hellinger et al. 2006), where $S_p$, $\beta_p$, and $\alpha_p$ are fitting parameters determined numerically for a given $\gamma_{mi}$ and $\gamma_{ci}$. The microinstabilities growth rate increases with larger $P_m^{\perp}/P_m^{\parallel}$, and is stabilized for smaller $\beta_m^{\parallel}$. While both

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Distribution (simulation cell count) of the ratio $P_m^{\perp}/P_m^{\parallel}$ as a function of $\beta_m^{\parallel} = P_m^{\parallel}/W_B$, obtained in 1D collisionless simulations with $\beta_0 = 2$ (left panels) and $\beta_0 = 10$ (right panels), at different times during the growth of the NR instability. The solid red line shows the anisotropy expected from incompressible CGL theory $P_m^{\perp}/P_m^{\parallel} = (\beta_m^{\parallel}/\beta_0)^{-3/4}$. The solid and dashed black lines indicate the thresholds for the mirror $\gamma_{mi} = 10^{-1}\Omega_0$ and ion-cyclotron $\gamma_{ci} = 10^{-1}\Omega_0$ modes respectively.</figcaption></figure>

instabilities feed on pressure anisotropies to develop, the mirror mode can be described at first order within a MHD framework, whereas the ion-cyclotron mode requires a kinetic description to correctly describe the wave-particle resonance leading to its growth (Southwood & Kivelson 1993). To understand the role of microinstabilities on the growth of the NR mode, we display in Fig. 2 the distribution (cell count) of the ratio $P_m^{\perp}/P_m^{\parallel}$ as a function of
$β_m^‖ = P_m^‖/W_B$, for two collision-less simulations with $β_0 = 2$ and $β_0 = 10$. Fig. 2 also shows the threshold anisotropy for the ion-cyclotron and mirror instabilities, obtained from linear kinetic theory assuming a homogeneous plasma with bi-Maxwellian populations (Hellinger et al. 2006) and for a growth rate comparable to that of the collision-less NR mode ($γ_0 = 0.15Ω_0$). Using the approximated forms of equations (5) and (6), the anisotropy within the CGL theory can be described with a power law as

$$ \frac{P_m^\perp}{P_m^\parallel} = \left( \frac{\beta_m^\parallel}{\beta_0} \right)^{-3/4} \qquad (7) $$

which is well recovered in the simulations with $β_0 = 2$ (Fig. 2 left panels), indicating that the anisotropies are not constrained by microinstabilities. In particular the mirror mode growth rate $γ_{mi} \sim 10^{-2}Ω_0$ always remains below that of the NR mode measured in the simulation $γ_{nr} = 0.15Ω_0$.

A different picture emerges for the simulation with $β_0 = 10$, obtained by increasing the main protons initial temperature, and shown in the right panels of Fig. 2. In this case the growth of the mirror mode competes with the NR mode ($γ_{mi} \sim γ$) for $P_m^\perp/P_m^\parallel \gtrsim 2$ and keeps the pressure anisotropies at the threshold values, below the expected power law. Similar results were obtained in 2D simulations. Although in both low and high-$β$ cases the ion-cyclotron instability has nominally a larger predicted growth rate than the mirror mode, and comparable or larger to that of the NR mode for $P_m^\perp/P_m^\parallel \gtrsim 3$, it does not appear to be limiting the anisotropy; a similar behaviour has also been observed in the solar wind where pressure anisotropies above the expected ion-cyclotron threshold are frequently measured (Hellinger et al. 2006). In our case it may be due to the spatial inhomogeneity of the magnetic field amplification (see Section 4), which renders the cyclotron resonance position dependent and impairs the growth of the ion-cyclotron mode (Southwood & Kivelson 1993). Studies of the growth of microinstabilities in the context of amplified magnetic field have also shown that the ion-cyclotron mode subdominance can be attributed to local departures from the Maxwellian distribution function (Riquelme, Quataert & Verscharen 2015; Walters et al. 2023) that is assumed when computing the ion-cyclotron mode linear properties, however we did not observe such departure in the simulations.

## 4 COLLISIONAL NON-RESONANT MODE

In a collisional plasma, pressure anisotropies may be mitigated if collisions are sufficiently frequent to redistribute the energy in all directions of space. Expanding on previous works (Marret et al. 2022), we investigate the two cases of first, a poorly ionized background where collisions of ions with a population of neutrals are dominant, relevant for example in the warm neutral medium (Recchia et al. 2021), and secondly, a fully ionized background plasma where ion Coulomb collisions are dominant, relevant to the interstellar medium. While intraspecies Coulomb collisions conserve the momentum and energy (Trubnikov 1965) of the population of particles of the same species, proton-neutral collisions do not, as protons loose their energy to the denser and colder background neutral population. This leads to important damping of the electromagnetic waves, including the NR mode (Forteza et al. 2007; Reville et al. 2007).

The proton-proton Coulomb collisions are implemented numerically using a Monte Carlo method which solves the Landau collisions operator by randomly pairing macroparticles in each cells, and calculating at each time-step the associated scattering angle and post-collision velocities (Takizuka & Abe 1977). The proton-

neutral collisions are implemented following a hard-sphere model (Nanbu 2000). We considered proton-Hydrogen elastic collisions, and a small ionization fraction such that the neutrals density and temperature are supposed to remain constant and uniform. At each time-step a collision between a macroparticle and an hydrogen atom occurs for $r < n_nσ_{in}ΔvΔt$, where $r \in [0, 1]$ is a randomly generated number, $n_n$ is the neutral density, $σ_{in}$ the ion-neutral collision cross-section, $Δv$ the relative velocity, and $Δt$ the time-step. The cross-section collision energy dependence is obtained from Krstic & Schultz (1999) assuming a neutral population temperature $T_n = 1$ eV as typically observed in the warm atomic phase of the interstellar medium (Ferriere 2001). The Coulomb and neutral collision frequencies, noted $v_0$, can be scaled with respect to the instability growth rate $γ$, allowing us to probe the weakly ($γ_{nr} > v_0$) and strongly ($γ_{nr} < v_0$) collisional regimes of the NR mode in weakly or fully ionized plasmas. The electron resistivity is negligible for the range of Coulomb collision frequencies investigated from $v_0 = 10^{-2}$ to $v_0 = 10^2 Ω_0$. Larger values were not investigated because of the prohibitively small numerical time-steps required to resolve low energy collisions.

Fig. 3 presents the magnetic field growth rate $γ$ as a function of the reference collision frequency $v_0$, for simulations including (a) neutral collisions and (b) Coulomb collisions. We find that the neutral collisions reduce the background fluid velocity perturbations, leading to a weaker induced electric field and consequently smaller growth rate and magnetic field amplification. The damping of the instability by proton-neutral collisions can be studied by introducing a collisional drag term in the background fluid momentum conservation equation of the form $-ρv_0(u - u_n)$ where $u_n$ is the neutral fluid velocity. The growth rate in the case of a weakly ionized plasma is then given by $γ_{in}(k) = -\frac{v_0}{2} + \frac{1}{2}(v_0^2 + 4Ω_0 \frac{n_{cr}}{n_m} ku_{cr})^{1/2}$ (Reville et al. 2007). This is plotted in Fig. 3(a) considering $k = k_{max}/4$ such that $γ_{in}(v_0 = 0) = γ_{nr}$, where $γ_{nr} = k_{max}v_{A0}/2$. The growth rate dependence with collision frequency is well recovered in the simulations (red triangles). The magnetic field intensity has been integrated over the $k$ spectrum before measuring the growth rate, in order to reduce fluctuations due to the dynamic nature of the range of unstable wavenumbers inherent to the NR mode (Bell 2004; Marret et al. 2021). This gives an overall smaller growth rate than if only the fastest growing mode was observed and is seen with the offset by a factor 2 in the figure.

The simulations of a fully ionized collisional background [Fig. 3(b) black dots] give an enhanced growth rate compared to the collisionless case $γ_0 = 0.15Ω_0$ for $v_0 > γ_0$. The increase is maximum for a collision frequency $v_0 = 27Ω_0$ two orders of magnitude larger than $γ_0$, yielding a growth rate $γ = 0.17Ω_0$. This case will be referred as the collisional simulation in the following. The saturated magnetic field energy density $W_{B,sat}$ is displayed in Fig. 3(c), and shows an increase up to $W_{B,sat} = 8.7l_0^{-3}m_p v_{A0}^2$ corresponding to $W_{B,sat}/W_{B,sat}^{cl} = 1.3$ with $W_{B,sat}^{cl}$ the saturated magnetic field energy density in the collision-less case. The fastest growing wavenumber is not modified by the collision frequency, with unstable waves growing on scales of the order $λ \approx 4πk_{max}^{-1} \approx 25l_0$ in agreement with the linear kinetic theory prediction for a negligible background plasma temperature (Winske & Leroy 1984). Because of the relatively large density of cosmic rays in the simulations, their collisions with the background protons become frequent for $v_0/Ω_0 \gg 10$, leading to a rapid loss of their drift kinetic energy and to a lower magnetic field amplification. This was verified by performing simulations where cosmic rays collisions were artificially suppressed (orange squares). In the case of completely collision-less cosmic rays, which is also more representative of the conditions found in space, the growth
<figure><img src="image_5.png" /></figure>

**Figure 3.** a) Measured magnetic field growth rate $\gamma$ ($C_{m/n}$, triangles) as a function of the collision frequency, in 1D runs including proton-neutral collisions where $v_0 = n_n\sigma_n v_{T0}$ is the main proton-neutral collision frequency, $n_n$ is the neutral density, and $\sigma_n$ is the neutral collision cross-section. The solid black line corresponds to the theoretical growth rate $\gamma_{in}/2$ considering $k = k_{max}/4$. The collision frequency for which $v_0/\gamma_0 = 1$, where $\gamma_0 = 0.15\Omega_0$ is the growth rate in the collisionless case, is indicated with the vertical dashed line and reported in the other figures. b) Magnetic field growth rate in 1D runs including Coulomb collisions between all protons populations ($C_{m/m;cr/m}$, dots) and including Coulomb collisions between main protons only ($C_{m/m}$, squares), where $v_0 = e^4 n_m \ln \Lambda / (4\pi m_\eta^2 \epsilon_0^2 v_{T0}^3)$ is the collision frequency among the main proton with $\ln \Lambda$ the Coulomb logarithm, and $v_{T0} = (T_m(t_0)/m_p)^{1/2}$ the thermal velocity. The maximum growth rate for simulations with Coulomb collisions is indicated with the second vertical dashed line at $v_0 = 27\Omega_0$, and reported in the other figures. The dashed black line (bottom right) corresponds to the initial cosmic ray-main proton collision frequency $v_{cr}$, normalized to $\gamma_0$. c) Magnetic field energy density $W_{B,sat} = B_{sat}^2/2\mu_0$ at saturation, normalized to the value in collisionless simulations $W_{B,sat}^{cl} = 6.8 l_0^{-3} m_p v_{A0}^2$, and d) mean value of the ratio $P_m^{\perp}/P_m^{\|}$ averaged over the exponential phase of growth, for simulations including Coulomb collisions.

rate and magnetic field energy at saturation remain at the same level. Fig. 3(d) shows the ratio $P_m^{\perp}/P_m^{\|}$ as a function of the collision frequency $v_0$. We find that the observed increase in the amplification of the magnetic field with $v_0$ is correlated to the gradual suppression of the pressure anisotropies for $v_0 \gtrsim 0.1\Omega_0 \sim \gamma_0$.

We interpret this result as follows. The pressure anisotropies generate spatial gradients of the pressure tensor along $B_0$, due to the helical spatial structure of the electromagnetic wave. Those pressure gradients affect the background plasma dynamics in the plane perpendicular to $B_0$ and oppose the magnetic force driving the NR mode. Suppressing these anisotropies promotes the growth of the magnetic field, however, we do not expect that a strong enhancement (> 100 per cent) is possible since it would require the pressure gradients to overcome the cosmic rays magnetic force, which would prevent the growth of the instability altogether. The competition between the magnetic and the pressure gradient forces is illustrated in Fig. 4 for 1D collision-less and collisional simulations. We compute the main protons fluid velocity component $u_m^\times$ in the direction of the local $-j_{cr} \times B$ force, that is opposed by the pressure tensor gradients of the background protons $-\nabla \cdot P_m$. We

<figure><img src="image_6.png" /></figure>

**Figure 4.** Perturbed magnetic field intensity $B_1$ in units of the initial magnetic field $B_0$ (blue) and background ion fluid velocity in the direction of the local $-j_{cr} \times B$ force $u_m^\times$ in units of $v_{A0}$ (red), for 1D collision-less (dashed line) and collisional (solid line) simulations.

find $|j_{cr} \times B|/|\nabla \cdot P_m| \sim 3$ during the exponential phase of growth in the collision-less case, which leads to a less efficient acceleration of the background fluid and to a proportionally smaller magnetic field amplification (equation 2) compared to the collisional case. We verified that collisional viscous forces are negligible in the simulations. We note that in both cases the instability saturation occurs when the fluid velocity $u_m^\times$ changes sign, corresponding to the reversal of the longitudinal electric field $E = -u_m^\times \times B$ that slows down the cosmic rays. This effectively halts the conversion of drift kinetic energy into magnetic field energy and saturates the NR mode (Marret et al. 2021).

In 2D simulations, we find a growth rate of the magnetic field intensity averaged over space marginally larger ($\sim 1$ per cent) in the collisional case. The magnetic field is amplified similarly to the 1D runs, with an average magnetic field energy density ratio between the collisional and collisionless simulations of $W_{B,sat}/W_{B,sat}^{cl} = 1.3$. This is illustrated in Fig. 5, which displays maps of the perturbed magnetic field intensity, together with maps of the pressure anisotropies $P_m^{\perp}/P_m^{\|}$, for 2D collision-less and collisional ($v_0 = 27\Omega_0$) runs. We find that regions of magnetic field amplification are well correlated with the regions of large pressure anisotropies, further confirming our interpretation of the 1D simulations results. In the collisional case, the anisotropies are suppressed ($P_m^{\perp}/P_m^{\|} \approx 1$) and the magnetic field is further amplified compared to the collision-less case. We note that the magnetic field amplification is spatially inhomogeneous, which may explain the absence of the ion-cyclotron mode in the simulations as mentioned in Section 3.

## 5 DISCUSSION

We have presented 1D and 2D hybrid-PIC simulations of the non-resonant streaming instability including Coulomb collisions between protons and elastic collisions with neutrals. We have shown that the instability leads to an important anisotropic heating of the background plasma, well described by the adiabatic CGL theory and, depending on the initial plasma-$\beta$, constrained by the mirror microinstability. The pressure gradients that are generated partially oppose the magnetic force and therefore the growth of the unstable waves. In fully ionized collisional plasmas, proton-proton Coulomb collisions reduce the pressure anisotropies, increasing the growth rate of the unstable waves and the amplification of the magnetic field. Simulations of poorly ionized plasmas, where proton-neutral collisions dominate, instead confirm the strong damping of the NR mode predicted by linear theory calculations (Reville et al. 2007).

One may evaluate the Coulomb collision frequency required to mitigate the pressure anisotropies by comparing the NR mode anisotropic heating rate to the isotropization rate by collisions. The
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Maps of the perturbed magnetic field intensity $B_1$ and of the pressure anisotropy $P_m^\perp/P_m^\parallel$ spatial distribution, taken during the exponential phase of growth at $t = 33\Omega_0^{-1}$ in 2D collision-less ($v_0 = 0\Omega_0$, two upper panels) and collisional ($v_0 = 27\Omega_0$, two lower panels) simulations.</figcaption></figure>

evolution of the anisotropy, $P_m^\perp - P_m^\parallel$, within incompressible CGL theory and including collisions can be expressed as

$$ \frac{\partial}{\partial t} (P_m^\perp - P_m^\parallel) = \gamma (P_m^\perp + 2P_m^\parallel) - v_0 \kappa P_0^{3/2} \frac{1}{\sqrt{P_m^\parallel}}, \quad (8) $$

where the first term on the right hand side is the anisotropic heating rate due to the amplification of the magnetic field (Hunana et al. 2019) by the NR mode, and the second term is the competing pressure isotropization rate due to Coulomb collisions (Trubnikov 1965). The magnetic field growth rate is defined as $\gamma \equiv (dB/dt)/B$, with $B = |B|$; $v_0$ is the fundamental Coulomb collision frequency, and

$$ \kappa = 3 (2\pi^{1/2} A)^{-1} (-3 + (A + 3) [\tan^{-1}(A^{1/2}) / A^{1/2}]) \quad (9) $$

is a decreasing function of $A = P_m^\perp/P_m^\parallel - 1$. Inserting the simplified equations (5) and (6) in equation (8), one then obtains the level of collisionality necessary to stop the growth of the pressure anisotropies and reach steady state

$$ \frac{v_0}{\gamma} = \frac{1}{\kappa}. \quad (10) $$

Larger collisions frequencies, i.e. $v_0/\gamma > \kappa^{-1}$, will start to strongly reduce the development of pressure anisotropies. The function $\kappa$ requires to calculate the parameter $A$, which may be inferred from the saturated magnetic field energy density prediction $W_B \approx W_{cr}/2$ obtained from quasi-linear theory (Winske & Leroy 1984), such that $A = (W_{cr}/2W_{B0})^{3/2} - 1$. The above estimates is valid if the anisotropy is controlled only by the magnetic field amplification. However, pressure anisotropies with $P_m^\perp/P_m^\parallel > 1$ can also drive the growth of the mirror mode. In the case where it competes with the NR mode, the maximum level of anisotropy will then be determined by the mirror mode threshold pressure anisotropy as presented in Section

3. The rapid growth of the mirror mode also leads to important density fluctuations (Southwood & Kivelson 1993) which would invalidate the incompressibility assumption in equation (8), and thus make analytical estimates of $P_m^\perp$ and $P_m^\parallel$ near saturation unreliable, especially in environments with $\beta_0 \gg 1$ where the mirror mode growth rate is maximum (Gary et al. 1997). For this reason, we restrict our analysis to regimes where $\beta_0 \lesssim 1$, such that the parallel plasma beta in the amplified magnetic field is small $\beta_m^\parallel \ll 1$, and the mirror mode remains subdominant.

For our simulation parameters, $\kappa^{-1} = 7.3$, which agrees well with the range of collision frequencies, $v_0/\gamma > 1$, for which pressure anisotropies are seen to be strongly reduced. Keeping the anisotropy small, say $P_m^\perp/P_m^\parallel - 1 < 0.1$, requires collision frequencies $v_0/\gamma > 10^2$, which is again consistent with the values obtained in the simulations. The parameter $\kappa^{-1}$ can also be used to assess the importance of Coulomb collisions and of pressure gradients in various environments. Considering the situation of a supernova remnant (SNR) shock propagating at a velocity $u_{cr} = 5 \times 10^3$ km s$^{-1}$ in a fully ionized interstellar medium with $n_m = 1$ cm$^{-3}$, $B = 5$ µG, $T_m = 10^4$ K, and a cosmic rays flux $n_{cr}u_{cr} = 5 \times 10^4$ cm$^{-2}$ s$^{-1}$ (Zweibel & Everett 2010), we obtain $\beta_0 = 1.4$, $\gamma/\Omega_0 = 2.3 \times 10^{-2}$, $A = 33.0$, and $\kappa^{-1} = 6.8$, which is larger than the normalized proton Coulomb collision frequency (Trubnikov 1965) $v_0/\gamma = 6.2 \times 10^{-3}$. Under such conditions, pressure anisotropies will develop unimpeded and act to reduce the NR mode growth rate and saturated magnetic field. Proton-neutral collisions are not relevant in this case, however one may also consider the situation of a SNR propagating in a molecular cloud. One then finds that the NR mode is strongly damped by the collisions with the weakly ionized background neutrals, preventing magnetic field amplification by the instability (Reville et al. 2007).

In the context of laboratory experiments, and considering parameters typical of high power laser-plasma interactions $u_{cr} = 10^3$ km s$^{-1}$, $n_m = 10^{19}$cm$^{-3}$, $B = 0.2$MG, and $T_m = 10^6$K, and a proton flux $n_{cr}u_{cr} = 10^{26}$cm$^{-2}$s$^{-1}$, one finds $\beta_0 = 8.7 \times 10^{-1}$, $\gamma/\Omega_0 = 3.6 \times 10^{-1}$, $A = 3.3$, and $\kappa^{-1} = 5.6$, smaller than the collision frequency $v_0/\gamma = 29.7$. In this case pressure anisotropies will be mitigated, enabling a faster growth of the magnetic perturbations.

Recent 1D numerical simulations of shocks including particle collisions have shown that the whistler waves excited upstream of the shock are strongly damped by Coulomb collisions (Nakanotani et al. 2022). Although this may hinder the diffusive shock acceleration process by degrading the confinement of particles at the shock front, the opposite effect of enhancing the NR mode may point towards a different conclusion. Large scale hybrid-PIC simulations with a self-consistent kinetic description of the ion populations and including particle collisions will be necessary to assess the resulting cosmic rays acceleration efficiency.

## ACKNOWLEDGEMENTS

This work was performed using HPC resources from Grand Équipement National de Calcul Intensif- [Très Grand Centre de Calcul] (Grant 2019- [DARI A0060410819]), and was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project EquipMeso (reference ANR-10-EQPX29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche.

## DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.
REFERENCES

Akimoto K., Winske D., Gary S. P., Thomsen M. F., 1993, *J. Geophys. Res.: Space Phys.*, **98**, 1419

Amato E., 2014, *Int. J. Mod. Phys. D*, **23**, 1430013

Amato E., Blasi P., 2009, *MNRAS*, **392**, 1591

Bai X.-N., Caprioli D., Sironi L., Spitkovsky A., 2015, *ApJ*, **809**, 55

Bell A. R., 2004, *MNRAS*, **353**, 550

Bell A. R., 2013, *Astropart. Phys.*, **43**, 56

Boris J. P., 1970, Acceleration Calculation from a Scalar Potential. Princeton University Plasma Physics Laboratory edn, Princeton

Chew G. F., Goldberger M. L., Low F. E., Chandrasekhar S., 1956, *Proc. R. Soc. Lon. Ser. A. Math. Phys. Sci.*, **236**, 112

Crumley P., Caprioli D., Markoff S., Spitkovsky A., 2019, *MNRAS*, **485**, 5105

Cui Y., Pühlhofer G., Santangelo A., 2016, *A&A*, **591**, A68

Ferriere K. M., 2001, *Rev. Mod. Phys.*, **73**, 1031

Forteza P., Oliver R., Ballester J. L., Khodachenko M. L., 2007, *A&A*, **461**, 731

Gary S. P., Montgomery M. D., Feldman W. C., Forslund D. W., 1976, *J. Geophys. Res.*, **81**, 1241

Gary S. P., Smith C. W., Lee M. A., Goldstein M. L., Forslund D. W., 1984, *Phys. Fluids*, **27**, 1852

Gary S. P., Wang J., Winske D., Fuselier S. A., 1997, *J. Geophys. Res.: Space Phys.*, **102**, 27159

Haggerty C. C., Caprioli D., 2019, *ApJ*, **887**, 165

Hellinger P., Trávníček P., Kasper J. C., Lazarus A. J., 2006, *Geophys. Res. Lett.*, **33**, L09101

Hunana P. et al., 2019, *J. Plasma Phys.*, **85**, 205850602

Krstic P. S., Schultz D. R., 1999, *Atom. Plasma-Material Inter. Data Fusion* (Sup. Nuclear Fusion), **8**, 0

Kulsrud R. M., 1983, Basic Plasma Physics: Selected Chapters, Handbook of Plasma Physics (Vol. 1). North-Holland Publishing Company, Amsterdam

Kulsrud R., Pearce W. P., 1969, *ApJ*, **156**, 445

Marcowith A. et al., 2016, *Rep. Prog. Phys.*, **79**, 046901

Marret A., Ciardi A., Smets R., Fuchs J., 2021, *MNRAS*, **500**, 2302

Marret A., Ciardi A., Smets R., Fuchs J., Nicolas L., 2022, *Phys. Rev. Lett.*, **128**, 115101

Matthews J. H., Bell A. R., Blundell K. M., Araudo A. T., 2017, *MNRAS*, **469**, 1849

Mignone A., Bodo G., Vaidya B., Mattia G., 2018, *ApJ*, **859**, 13

Nakanotani M., Camata R. P., Arslanbekov R. R., Zank G. P., 2022, *Phys. Rev. E*, **105**, 045209

Nanbu K., 2000, *IEEE Trans. Plasma Sci.*, **28**, 971

Ohira Y., Reville B., Kirk J. G., Takahara F., 2009, *ApJ*, **698**, 445

Onsager T. G., Winske D., Thomsen M. F., 1991, *J. Geophys. Res.: Space Phys.*, **96**, 1775

Recchia S., Galli D., Nava L., Padovani M., Gabici S., Marcowith A., Ptuskin V., Morlino G., 2022, *A&A*, **660**, A57

Reville B., Kirk J. G., Duffy P., O'Sullivan S., 2007, *A&A*, **475**, 435

Reville B., Kirk J. G., Duffy P., O'Sullivan S., 2008, *Int. J. Mod. Phys. D*, **17**, 1795

Riquelme M. A., Spitkovsky A., 2009, *ApJ*, **694**, 626

Riquelme M. A., Quataert E., Verscharen D., 2015, *ApJ*, **800**, 27

Sentman D. D., Edmiston J. P., Frank L. A., 1981, *J. Geophys. Res.: Space Phys.*, **86**, 7487

Smets R., Aunai N., 2011, Heckle. Available at: https://github.com/rochSme ts/heckle

Southwood D. J., Kivelson M. G., 1993, *J. Geophys. Res.*, **98**, 9181

Takizuka T., Abe H., 1977, *J. Comput. Phys.*, **25**, 205

Trubnikov B. A., 1965, *Rev. Plasma Phys.*, **1**, 105

van Marle A. J., Casse F., Marcowith A., 2018, *MNRAS*, **473**, 3394

Walters J., Klein K. G., Lichko E., Stevens M. L., Verscharen D., Chandran B. D. G., 2023, *ApJ*, **955**, 97

Winske D., Leroy M. M., 1984, *J. Geophys. Res.: Space Phys.*, **89**, 2673

Zirakashvili V. N., Ptuskin V. S., Völk H. J., 2008, *ApJ*, **678**, 255

Zweibel E. G., Everett J. E., 2010, *ApJ*, **709**, 1412

This paper has been typeset from a *TEX/LATEX* file prepared by the author.