# How large can the electron to proton mass ratio be in Particle-In-Cell simulations of unstable systems?

###### Abstract

Particle-in-cell (PIC) simulations are widely used as a tool to investigate instabilities that develop between a collisionless plasma and beams of charged particles. However, even on contemporary supercomputers, it is not always possible to resolve the ion dynamics in more than one spatial dimension with such simulations. The ion mass is thus reduced below 1836 electron masses, which can affect the plasma dynamics during the initial exponential growth phase of the instability and during the subsequent nonlinear saturation. The goal of this article is to assess how far the electron to ion mass ratio can be increased, without changing qualitatively the physics. It is first demonstrated that there can be no exact similarity law, which balances a change of the mass ratio with that of another plasma parameter, leaving the physics unchanged. Restricting then the analysis to the linear phase, a criterion allowing to define a maximum ratio is explicated in terms of the hierarchy of the linear unstable modes. The criterion is applied to the case of a relativistic electron beam crossing an unmagnetized electron-ion plasma.

## I Introduction

The dynamics of collision-less plasma far from its equilibrium is frequently examined with particle-in-cell (PIC) simulations. The unique capability of PIC codes to model such systems from first principles on macroscopic scales implies that they can bridge the gap between theory and experiment. For example, just a few years ago it was still unclear if relativistic shocks exist. It was not known whether the motion energy could be dissipated rapidly enough to sustain the shock discontinuity Brainerd2000 (); Waxman2006 (). Such shocks have not yet been observed directly, because they do not exist in solar system plasma. Recent PIC simulations could shed light on how they develop in response to the filamentation (Weibel) instability SilvaApJ (); Spitkovsky2008 (); Chang2008 (). The particle acceleration and the generation of electromagnetic radiation within the context of active galactic nuclei GALLANT1992 (); FrederiksenApJ2004 (); DieckmannMNRAS2006 (); Spitkovsky2008a (), supernova remnants HOSHINO1992 (); DieckmannApJ2009 () or gamma ray bursts Jaroschek2004 (); Jaroschek2005 (); MedvedevApJ2005 (); DieckmannMNRAS2006a () have also been investigated. PIC simulations are now instrumental in investigating the plasma thermalization within solar flares KarlickyApJ (); Karlicky2009 () and the dynamics of magnetic reconnection CheYoon (); Lapenta (). On a completely different length and density scale, the fast ignition scenario for inertial confinement fusion Tabak () has prompted within the last decades many numerical works focusing on the propagation of charged particle beams in a collisionless plasma Pukhov96 (); Ren2004 ().

Current simulations employ billions of computational particles, placing physically realistic PIC simulations within our reach MedvedevSpit2009 (). The inclusion of ions in PIC simulations nevertheless remains a formidable challenge. As long as the system under scrutiny involves only electrons and positrons with the mass , the time scale that must be resolved is typically the inverse electronic plasma frequency . Running the simulation for hundreds or thousands captures the evolution of the system way beyond its linear phase. Mobile protons or ions in the simulation result in an additional and much longer timescale. A PIC simulation must then resolve many inverse proton plasma frequencies and cover a time interval that is times longer, if the plasma is unmagnetized and if protons are the only ion species. A further penalty is introduced by the larger spatial scales of the ion structures. The size of the ion filaments is, for example, comparable to the ion skin depth , while that of the electron filaments is . In principle, the simulation box size that is necessary to model electron-proton plasmas increases compared to that required by leptonic plasmas by a factor , where is the number of resolved spatial dimensions.

For this reason and until now, PIC simulations that use the correct electron-to-proton mass ratio are restricted to 1D systems DieckmannMNRAS2006 () and to 2D simulations that resolve a limited spatio-temporal domain HondaPRL (), while 2D PIC simulations that cover a large domain with regard to the ion scales or even 3D PIC simulations normally resort to reduced ion masses between 10-100 electron masses Spitkovsky2008 (); Niemiec2008 (); Martins2009 (). Ion masses of up to 1000 electron masses have been used in a 2.5D simulation Spitkovsky2008 () thanks to a low number of particles per cell, which is beneficial for the scalability of a domain-decomposed PIC code. Multidimensional plasma simulations that employ the correct mass ratio and capture the largest ion scales are possible, if the electron dynamics does not have to be resolved accurately. Implicit PIC schemes can dissipate away the energy contained in the smallest scales in a form of Landau damping Lapenta2006 (); Noguchi2007 (). The cell size can then be increased beyond the plasma Debye length, without restricting the physical accuracy of the large-scale dynamics. However, if both the electron and the ion dynamics must be resolved simultaneously, the implicit PIC codes are equally costly as the explicit ones. The speeding up of the plasma evolution through a reduced ion mass will remain a necessity in particular for 3D simulations. It is thus important to assess how the plasma dynamics changes with this parameter.

Various studies exist that demonstrate the importance of the electron-to-ion mass ratio for the plasma dynamics in several types of plasma processes. Parametric studies of plasma shocks have addressed this issue both in the non-relativistic Scholer2003 (); Scholer2004 (); Shimada2005 () and relativistic Spitkovsky2008 () regimes. Simulation studies of the interplay between electron phase space holes with the ions and its dependence on the mass ratio can be found in califanoIOns2007 (); EliassonPRL2004 (). The impact of the mass ratio on the reconnection of magnetic field lines and the associated particle acceleration has been investigated in Ref. Ricci2004 (); GuoandLu2007 (). However, these studies related to the effects of a reduced ion mass focus primarily on the nonlinear evolution of the simulation.

This article is a first systematic study of the consequences of a reduced ion mass within a theoretical framework. The impact of the mass ratio on the nonlinear coupling of the plasma dynamics across the different scales is not considered here; it is too complex and multifaceted. An example would be the enlargement of the foreshock of a perpendicular shock with an increasing ion mass, which influences the resulting instabilities and the thermalization of the shock-reflected ion beam Scholer2003 (); Scholer2004 (); Shimada2005 (); Lee2004 (). We study here the spectrum of linearly unstable waves, which should depend on the mass ratio between the ions and the electrons. Ions only one time “heavier” than electrons are obviously too light as they behave like positrons, not like protons. Is it therefore possible to draw a line, from which ions will start being “too light” to represent protons?

Even before answering this question, one could ask whether the PIC simulation plasmas could be governed by some similarity laws involving the mass ratio. Similarity theory has been applied successfully in hydrodynamics. It allows us to predict certain properties of an object from experiments performed with its miniaturised model. Well-known cases of such experiments involve pumps, turbines or aircrafts BrennenPumps (); Birkhoff (). Similarity laws have also been derived for magnetic confinement fusion (see petty2008 (), and references therein) or relativistic laser-plasma interactions and laboratory astrophysics Drury2000 (); Ryutov2001 (); gordienko2005 (). Similarity laws would allow us to compensate a reduced mass ratio with some other parameter, by which the computational efficiency can be altered.

We show in section 2 that it is not possible to derive an universal description of the growth rate spectrum, which is not explicitly dependent on the mass ratio. Section 3 will therefore aim at providing a restricted solution to the problem. PIC simulations typically probe the long term nonlinear evolution of an unstable beam-plasma system. The unstable spectrum usually contains more that one unstable mode, and these modes grow at different exponential rates. A hierarchy of unstable modes can be established in terms of their growth rate, and a criterion can be imposed on the mass ratio by demanding that this hierarchy be preserved. The consequences of this condition are then explicitly calculated for the case of a cold relativistic electron beam passing through a un-magnetized and cold plasma califano3 (). Section 4 is the discussion, which brings forward a possible explanation why the shock formation in Ref. Spitkovsky2008 () does apparently not depend on the mass ratio.

## Ii Absence of a similarity law

Consider a problem that involves independent dimensional variables , and fundamental dimensions such as meter, second, etc. The so-called Buckingham’s method Buckingham () to reduce the number of variables and to derive similarity laws is the following.

Identify the pairs of ’s that share the same physical unit. If this is the case for the variables and , then replace by in the list of variables. After iterating this process for any such pairs, we are left with the modified set of variables where . Buckingham’s “ theorem” then states that any unknown function of the form

(1) |

can indeed be cast under the form,

(2) |

where the variables are dimensionless products of the initial ’s, and the number of fundamental dimensions among .

Consider as a simple illustration a swinging pendulum with the mass (kg) and the length (m), which oscillates with the constant period (s) in a gravitational field (m/s). The first step of Buckingham’s method, namely the pairing of variables that share the same dimension, can be skipped here since all 4 variables have different dimensions (kg, m, s, m/s). We thus have here (4 variables) and only (kg, m, s) as does not add any extra fundamental dimension to the problem.

Buckingham’s theorem states in this case that any function can be expressed as , so that the problem is eventually a function of one single dimensionless parameter. A subsequent dimensional analysis shows that the universal parameter must be a power of . Clearly, the mass cannot participate in the dimensionless parameter, because no other variables could cancel its physical unit. The period is thus independent of the mass and only a function of the ratio . Note that the theorem does not distinguish between “input” variables (what is known, e.g. ) and “output” variables (what is looked for, e.g. ). Each quantity is treated in the same way and all contribute to and to .

Turning now to the present problem, we see from the first step of Buckingham’s method that a similarity law without an explicit dependence on the mass ratio cannot exist. Whatever the list of variables describing the problem may be, the mass parameters will be a part of it. The first step of the process will just replace by , and Buckingham’s theorem reduces the number of variables left once all the dimensionless trivial ratios have been formed. At any rate, Buckingham’s theorem states that the mass ratio remains as an explicit parameter in the final reduced set of parameters.

There is therefore no hope of unraveling similarity laws connecting two systems (A) and (B) with different mass ratios. A change of the mass ratio cannot be “compensated” by a shift of the other variables. Buckingham’s analysis of the problem proves an intuitively simple reasoning: electrons and ions define different time scales in terms of their respective mass. The time evolution of the system can be normalized to any one of them, but it cannot fit both at the same time. Note that although the rest of the article focuses on the linear regime of an electron beam plasma system, the present conclusion is very general, and valid for the overall evolution of any kind of system comprised of two species.

Let us initially assume that the ions are immobile, yielding an electron-to-ion mass ratio of zero. Electrons are therefore the only population bringing a mass into the parameter list. The electron mass will therefore appear among the in Eq. (1). Buckingham’s theorem here states that these dimensional variables can be replaced by dimensionless variables . Because no mass ratio can appear among the ratios which are arguments of the function in Eq. (2), the underlying equations can not rely explicitly on the electron mass. Equation (2) shows that must have been “absorbed” by one of the dimensionless ’s. This is precisely what is observed when dealing with such questions: the time parameter is frequently normalized to the electronic plasma frequency which includes the electron mass.

## Iii Preserving the growth rate hierarchy of the unstable modes

In ultra-relativistic laser-plasma interaction, similarity theory states that laser-plasma interactions with different and are similar, as soon as the similarity parameter is the same gordienko2005 () ( is the laser amplitude, is the plasma electron density, and is the critical density for a laser with frequency .

The previous Buckingham analysis has demonstrated that there cannot be any such similarity parameter in the problem we consider here. As soon as the ions are allowed to move, the electron to ion mass ratio must appears explicitly in any list of dimensionless parameters describing the system. Two systems differing only by their mass ratio will not evolve similarly.

A deviation of the simulation dynamics from the true plasma dynamics is acceptable, as long the modifications are only quantitative. The simulation can in this case still provide important qualitative insight into the plasma evolution, which can not be obtained by any other means. However, somewhere in between the mass ratio of 1/1836 and the (positron) mass ratio of 1, a line must be crossed when even this is not the case any more.

For an unstable system with ratio between the ion and the electron mass, the unstable spectrum is comprised of all the modes with wavenumbers with an amplitude that grows at the exponential (positive) growth rate . Among these modes, the most unstable mode defined by,

(3) |

plays a peculiar role because it is the one which growth determines the outcome of the linear phase. The evolution of , as a function of , may be continuous, or not. For clarity, let us consider the example studied in Sec. III.1, of a one dimensional beam-plasma system. The dispersion equation in this case gives two kinds of unstable modes: the two-stream modes and the Buneman modes. Let us assume that we have plasma parameters that are such that the dominant mode is a two-stream mode. It is possible to find a range of mass ratios , for which the two-stream mode always grows fastest. In this case, evolves continuously with . A variation of the mass ratio may trigger in another case a transition from the two-stream regime to a Buneman regime where is the wavenumber of a Buneman mode. Here, the evolution of will be discontinuous and we will talk about an altered mode hierarchy.

Changing the mass ratio in such a way that the mode hierarchy is altered will thus result in a different plasma evolution during the linear growth phase of the instabilities. In our 1D example, the typical size of the patterns generated in the early evolution will change abruptly by a factor , where are the beam and plasma electronic densities respectively. For the 2D system considered in Sec. III.2, a switch from a two-stream to a filamentation regime results in the generation of magnetized filaments instead of electrostatic stripes.

Although we focus in what follows on a specific setup, most kinds of beam-plasma systems encountered in the literature also exhibit more than one type of unstable modes califano3 (); Bell2004 (); bretApJ2009 ().

The criterion we propose is that the modified mass ratio must not alter the mode hierarchy. Note that this is a necessary but not a sufficient condition. If the mode hierarchy is altered, then the evolution of the system should be affected as well. However, even a similar linear phase could result in a different non-linear long-term evolution prompted by a different mass ratio.

Even if thermal effects are neglected, the analysis of the full spectrum of unstable waves is involved for energetic astrophysical plasmas bretApJ2009 (). The identification of the fastest-growing wave mode requires the evaluation of the full three-dimensional spectrum of wave vectors califano1 (); califano2 (); califano3 (); BretPRL2005 (). In what follows, the proposed criterion for the mass ratio is applied to the simple and generic, yet important, system formed by a relativistic electron beam that passes through a plasma with an electronic return current. The return current initially cancels the beam current and the ion charge density cancels the total electronic one. Since some PIC simulations are still performed in a 1D geometry, we start by analyzing the 1D case before we turn to the more realistic 2D and 3D ones.

### iii.1 Relativistic electron beam - 1D simulation

We consider a relativistic electron beam with the density , the velocity and the Lorentz factor , which passes through a plasma with the ion density and the electron density with . The drift velocity of the electrons of the background plasma fulfills . For the stability analysis, we consider the response of the system to harmonic perturbations . We assume that the particles move along the -axis and we consider only wavevectors with . All plasma species are cold. The dispersion equation is readily expressed as the sum of the contributions by the beam, by the return current and by the ions with mass .

(4) |

where is the Lorentz factor of the return current. We introduce with the dimensionless variables

(5) |

The dispersion equation expressed in these variables is

(6) |

As long as the return current remains non-relativistic with . For the strictly symmetric case with , the return current becomes relativistic with . The dispersion equation (6) defines two kinds of unstable modes BretPoPIons (). The two-stream instability is driven by the two electron beams. In the limit this instability has its maximum growth rate at the wavenumber

(7) |

The unstable Buneman modes arise from the interaction of the electronic return current with the ions. These additional modes grow for at the wavenumber

(8) |

Figure 1 displays the growth rate curves obtained from Eq. (6). Both wave branches share the same -interval if . Equations (7,8) show how the mode hierarchy relies explicitly on the mass ratio. Only the growth rate of the Buneman modes scales like . Changing the mass ratio can thus change the mode hierarchy.

Figure 2 depicts the range of parameters where two-stream and Buneman modes govern the spectrum for various mass ratios . The separatrix between the domain is plotted for , 1/100, 1/400 and 1/1836. The Buneman modes grow faster below the curve, while the two-stream modes are dominant above. The separatrix between the two domains is given by

(9) |

Let us assume that we run a 1D PIC simulation from the parameters pictured by the circle labeled “PIC1”. For , the corresponding system lies in the two-stream region. As we increase , the growth rate of the Buneman instability increases relative to that of the two-stream instability. For sufficiently light ions, the Buneman instability can even outgrow the two-stream instability. Given some simulation parameters , the largest mass ratio that leaves the mode hierarchy unchanged is readily calculated from Eq. (9) if . One can see that a mass ratio as high as 1/30 is allowed only for weakly relativistic systems. If a PIC simulation uses the parameter values denoted by “PIC2”, the present criterion does not restrict the mass ratio. Any value larger than would be in favor of the Buneman modes, which are already governing the system.

### iii.2 Relativistic electron beam - 2D and 3D simulations

The previous reasoning is now expanded to a 2D geometry. It is equivalent to a full 3D geometry with regard to a linearized theory, as long as the system is cold and does not have two distinct symmetry axes. An example is a beam velocity vector that is not aligned with the magnetic field direction. Here, the forms the sole symmetry axis and it defines one direction. A second dimension takes into account the unstable modes with wave vectors that are not parallel to . These modes compete with the two-stream and Buneman modes. One finds the Weibel (or filamentation) modes for , which could play a major role in the magnetic field generation that is necessary to explain Gamma Ray Bursts Medvedev1999 (); Brainerd2000 (). For obliquely oriented wave vectors, the so-called “oblique modes” are likely to govern parts of the relativistic regime BretPRL2008 ().

The dispersion equation is more involved in 2D than in 1D, because unstable modes are generally not longitudinal (i.e. electrostatic with ). While oblique unstable modes have been known to exist for some decades now Watson (); Bludman (); fainberg (), the first exact cold fluid analysis of the full unstable spectrum was only recently performed by Califano et. al. califano3 (). The dielectric tensor is computed exactly from the Maxwell’s equations, the continuity equation and the Euler equation for the three species involved. We choose and a wave vector in the () plane. The normalized wave vector from Eq. (5) is now extended to two dimensions,

(10) |

The dielectric tensor has been computed symbolically using a Mathematica Notebook designed for this purpose BretCPC (). The dispersion equation reads now

(11) |

where the tensor is specified in the Appendix. The dispersion equation is an 8th degree polynomial. Its numerical solution is straightforward. Figure 3 shows a plot of the growth rate in terms of for the same parameters as Fig. 1. The Weibel or filamentation modes are characterized by , the two-stream and Buneman modes by , while the oblique modes constitute the remaining spectrum. The ridge in the growth rate map at low stems from the interaction of the two electron beams, while the interaction of the ions with the bulk electrons is responsible for that at larger ’s.

The growth rates of the Weibel modes and of the oblique modes can be estimated for immobile ions with and for low as fainberg (),

(12) |

These expressions must be corrected in a nontrivial way in the ultra-relativistic limit as approaches unity. Prior to the formulation of our criterion for the mass ratio , we elucidate the hierarchy map and how it evolves with . Figure 4 pictures the separatrices of the domains in the parameter space for and for . Weibel modes tend to govern the high beam density regime, and slightly expand the mildly-relativistic (around ) part of their domain when grows. Buneman modes modes govern the lower-right corner of the graph, and being scaled like , increase their domain as well with . As a result, the domains governed by the oblique modes shrinks with a growing . The two-stream modes actually never govern the system, because they are outgrown by the oblique modes as soon as .

The parameter space diagram reveals a triple point, at which the separatrices merge. For , its coordinates are . Figure 5 shows how the triple point location evolves towards the ultra-relativistic regime as the ion mass is increased to that of a proton. This triple point is not likely to be important in astrophysical flows, since even the Lorentz factors of GRB jets do not reach such high values. However, it can become an issue in PIC simulations, where mass ratios of 30 and Lorentz factors of a few tens are not uncommon.

We now compute the largest that leaves unchanged the mode hierarchy for a given parameter set , and display the result on Fig. 6 (in fact, the smallest inverse mass ratio is plotted, for better clarity). If, for , a system lies in the Buneman region, then increasing will not change the dominant mode. The same holds for the systems pertaining initially to the Weibel zone. But the system represented by “PIC1” on Fig. 4 remains governed by the oblique modes only up to a certain value of the mass ratio, beyond which it goes over into the Buneman domain. The same can be said for “PIC2”: initially lying in the oblique domain, it goes over into the Weibel domain beyond a critical value of the mass ratio. Only systems already located in the oblique domain for continue to do so as we alter from to . Of course, we speak here only about the lower part of the graph that corresponds to small values of . The dominant mode depends through critically and in a nontrivial way on both, and on in the upper part of Fig. 6.

We thus find a significantly extended region of the parameter space, namely, the uniform white domain on Fig. 6, where the criterion that the mode hierarchy be unchanged does not restrict the value of the mass ratio. For a system lying in this region, the dominant mode is the same, regardless of whether or . In the lower-right corner (i.e., diluted ultra-relativistic beams), the border is defined by the equality of the oblique mode (see Eq. III.2) with the Buneman one for ,

(13) |

In the lower-left corner (i.e, diluted, weakly relativistic regime), the border is determined by equating the oblique growth rate with the Buneman one, but now for ,

(14) |

The upper-border of the uniform Weibel domain is analytically more intricate. Let us just mention that the particular shape exhibited for arises from the Weibel growth rate which reaches a maximum around this value. Expression (III.2) for this quantity makes it clear that , while . As a consequence, reaches a maximum for an intermediate Lorentz factor , which is easily calculated from Eq. (III.2). Although this value is not exact for close to unity and for , the “Weibel optimum” for stays close to , explaining the “bump” at this location.

## Iv Conclusion

The importance of the electron to ion mass ratio for the realism of PIC simulations has been addressed here from an analytical point of view. Since there cannot be any rigorous similarity theory encompassing this quantity, an attempt has been made to identify a threshold , beyond which a given simulation can no longer be trusted to be physically accurate during the initial exponential growth phase of the instability. This initial wave growth can be addressed by a linearized theory.

Whether it be relevant for astrophysical plasmas or for inertial fusion, many systems investigated through PIC simulations give rise to the growth of waves that can be addressed by an analysis of the linear dispersion relation. The idea is therefore to find the maximum value of , which leaves unchanged the hierarchy of the linearly unstable modes. The condition we propose is necessary but not sufficient: for the system evolution to be preserved, the linear evolution and, more specifically, the type of the fastest growing mode must remain unchanged as we change . However, two similar linear growth phases can eventually result in a different non-linear state.

Because the application of the criterion depends on the linear unstable spectrum, and therefore on the system under scrutiny, we have focused on the generic system formed by a relativistic electron beam passing through a plasma with return current. For a 1D simulation, the competing modes are the two-stream mode and the Buneman mode (see Fig. 2). The criterion of the preserved mode hierarchy does provide a value of , if the spectrum is governed by the two-stream instability for . As a result, for example, the simulation of a 10 times diluted beam with cannot be performed with ions that have a mass below times heavier than that of the electrons. For systems governed by the Buneman instability when , our criterion does not give any upper value of .

The 2/3D case is even more interesting as more modes compete in the linear phase. Here, the Buneman, the oblique and the Weibel instabilities can dominate the linear phase, while the two-stream instability is unimportant for relativistic beam speeds. For and 1/1836, the hierarchy map is plotted on Fig. 3 in terms of the density ratio and the beam Lorentz factor . One can notice how the upper part () does not vary with . This could explain why PIC simulations of collisions between equally dense plasma shells did not show much difference as the mass ratio has been altered Spitkovsky2008 (). The dominant mode is definitely the Weibel (filamentation) one in this region. Things should be different when simulating collisions of shells with a different density, like in DieckmannApJ (); Bessho (), and varying the mass ratio.

The value of in terms of ) is predictable at low , and more involved for (see Fig. 6). A few points can be emphasized at this junction: First, the most sensitive points are the ones located near a border between two modes for . When departs from this value, the border moves, say from mode (A) domain to mode (B) domain, so that (A) increases to the expenses of (B). If the point was initially in the (A) domain, it remains there and the criterion is not binding. But if the point was close to the border, yet in the (B) region, then a slight increase of transfers it to the (A) region. This is why on Fig. 6, the white region of unconstrained always borders .

Second, some alteration of the mode hierarchy are more dramatic than others. Figures 4 & 6 show that 3 kinds of transitions can be triggered when increasing : oblique to Buneman (OB - for diluted beam), oblique to Weibel (OW - high density, weakly relativistic) and Buneman to Weibel (BW - high density, ultra-relativistic). For diluted beams, the OB transition switches the wavelength of the dominant mode from to , resulting in generated structures times smaller. Furthermore, oblique modes generate partially electromagnetic transverse structures whereas the electrostatic Buneman modes do not. More dramatic can be the OW transition as we now switch from a quasi-electrostatic dominant modes to an electromagnetic one. But the BW transition is by far the most powerful as the generated patterns switch from stripes (Buneman) to Filaments (Weibel).

Note however that transitions are smoother than they appear because the switch from one mode regime to another is not immediate when a border is crossed. Suppose we move from domain (A) to (B). As we approach the border, mode (A) keeps growing faster, but mode (B) grows almost as fast, until the growth rates are strictly equal right on the border. There is therefore a zone extending on both side on the line where a proper interpretation of the linear regime needs to account for the growth of (A) plus (B), thus smoothing out the transition.

The present article is a first step towards a systematic search. The method proposed has been applied to a generic beam-plasma system, evidencing non-trivial values of . A similar analysis can be easily conducted varying the set-up: one needs first to evaluate the growth-rate map (the counterpart of Fig. 3) as a function of for the system under scrutiny with . The same plot is then evaluated for the desired value of . If the dominant mode remains the same, then the present criterion is met.

For example, when dealing with the problem of magnetic field amplification and particles acceleration in Supernova Remnants, a typical PIC setup consists in a non-relativistic beam of protons passing through a plasma with a guiding magnetic field Niemiec2008 (); Riquelme (). Due to the magnetization, unstable modes such as the Bell’s ones Bell2004 () enrich the spectrum, and it would be interesting to apply our criterion also to these cases.

###### Acknowledgements.

A Bret acknowledges the financial support of project PAI08-0182-3162 of the Consejeria de Educacion y Ciencia de la Junta de Comunidades de Castilla-La Mancha. ME Dieckmann acknowledges financial support by the Science Foundation Ireland grant 08/RFP/PHY1694 and by Vetenskapsrådet. Thanks are due to Martin Pohl and Jacek Niemiec for useful discussions.## Appendix A 2D and 3D Tensor

Choosing the axis for the flow direction, and , the tensor involved in the dispersion equation (11) is symmetric, and reads,

(15) |

where,

## References

- (1) J. Brainerd, Astrophys. J. 538, 628 (2000)
- (2) E. Waxman, Plasma Phys. Controlled Fusion 48, B137 (2006)
- (3) L. Silva, R. Fonseca, J. Tonge, J. Dawson, W. Mori, and M. Medvedev, Astrophys. J. 596, L121 (2003)
- (4) A. Spitkovsky, Astrophys. J. Lett. 673, L39 (2008)
- (5) P. Chang, A. Spitkovsky, and J. Arons, Astrophys. J. 674, 378 (2008)
- (6) Y. Gallant, M. Hoshino, A. Langdon, J. Arons, and C. Max, Astrophys. J. 391, 73 (1992)
- (7) J. T. Frederiksen, C. B. Hededal, T. Haugbolle, and A. Nordlund, Astrophys. J. 608, L13 (2004)
- (8) M. Dieckmann, B. Eliasson, M. Parviainen, P. Shukla, and A. Ynnerman, Mon. Not. R. Astron. Soc. 367, 865 (2006)
- (9) A. Spitkovsky, Astrophys. J. Lett. 682, L5 (2008)
- (10) M. Hoshino, J. Arons, Y. Gallant, and A. Langdon, Astrophys. J. 390, 454 (1992)
- (11) M. Dieckmann and A. Bret, Astrophys. J. 694, 154 (2009)
- (12) C. H. Jaroschek, H. Lesch, and R. A. Treumann, Astrophys. J. 616, 1065 (2004)
- (13) C. H. Jaroschek, H. Lesch, and R. A. Treumann, Astrophys. J. 618, 822 (2005)
- (14) M. Medvedev, M. Fiore, R. Fonseca, L. Silva, and W. Mori, Astrophys. J. 618, L75 (2005)
- (15) M. Dieckmann, P. Shukla, and L. Drury, Mon. Not. R. Astron. Soc. 367, 1072 (2006)
- (16) M. Karlický, Astrophys. J. 690, 189 (2009)
- (17) M. Karlický and M. Bárta, Nonlinear Processes in Geophysics 16, 525 (2009)
- (18) H. Che, J. Drake, M. Swisdak, and P. Yoon, Phys. Rev. Lett. 102, 145004 (2009)
- (19) T. Intrator, X. Sun, G. Lapenta, L. Dorf, and I. Furno, Nature Phys. 5, 521 (2009)
- (20) M. Tabak, J. Hammer, M. E. Glinsky, W. L. Kruer, S. C. Wilks, J. Woodworth, E. M. Campbell, M. D. Perry, and R. J. Mason, Phys. Plasmas 1, 1626 (1994)
- (21) A. Pukhov and J. Meyer-ter-Vehn, Phys. Rev. Lett. 76, 3975 (1996)
- (22) C. Ren, M. Tzoufras, F. S. Tsung, W. B. Mori, S. Amorini, R. Fonseca, L. O. Silva, J. C. Adam, and A. Heron, Phys. Rev. Lett. 93, 185004 (2004)
- (23) M. Medvedev and A. Spitkovsky, Astrophys. J. 700, 956 (2009)
- (24) M. Honda, J. Meyer-ter-Vehn, and A. Pukhov, Phys. Rev. Lett. 85, 2128 (2000)
- (25) J. Niemiec, M. Pohl, T. Stroman, and K.-I. Nishikawa, Astrophys. J. 684, 1174 (2008)
- (26) S. Martins, R. Fonseca, L. Silva, and W. Mori, ApJL 695, L189 (2009)
- (27) G. Lapenta, J. U. Brackbill, and P. Ricci, Phys. Plasmas 13, 055904 (2006)
- (28) K. Noguchi, C. Tronci, G. Zuccaro, and G. Lapenta, Phys. Plasmas 14, 042308 (2007)
- (29) M. Scholer, I. Shinohara, and S. Matsukiyo, J. Geophys. Res. 108 (2003)
- (30) M. Scholer and S. Matsukiyo, Annales Geophysicae 22, 2345 (2004)
- (31) N. Shimada and M. Hoshino, J. Geophys. Res. 110, A02105 (2005)
- (32) F. Califano, L. Galeotti, and C. Briand, Phys. Plasmas 14, 052306 (2007)
- (33) B. Eliasson and P. K. Shukla, Phys. Rev. Lett. 93, 045001 (2004)
- (34) P. Ricci, J. U. Brackbill, W. Daughton, and G. Lapenta, Physics of Plasmas 11, 4102 (2004)
- (35) J. Guo and Q. M. Lu, Chinese Physics Letters 24, 3199 (2007)
- (36) R. Lee, S. Chapman, and D. R.O., Astrophys. J. 604, 187 (2004)
- (37) C. Brennen, Hydrodynamics of Pumps (Oxford University Press, Oxford, England, 1994)
- (38) G. Birkhoff, Hydrodynamics (Princeton University Press, Princeton, NJ, 1960)
- (39) C. C. Petty, Phys. Plasmas 15, 080501 (2008)
- (40) L. Drury and J. Mendonca, Phys. Plasmas 7, 5148 (2000)
- (41) D. Ryutov, B. Remington, H. Robey, and R. Drake, Phys. Plasmas 8, 1804 (2001)
- (42) S. Gordienko and A. Pukhov, Phys. Plasmas 12, 043109 (2005)
- (43) F. Califano, R. Prandi, F. Pegoraro, and S. V. Bulanov, Phys. Rev. E 58, 7837 (1998)
- (44) E. Buckingham, Phys. Rev. 4, 345 (1914)
- (45) A. Bell, Mon. Not. R. Astron. Soc. 353, 550 (2004)
- (46) A. Bret, Astrophys. J. 699, 990 (2009)
- (47) F. Califano, F. Pegoraro, and S. V. Bulanov, Phys. Rev. E 56, 963 (1997)
- (48) F. Califano, F. Pegoraro, S. V. Bulanov, and A. Mangeney, Phys. Rev. E 57, 7048 (1998)
- (49) A. Bret, M.-C. Firpo, and C. Deutsch, Phys. Rev. Lett. 94, 115002 (2005)
- (50) A. Bret and M. Dieckmann, Phys. Plasmas 15, 012104 (2008)
- (51) M. Medvedev and A. Loeb, Astrophys. J. 526, 697 (1999)
- (52) A. Bret, L. Gremillet, D. Bénisti, and E. Lefebvre, Phys. Rev. Lett. 100, 205008 (2008)
- (53) K. M. Watson, S. A. Bludman, and M. N. Rosenbluth, Phys. Fluids 3, 741 (1960)
- (54) S. A. Bludman, K. M. Watson, and M. N. Rosenbluth, Phys. Fluids 3, 747 (1960)
- (55) Y. B. Faĭnberg, V. D. Shapiro, and V. Shevchenko, Soviet Phys. JETP 30, 528 (1970)
- (56) A. Bret, Comp. Phys. Com. 176, 362 (2007)
- (57) M. E. Dieckmann, P. K. Shukla, and L. Drury, Astrophys. J. 675, 586 (2008)
- (58) N. Bessho and Y. Ohsawa, Phys. Plasmas 6, 3076 (1999)
- (59) M. Riquelme and A. Spitkovsky, Astrophys. J. 694, 626 (2009)