Kinetic effects in the transverse filamentation instability of pair plasmas

The evolution of the filamentation instability produced by two counter-streaming pair plasmas is studied with particle-in-cell (PIC) simulations in both one (1D) and two (2D) spatial dimensions. Radiation friction effects on particles are taken into account. During the nonlinear stage of the instability, a strong broadening of the particle energy spectrum occurs accompanied by the formation of a peak at twice their initial energy. A simple theory of the peak formation is presented. The presence of radiative losses does not change the dynamics of the instability but affects the structure of the particle spectra.


Introduction
The study of a model problem characterized by two colliding electron-positron plasma clouds at relativistic energies can be relevant to various astrophysical scenarios including the fireball model of Gamma Ray Bursts [1], pulsar wind outflows in Pulsar Wind Nebulae [2], and relativistic jets from Active Galactic Nuclei [3].
In this paper we examine the unstable modes generated by two counter-streaming neutral beams of pair plasmas in the ultra-relativistic regime.In particular, we address kinetic effects, such as particle acceleration, taking place during the nonlinear phase of the instability, and we take radiation friction (RF) effects into account.
We use an electromagnetic, fully relativistic particlein-cell (PIC) code to perform simulations in one (1D) and two (2D) spatial dimensions.In both cases, the transition from the linear to the nonlinear phase is characterized by the coalescence of the current filaments.Subsequently, magnetic field and current filaments reach a quasi-stationary regime, with typical scales of several skin depths.Most particles are magnetically confined inside the current filaments.A group of particles are accelerated at twice their initial momentum, forming a peak in the energy spectrum of particles.A simple model for the formation of the spectral peak is presented.
In the simulations performed including the RF we observe a "cooling" effect in the high energy tail of the distribution function of each species without a substantial modification in the dynamics of instability and in the temporal evolution of fields.

One-dimensional results
The initial configuration consists of two neutral beams of electron-positron pairs which propagate in opposite a e-mail: marta.dangelo@gssi.infn.itb e-mail: luca.fedeli@for.unipi.itdirections (corresponding to p z in the momentum space) and fill the entire simulation box.The system is symmetric, with the populations of the two beams having the same initial density n T /4, where n T is the total initial density, and the same momentum absolute value p 0 .Consistently the initial values of charge and current densities and of the electric and magnetic fields are zero.A very small temperature is introduced to seed the instability.
We performed simulations with different Lorentz factors γ 0 from 1 to 10 3 .Here we describe the case with p 0 /m e c = 200 as it is representative of the most relevant effects observed.The simulation box is aligned along the x-direction and it is divided into 15 000 grid cells of equal length x = 0.01 λ p with λ p = c/ω p the skin depth and ω p = 4π e 2 n T /m e 1/2 .Each of the four plasma species is represented by 3 × 10 6 computational particles (200 particles per cell).The total simulation time is t sim = 1000 T p , where T p = 2π/ω p , and the temporal resolution is t = 0.01 T p .At the beginning of the simulation, microfilaments of opposite current start to develop from the thermal noise of the plasma.Each of these fluctuations produces a magnetic field perturbation, thanks to the tendency of a microfilament to shrink and to the repelling interaction between two opposite fluctuations.This magnetic perturbation acts on the fluctuations enhancing the current inhomogeneities.This positive feedback produces an exponential growth of all the quantities, characterizing the linear phase of the instability.The structure of the current density J z as a function of (x, t) is shown in Fig. 1.The development of the instability can be divided into three phases: a linear phase for t < 50 T p , a transition phase for 50 T p < t < 200 T p and a nonlinear quasi-stationary phase for t > 200 T p .
During the exponential growth of the perturbations, J z assumes a filamentary structure with a very small scale length ( λ p ).At t ∼ 50 T p separate filaments of opposite EPJ Web of Conferences current, having a typical scale close to the electron skin depth, become distinguishable and start to merge.This coalescence characterizes the transition of the instability from the linear to the nonlinear quasi-stationary regime.For t > 200 T p the merging phase of the filaments finishes and the size of each filament is constant, so that the configuration can be described as stationary except for some "vibration" which is observable in Fig. 1.
During the initial phase of the linear regime particles can cross field lines, because the amplitude of the magnetic field is not sufficient to strongly deflect their trajectories.As time goes on, the amplitude of transverse magnetic field increases while the Larmor radius of particles decreases.The linear phase ends when the characteristic scale length of the field becomes comparable to the Larmor radius and particles become magnetically trapped into filaments.The analysis of the particle phase space gives an additional insight into the confinement dynamics.
The system presents a two-fluid symmetry, as shown in Refs.[4] and [5].Therefore, it is possible to describe particle behavior considering only two populations, i.e. f − , composed by electrons with p 0 > 0 and positrons with p 0 < 0, and f + , composed by electrons with p 0 < 0 and positrons with p 0 > 0.
Figure 2 (Top panel) represents the projection of the phase space on the (x, p z ) for electrons belonging to f − population at t = 800 T p .Thanks to the two-fluid symmetry, particles of the f + population show a pattern analogous to Fig. 2 with their spatial distribution in space being complementary to that of the f − population, i.e. corresponding to oppositely directed current filaments.
In the regions of peak density, i.e. in the inner part of each filaments, p z ≈ 0. Outside the filaments there is a small number of electrons which have p z ≈ 2 p 0 = 400.The mechanism of acceleration of these high-energy particles will be discussed in the next Sect.2.1.
In the nonlinear, quasi-stationary regime the spatial distribution of particles may be described in terms of an effective potential as follows.The conserved canonical momentum is z = p z ± a z , where a z = (e/m e c 2 )A z is the dimensionless vector potential and ± refers to the sign of the particle charge.At t = 0 we have a z = 0, so z = ± p 0 for the two beams, respectively.The normalized energies of two electrons belonging to the two populations are Equation ( 1) is similar to the energy equation of a particle moving in a one-dimensional configuration under the presence of an effective potential Figure 2 (Bottom panel) shows P− (x) at t = 800 T p .The asymptotic state of the system may thus be described as a state in which the particles cluster into the minima of the effective potential.

Energy spectra and particle acceleration
For each of the two populations ( f + and f − ), the highenergy particles having p z = ±2 p 0 are localized outside the filaments where most of the particles belonging to the other population are localized, as shown in Fig. 2.This feature is also visible in the energy spectrum, reported in Fig. 3, where a peak is present at twice the initial energy The acceleration mechanism may be described using conservation principles and a single particle approach.Energy conservation for two electrons belonging to the two fluid populations gives the relations for the f + and f − fluids, respectively; C + and C − are constants.Suppose now that there are some positions x 0 where an electron of the f − population has | p z | = 0 while p x is maximum, and an electron of the f + population has | p z | maximum while p x = 0 (this hypothesis is consistent with the simulation results, see Fig. 2 and Ref. [4]).This means that in x 0 the x component of the force acting on f − ← electrons ( p 0 < 0) and the z component of the force acting on f − → electrons ( p 0 > 0) have to be zero.The two integration constants C − and C + can be calculated.Independently of C − and assuming that p + x = 0 at x 0 , p x assumes its maximum value where p − z = 0, that is 02005-p.2where a(x 0 ) = − p 0 (due to conservation of the canonical momentum).Then from Eq. ( 3), we obtain C + = 1 + 4 p 2 0 .Due to the symmetry of the system, the vector potential a(x) assumes the same values for its maxima and minima as a function of x, so there will be another point x 0 such that a(x 0 ) = p 0 .Substituting this value of a(x 0 ) in the expression for z , we obtain | p + z | = 2 p 0 .This procedure can be inverted, starting from x and then going to x 0 , obtaining C − = 1 + 4 p 2 0 .Thus, we have demonstrated that in a one-dimensional geometry the energy of the particles reaches a maximum The formation of a peak in the energy spectrum (Fig. 3) around E max is due to the formation of a caustic analogous to that at the base of the sharpness and intensity of a rainbow: the internally reflected rays cross and cluster to form a caustic sheet (see [6] Chap. 3, p. 127), where classically the intensity is infinite.We can apply this principle to explain the peaks in the spectrum.These particles find themselves in positions where |a z | is maximum.Because of the conservation of the canonical momentum, in these positions we have dp z /dx = ±da z /dx = 0. Thus, these particles all gain the same momentum to first order in their distance from the maximum of |a z |.

Two-dimensional results
We now describe the results of 2D simulations, in which the simulation plane (x, y) is perpendicular to the direction of the beams.The simulation box is divided into 2000 × 2000 square cells.The grid dimensions are L g,x × L g,y = 100 λ p × 100 λ p so x = y = 0.05 λ p .Each of the four species is represented by 2 × 10 8 computational particles.The total simulation time is t sim = 200 T p with t = 0.0325 T p .
For these simulations we have used PICCANTE [7], an open-source, massively parallel particle in cell code.
Figure 4 shows the 2D distribution of current density, i.e.J z , at two different times t = 75 T p (left panel) and t = 200 T p (right panel).As a general trend, we observe the formation of isolated structures with a scale length of several electron skin depths at the earliest time (t = 75 T p ). Later, these structures eventually merge forming a structure ("island") of size close to the numerical box in the saturation regime (t = 200 T p ).The distribution of J z inside the island at t = 200 T p is similar to that observed in 1D: the current peaks near the boundary of the island and has much weaker values well inside the island; locally, small scale filaments where the current changes sign are also observed.The island profiles assumed by J z during the non-linear phase are similar to those observed as asymptotic numerical solutions of 2D Navier-Stokes and magnetohydrodynamic equations [8].
The energy spectrum presents a peak at K ∼ 2γ 0 m e c 2 , as in the 1D case (see Ref. [4]), but it is probably smoothed out over an additional degree of freedom.

Effects of radiation friction
The friction effect of the RF force physically corresponds to the incoherent emission of high-frequency radiation by ultra-relativistic electrons and positrons, see [9].Most of the emitted radiation has a frequency high enough to escape the plasma.From a numerical point of view, it is unfeasible to perform simulations with a spatial resolution high enough to resolve such a small wavelength radiation.Thus, it is assumed that such radiation escapes from the system without re-interacting with other electrons or positrons, and the RF acts as a loss term.To evaluate the amount of energy loss due to RF it is thus necessary to compare simulations with and without the inclusion of the RF force.
We checked that, for a fixed initial energy of the beams, the RF effects increase with the initial density n T .This can be explained by the scaling of the field amplitudes as ω p ∼ n 1/2 T .We show results obtained for the highest value considered n T = 3 × 10 21 cm −3 .RF effects have been found to be quite weak for densities < 10 20 cm −3 .
The system organizes itself in filamentary structures for the current density which have almost the same size and features of the filaments obtained in the non-RF simulation.This behavior can be simply understood by noticing that the EM fields have to grow in order for RF to be important, so that the RF plays little role before the 02005-p.3saturation phase.Moreover, in the ultra-relativistic case the dominant term of RF force (see [10]) is ∝ γ 2 .Thus the RF contribution is strongly increased by the acceleration of some particles to higher energy, which is maximized at the instability saturation stage.Indeed, RF effects are more evident in the particle spectra.
Figure 5 shows the kinetic energy distribution of particles at three different times.We observe that the distribution has a peak, i.e. a "bump on tail" as in the case without RF.However, with RF the energy corresponding to the peak decreases with time, and the peak smooths out.This means that the RF force is much stronger for particles belonging to the high-energy tail of the energy distribution, and its effect is thus to "cool down" such high-energy tail.In particular, RF effects play a much stronger role after the initial development of the FI because of the generation of both strong magnetic fields (which lead to synchrotron emission) and the acceleration of particles to high energy.Similar features are observed in the kinetic energy spectrum of 2D simulations reported in Ref. [4].

Conclusions
In this work we have studied the evolution of the filamentation instability produced by two counter-streaming pair plasmas using PIC simulations in both 1D and 2D spatial dimensions, with and without radiation friction effects.The instability development leads to the acceleration of a group of particles to high energy, forming a spectral peak in correspondence to twice the initial kinetic energy and a sharp cutoff at twice the initial drift momentum in the momentum distribution.A simple model has been outlined to account for the acceleration of such particles during the instability development, using conservation principles and a single particle approach.Radiation friction effects have been found to be significant only for relatively high density (∼ 10 20 cm −3 ) and to affect strongly the particle spectra, cooling down the distribution functions, while the instability development is weakly affected.

Figure 1 .
Figure 1.1D complete evolution of the current density J z as a function of position (x) and time (t) over the whole simulation box.The parameters λ p = c/ω p and time T p = 2π/ω p .

Figure 2 .
Figure 2. 1D simulation: (top panel) phase space contours (x, p z ); (bottom panel) effective potential P− (x) = ( p 0 + a z ) 2 , defined in the Eq.(2).These plots show the behavior of the electrons with p 0 /m e c = 200.All these quantities are plotted at t = 800 T p .The color bar indicates number density.

Figure 3 .
Figure 3. 1D simulation: kinetic energy spectrum of the f → at three different times, t = 0 T p (red line), t = 200 T p (green line) and t = 800 T p (blue line) .The spectrum presents a peak at twice the initial energy K ∼ 2γ 0 m e c 2 both at t = 200 T p and t = 800 T p .

Figure 4 .
Figure 4. 2D simulation: contour of current density J z as a function of (x, y) and for times t = 75 T p (left) corresponding to the end of the linear phase and t = 200 T p (right) corresponding to the saturation regime.

Figure 5 .
Figure 5. 1D simulation: energy with RF included of the f − → population at three different times, t = 200 T p (red line), t = 400 T p (green line) and t = 1000 T p (blue line).