{'fraction_non_alphanumeric': 0.034090497409701846, 'fraction_numerical': 0.019200811505995726, 'mean_word_length': 4.284073506891271, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The nonlinear saturation of the non-resonant kinetically driven streaming instability

A non-resonant instability for the amplification of the interstellar magnetic field in young Supernova Remnant (SNR) shocks was predicted by Bell (2004), and is thought to be relevant for the acceleration of cosmic ray (CR) particles. For this instability, the CRs streaming ahead of SNR shock fronts drive electromagnetic waves with wavelengths much shorter than the typical CR Larmor radius, by inducing a current parallel to the background magnetic field. We explore the nonlinear regime of the non-resonant mode using Particle-in-Cell (PIC) hybrid simulations, with kinetic ions and fluid electrons, and analyze the saturation mechanism for realistic CR and background plasma parameters. In the linear regime, the observed growth rates and wavelengths match the theoretical predictions; the nonlinear stage of the instability shows a strong reaction of both the background plasma and the CR particles, with the saturation level of the magnetic field varying with the CR parameters. The simulations with CR-to-background density ratios of n_CR/n_b=10^(-5) reveal the highest magnetic field saturation levels, with energy also being transferred to the background plasma and to the perpendicular velocity components of the CR particles. The results show that amplification factors>10 for the magnetic field can be achieved, and suggest that this instability is important for the generation of magnetic field turbulence, and for the acceleration of CR particles.


INTRODUCTION
Very energetic CRs (∼ 10 14 eV to 10 15 eV) are thought to be accelerated in SNR shocks. There is direct evidence that electrons are accelerated up to energies of 10 14 eV at SNR sites (Koyama et al. 1995;Allen et al. 1997;Tanimori et al. 1998;Aharonian et al. 2001;Naito et al. 1999;Aharonian 1999;Berezhko et al. 2003;Vink & Laming 2003), and the measured power law spectra of the CRs indicates Diffusive Shock Acceleration (DSA) as the most likely mechanism responsi-ble for the acceleration (Axford et al. 1977;Bell 1978;Blandford & Ostriker 1978). The acceleration of these particles up to energies of ∼ 10 15 eV through the DSA mechanism requires the existence of magnetic fields much stronger than the typical B 0 ∼ 3 µG. These strong fields have also recently been inferred from observations (Longair 1994;Berezhko et al. 2003;Vink & Laming 2003); their existence, along with the requirement of stronger fields for the DSA mechanism, suggests that a magnetic field amplification mechanism is in operation.
A possible amplification mechanism through a nonresonant instability was suggested in Bell (2004), following previous work on the resonant mode in Lucek & Bell (2000), and later extended to include multidimensional effects in Bell (2005). The non-resonant streaming instability, part of a class of streaming instabilities derived in Winske & Leroy (1984), can be described by considering a MHD model for the background plasma, and a CR induced current imposed externally; the feedback of the electromagnetic fields on the CR particles is thus neglected. Although in the linear stage of the instability λ max ≪ r LCR (λ max the fastest growing wavelength and r LCR the typical Larmor radius of the CRs), recent full PIC simulations by Niemiec et al. (2008), Riquelme & Spitkovsky (2009), Ohira et al. (2009) and Stroman et al. (2009) indicate that the feedback mechanism of the fields on the CRs is important in the nonlinear stage, and suggest that a careful study of the instability in this regime is important to determine the saturation levels of the magnetic field. Recent works by Amato & Blasi (2009) and Luo & Melrose (2009), using kinetic theory, have also shown how the saturation of this non-resonant mode depends on the details of the particle distributions.
The analysis of the behavior of the non-linear stage of the instability is very complex; full assessment of the saturation level of the magnetic field would imply a firstprinciple calculation of the shock formation, and the long term evolution of the shock precursor, including the self consistent acceleration of particles. The saturation of the instability was thus first assessed numerically in Bell (2004), where an external current was used to model the CRs, and an MHD model was used to simulate the background plasma. The saturation mechanism found was due to the tension of the magnetic field, which grows faster than the driving term; also, the size of the simulation box was seen to limit further growth of the instability (Bell 2004). More recent PIC simulation results in Niemiec et al. (2008) actually show a saturation level of δB/B ∼ 1, but consider parameters such that γ max /ω ci ≪ 1 (the ratio of the growth rate of the fastest growing mode to the CR cyclotron frequency) is not strictly maintained, which is a requirement for the development of the non-resonant parallel mode. The full PIC simulation results in Riquelme & Spitkovsky (2009) and Stroman et al. (2009), taking into account the feedback of the electromagnetic fields on the CR particles, and using m i /m e ratios up to 100 and n CR /n b down to 4 × 10 −3 , imply a saturation level of δB/B ∼ 10, which occurs when the relative drift between the CRs and the background plasma decreases. Also, in Ohira et al. (2009), similar saturation levels for the magnetic field are found.
For PIC simulations, therefore, the instability saturates when CRs loose some of their bulk momentum to the background plasma, which starts to drift more rapidly in the direction of CRs until the velocities of the two species converge and the driver for the instability is eliminated. As shown in kinetic modeling by Luo & Melrose (2009), the reduction of the CR streaming motion is due to CR resonant diffusion in nonresonant magnetic turbulence, which has been recently observed in Stroman et al. (2009). The nonlinear evolution of the turbulence observed in fully kinetic simulations is also in qualitative agreement with the predictions of quasi-linear calculations for non-resonant modes derived for non-relativistic beams in Winske & Leroy (1984). However, the accurate modeling of the late time evolution of the system, beyond the saturation, is limited in the kinetic simulations, when the assumption of plasma homogeneity is no longer valid. As argued in (Ohira et al. 2009), when the background plasma veloc-ity approaches the CR bulk velocity, the density in a real scenario also changes, modifying the underlying conditions for the simulation model. This reinforces the necessity for a thorough analysis of the behavior of this instability in its non-linear stage, covering a wide range of dynamical scales, and exploring the saturation limits in different regimes.
Here, we present multi-dimensional simulation results of the non-resonant streaming instability, using the kinetic ion fluid electron hybrid model implemented in dHybrid (see Gargaté et al. (2007) for numerical implementation details). An important advantage of the hybrid simulations is to enable the study of the instability on the ion time scale, neglecting the high-frequency modes associated with the electrons; realistic density ratios down to n CR /n b = 10 −5 can then be used, along with realistic ion-to-electron mass ratios, and simulations can be run to a point well beyond the linear stage, into the saturated state. The variation of the density ratio has an impact on the behavior of the nonresonant mode, with larger values implying the generation of different modes (Bell 2004;Niemiec et al. 2008;Riquelme & Spitkovsky 2009). Though the saturation mechanism is independent of the initial linear instability (Stroman et al. 2009), these modes lead to smaller saturation amplitudes and thus may not be directly relevant for the typical SNR shock scenarios. Also, in the hybrid simulation results shown, both the background and the CR ions are modeled kinetically, which enables the study of the feedback mechanism on the CR population. Results can then be easily compared with both the external current-driven MHD simulations, where feedback on the CR population is not modeled, and the fully kinetic simulations, which use small ion-to-electron mass ratios, and larger n CR /n b ratios. This paper is organized as follows. In section 2, the parameters used in the simulations are discussed, and the results concerning the evolution of the linear and nonlinear stages of the instability are presented. The saturation mechanism is discussed in section 3, and compared with the most recent results in the literature. Finally, in the last section, we present the conclusions and outline future research directions.

EVOLUTION OF THE NON-RESONANT INSTABILITY
The original theoretical model of the non-resonant streaming instability, as detailed in Bell (2004), considers that a magnetized background plasma with density n b is stationary in the SNR shock upstream reference frame. The CR particles stream along the background magnetic field lines B with a drift velocity v sh , similar to the shock velocity, generating a current that drives the growth of the perpendicular components of the magnetic field B ⊥ . For the canonical parameters in the literature (Bell 2004), of B = 3 µG, n b = 1 cm −3 , v sh = 10000 km/s, and r LCR ∼ 1.1 × 10 13 m, the Alfvén velocity is v A = 6.6 km/s, and the Alfvén Mach number is M A = v sh /v A = 1515 ≫ 1. Under these conditions, the reduced wave equation ω 2 − v 2 A k 2 ± B j k/ (n b m) = 0, where m is the mass of the background ions, j is the zeroth order current imposed by the CRs, k is the wavenumber and ω is the frequency of the mode, yields a purely growing mode with a , and a maximum growing wavenumber k max = 1/2 B j /(n b m v 2 A ) with a growthrate of γ max = k max v A . For the nonresonant mode to be predominant, then 1 < kr LCR < ζv 2 sh /v 2 A , with ζ = B j r LCR /(n b m v 2 sh ); for the canonical parameters it follows that k max r LCR = 3691, and ζv 2 sh /v 2 A = 7381. Moreover, γ max /ω ci ≪ 1 is maintained, which is an essential requirement for the development of the non-resonant parallel mode.
The results are presented in normalized units, with the magnetic field normalized to the background field B 0 = B , the density normalized to the background plasma density n 0 = n b , the velocities normalized to the Alfvén velocity v A = B 0 /(n 0 m p µ 0 ) 1/2 , the spatial dimensions normalized to the ion inertial length c/ω pi , and time normalized to the inverse maximum growth rate γ −1 max . The normalized background magnetic field in the simulations is B = B x = 1 B 0 , and the normalized background plasma density is n b = 1 n 0 , corresponding to SNR scenarios where B /n b = 0.4 µG cm 3 .
The main parameters for the simulations are described in Table 1. The CRs are modeled through a constant external current in run A, and with kinetic particles in all other runs (the kinetic-CR runs). The value of the initial CR current is always j = n CR e v sh = π/20 [n 0 e v A ] (where [n 0 e v A ] are the normalized current density units), which yields a maximum growing wavelength λ max = 80 c/ω pi (k max = π/40 ω pi /c). For the kinetic-CR runs, the total velocity in the CR population is initially set to v sh + v iso , with v iso the isotropic velocity component, so that the whole CR population drifts with the shock velocity, as observed in the upstream reference frame. For runs B 1 through B 3 the CR-to-backghround density ratio was varied, along with the drift velocity v sh , and for runs C 1 through C 3 , the v iso was varied, enabling the analysis of the dependence of the instability on the scale separation between the typical Larmor radius of the CRs, and λ max . Finally, in run D, the parameters from run B 1 were used in a 3D setup. For all the runs the time step was chosen to be ∆t = 0.001 ω −1 ci , as to account for the cavitation in the background density plasma, since low density zones can render the hybrid (and MHD) algorithms unstable (Lipatov 2002;Gargaté et al. 2007). Full convergence tests have been performed to assure that the saturation level of the magnetic field was not affected by the choice of the time step.
The simulation box size is L x = L y = 1280 c/ω pi for the 2D runs, and L z = 640 c/ω pi for the third dimension in the 3D run presented. The L x and L y dimensions are then ∼ 16 × λ max , corresponding to L x /L diff ∼ 1, where L diff ∼ D/v sh is the CR diffusion length scale in the Bohm diffusion limit for the typical CR velocities used in the simulations. This choice of parameters then allows for the CR Larmor radius to be fully resolved, and is such that the non-resonant mode wavelength can grow to be of the order of r LCR in the non-linear regime, enabling the analysis of the back-reaction of the generated electromagnetic fields on the CR particles.
The diverse parameters chosen for the four sets of runs (A through D) does not directly reflect typical physical conditions found in young SNR shock precursors for all the runs; instead a wide dynamical range of parameters was analyzed. The results of run A are most directly comparable to the simulations by Bell (2004), in which the CR current was imposed externally, thus excluding any feedback of magnetic turbulence on theCRs. For the B set of runs, the driving current is maintained constant while the driving energy of the CR beam is increased from run B 1 to run B 3 by increasing v sh and decreasing n CR ; as such v iso /v sh is smaller than the typical values found in young SNR shock precursors, especially for runs B 2 and B 3 . The C set of runs explores the behavior of the instability when progressing from the non-resonant to the resonant limit of the CR streaming instability by lowering the v iso velocity and thus the r LCR /λ max ratio.
Comparison of results between 3D and 2D runs showed no significant difference in the development of the instability. This motivated the analysis of the system using a 2D simulation setup. For run D, the measured wavelength of the instability was λ = 80 c/ω pi , and the growth rate γ ∼ 0.889 of the non-resonant mode is the same as for all the runs in the B set. Figure 1 shows the evolution of the B z magnetic field component and the density of the background plasma for run B 2 ; Fig. 1 a) and c) correspond to the linear stage and show the well defined wavelength λ = 80 c/ω pi in the magnetic field, and a low amplitude modulation of the background density. In the nonlinear stage, cavities are formed in the background plasma, due to turbulent plasma motions induced by the electromagnetic waves being generated. The well defined wave structure of the linear stage disappears in Fig. 1 b), with magnetic field compression zones being formed, which are correlated with the density enhance- ments at the same positions, Fig. 1 d), confirming the results of (Riquelme & Spitkovsky 2009) for lower mass ratios. Also, very similar results are observed for runs B 1 and B 3 .
In Fig. 2, the main characteristics of the instability that agree with expectations of a quasi-linear theory are observable. At t ∼ 9 γ −1 the instability becomes nonlinear, and saturation of the linear mode is observed a short time after, around t = 10 γ −1 (Fig. 2 b) and d).
In the linear stage, the growth rate from Fig. 2 d) yields γ = .889 which deviates from the theoretical value by 11%. It is also clear from Fig. 2 a) that the peak power is for k ∼ 0.078 ω pi /c up to t = 8 γ −1 (dotted line), and that at t = 16 γ −1 the peak shifts to k ∼ 0.044 ω pi /c (dashed line).
In the nonlinear stage, Fig. 2 b) and c), the increase in power over large k is consistent with the magnetic field structures in Fig. 1 b), and indicates an increase in the turbulence level in B z . Figure 2 c) also shows that the waves move in the x direction at a velocity v ∼ 1.8 v A in the nonlinear stage.

SATURATION MECHANISM
The maximum energy level in the perpendicular magnetic field components, at the end of the linear stage, is similar for the runs B 1 , B 2 (see Fig. 3), and B 3 . In the linear stage of the instability, the behavior of run A is identical to the kinetically-driven runs, B 1 and B 2 . In the nonlinear stage, however, the magnetic field energy peaks at t = 10 γ −1 for runs B 1 and B 2 , while for run A the magnetic field energy still grows for t > 10 γ −1 . The growth beyond t = 10 γ −1 for run A is due to the continuous injection of energy into the system; a similar growth is observed for run B 3 , at a slower rate, which is due to the excess of free energy in the CR population for this run. Finally, run C 1 (Fig. 3 b) represents an hybrid scenario, between the resonant and the nonresonant modes of the instability (r LRC /λ max ∼ 1.63), and thus the magnetic field growth rate is smaller in the linear stage. . Frame c) shows the current density components for run B 1 , and frame d) shows the current density components for run B 2 ; the vertical arrows, from left to right, indicate points where the background plasma current, Jx is: at a minimum, zero, at a maximum, and zero, respectively. Lines are labeled in the plot for ease of interpretation. The energy scale is normalized to the total energy in run B 3 .
The CR population velocity distribution is nearly isotropic at the end of all the simulations (as seen in the upstream / simulation frame). In the linear stage, the growth rates are identical between runs in the nonresonant regime, and the background plasma also gains energy at the same rate in the perpendicular direction (Fig. 3 c) and d). Unlike what is predicted from theory, v x of the background plasma increases, and the parallel energy becomes comparable to the perpendicular energy. With this energy increase, the plasma is no longer magnetized, and thus the nature of the non-resonant mode changes.
In order to understand the dependence of the nonresonant mode on the scale-separation between λ max and r LCR , the isotropic velocity was varied in runs C 1 through C 3 . The typical CR Larmor radius varies as r LCR = {130.6, 261.3, 408.2} c/ω pi (and 653.2 c/ω pi for run B 1 ). Figure 3 and Fig. 4 a) and b) show that for r LCR /λ max ratios greater than ∼ 3.2 (run C 2 ), the behavior is consistent with the non-resonant mode, and only run C 1 (r LCR /λ max = 1.63) results in an hybrid scenario. Analysis of Fig. 4 c) and d) also shows that for runs B 1 and B 2 , at t ∼ 9 γ −1 , the J x current density component for the CRs is not close to zero. This means that the saturation mechanism does not depend directly on the amount of free energy in the CR population. Instead, the saturation of the instability is due to the energy gain of the background plasma, which results in an average bulk flow in the CR propagation direction, and in the de-magnetization of the plasma.
The relation between J x for the background plasma, and both the onset and the peak in the magnetic field energy level of the instability is highlighted in Fig. 4 c) and d). The peak of J x , for the background plasma, marks the end of the linear stage of the instability. Likewise, when J x is at a minimum (t ∼ 3 γ −1 ), the instability is entering the initial linear stage. Finally, for run B 1 , Fig. 4 a) and c), the energy in the CRs at the saturation point is transferred to the background plasma and to the magnetic field; for run B 2 , Fig. 4 b) and d), some of the energy is converted into perpendicular energy in the CR population, as J ⊥ increases from t ∼ 10.8 γ −1 onwards.

DISCUSSION AND CONCLUSIONS
The identification of the nonlinear saturation mechanism for the non-resonant streaming instability is important for the determination of the maximum magnetic field amplification in SNR shock scenarios. Amplification factors of B ⊥ /B ∼ 10, with local temporary peaks of B ⊥ /B ∼ 25 are observable in the current hybrid simulations; similar results were also recently shown in the kinetic simulations of Riquelme & Spitkovsky (2009) Ohira et al. (2009) and Stroman et al. (2009, with reduced temporal and spatial scales, mass ratios and density ratios. The hybrid simulation results presented expand the current knowledge of the non-resonant streaming instability by extending the dynamical range under study, and by performing a detailed study of the instability under different driving regimes. Leveraging on the hybrid simulations presented here, it is possible to scan the parameter space and analyze the nonlinear stage of the instability in detail, using realistic n CR /n b and m/m e ratios. Close examination of Fig. 3 and Fig. 4 shows a number of important points, relating to the saturation mechanism. As the n CR /n b ratio is decreased towards 10 −5 , the energy in the CRs increases, as in our setup the current is initialized at the same level in all simulations; the increase in the free energy of the CRs, associated with the higher drift velocity v sh , does not change the peak energy level of the magnetic field. The excess free energy is instead transferred to the background plasma, and also to the perpendicular velocity components of the CR population (Fig. 3 b), run C 1 ).
The reaction of the background plasma is critical in the nonlinear stage of the instability. The average local Alfvén velocity (v A ∝ B / √ n calculated in each cell and averaged over the entire simulation box) is approximately constant over the linear development of the instability. In the nonlinear stage, v A increases with the local magnetic field B , and with the background density varia-tions (which are in phase), and thus v A ∝ δB/B . The background plasma accelerates in the CR propagation direction, and the current density peaks, marking the end of the linear stage (Fig. 4). At this point in time, the CRs v x velocity component is still v x ∼ 1500 v A , well above the Alfvén velocity and above the background plasma drift velocity, which is sub-Alfvénic. This shows that the saturation is not dependent on the amount of free energy in the CRs, but instead results from the reduction of the scale separation of the background plasma, and the CRs.
The CR's J ⊥ increases for run B 2 , Fig. 4 d), as it is substantially lower than J x initially, and the distribution becomes isotropic after the saturation, at t ∼ 15 γ −1 . For run B 1 in Fig. 4 c), the distribution is also isotropic at the end, but the CRs J ⊥ velocity component is not affected. Increasing the free energy in the CRs does not affect the saturation level of the magnetic field significantly, as long as the driving current is maintained.
The hybrid simulation results presented thus show that, for a broad range of parameters, the nonlinear growth of the streaming instability is independent on the energy of the driving CR population, and depends only on the current carried by the CRs. The amplification factor of ∼ 10 for the perpendicular magnetic field above the seed B feed indicates that the mechanism is relevant for the magnetic field amplification by CR particles in SNR shocks. Beyond the saturation point, for t > 16 γ −1 , the CR particle distribution becomes isotropic, with the instability saturating after 10 e-foldings. The de-magnetization of the background plasma, as it gains energy, hinders further growth of the non-resonant mode; when both the CRs and the background plasma are unmagnetized, the instability should behave more like the Weibel instability (Weibel 1959).
Our results thus indicate that the instability begins to saturate when the current carried by the background ions is similar to the current carried by the CR population (see Fig. 4). This results in a final average velocity for the background plasma of v b ∼ n CR /n b v sh so that V b ∼ 10 −5 v sh for realistic parameters. At this point, the instability is in the non-linear stage, and the magnetic field growth rate is much lower (t > 9 γ −1 in Fig. 3 c) and Fig 4 a) and c) for run B 1 ). In fact, even in the linear stage of the instability, energy is being transferred to the perpendicular components of the background plasma at a greater rate than into the magnetic field perpendicular components (see Fig. 3 c), and thus at some point the background plasma de-magnetizes and the instability saturates.
Previous simulation PIC results (Riquelme & Spitkovsky 2009;Ohira et al. 2009;Stroman et al. 2009) show similar saturation levels for the magnetic field δB/B ∼ 10: the instability saturates when the background plasma is moving with a bulk velocity comparable to the CR population drift speed, which is equivalent to saying that the currents carried by the background plasma and the CRs are similar (since n CR ∼ n b , and thus J CR = n CR e v sh ∼ J b = n b e v b ). It has been claimed (Ohira et al. 2009) that this saturation limit might be numerical rather than physical because in a real shock precursor a significant variation in the background plasma velocity should be concurrent with a variation in the background plasma density.However, here we show for m i ≫ m e and n CR ≪ n b that the instability saturates due to the de-magnetization of the background plasma; at that point the currents carried by the background ions and the CRs are comparable, although the velocities differ significantly (since in our hybrid simulations n CR ≪ n b implies that v b ∼ v sh n CR /n b from the equality of the currents at the saturation). Our simulations thus reinforce the result δB/B ∼ 10 at saturation, and indicate that the instability should be important in a young SNR shock precursor. The simulations present in the literature make assumptions about the physical scenarios, and therefore significant care should be taken when interpreting and comparing the results. For instance, the back-reaction of the fields in the CRs is not accounted for when using external currents as drivers, in full PIC simulations mass ratios and density ratios far from the physical conditions have been explored, and the dynamics of electrons is not accounted for in hybrid simulations. Thus, our results further indicate that first-principle simulations of SNR shocks and particle acceleration should be attempted in a future work, since the different simulation methods can capture the instability, but differences might still be found when the simulation setup and the models are refined in order to more closely match the experimental conditions.
Further analysis for SNR shocks will be possible by considering the saturation mechanism when unperturbed CRs are continuously injected into the simulation box. This will allow for an improved comparison of the SNR shock scenario, where CR particles are also thought to be continuously injected from the shock into the upstream medium. This will be explored in a future publication, leveraging on the unique characteristics of the hybrid model. Finally, with hybrid simulations, it will also be possible to determine the dependence of usual characteristics of the instability with the distance to the shock, and to study in detail the energy profile of the accelerated CR particles. The simulations presented in this paper were produced using the IST Cluster (IST/Portugal). The work of JN is supported by MNiSW research project N203 393034, and the Foundation for Polish Science through the HOMING program, which is supported by a grant from Iceland, Liechtenstein, and Norway through the EEA Financial Mechanism. The authors also want to thank KITP (UCSB) where part of this work was done, partially supported by NSF under Grant No. PHY05-51164, and A. R. Bell for the useful discussions on the instability mechanism.