Aspects of self-similar current distributions resulting from the plasma filamentation instability

Colliding plasmas can form current filaments that are magnetically confined and interact through electromagnetic fields during the nonlinear evolution of this filamentation instability. A nonrelativistic and a relativistic electron flow are examined. Two-dimensional (2D) particle-in-cell (PIC) simulations evolve the instability in a plane orthogonal to the flow vector and confirm that the current filaments move, merge through magnetic reconnection and evolve into current sheets and large flux tubes. The current filaments overlap over limited spatial intervals. Electrons accelerate in the overlap region and their final energy distribution decreases faster than exponential. The spatial power spectrum of the filaments in the flow-aligned current component can be approximated by a power-law during the linear growth phase. This may reflect a phase transition. The power spectrum of the current component perpendicular to the flow direction shows a power-law also during the nonlinear phase, possibly due to preferential attachment. The power-law distributed power spectra evidence self-similarity over a limited scale size and the wavenumber of the maximum of the spatial power spectrum of the filament distribution decreases linearly in time. Power-law correlations of velocity fields, which could be connected to the current filaments, may imply super-diffusion.


3
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT resolution is one purpose of this study. To achieve the higher resolution, we simplify the system. We assume infinitely long electron plasma slabs with antiparallel velocity vectors. The simulation plane cuts the velocity vectors orthogonally, as in the [31,35,42]. We thus neglect important contributions from electrostatic instabilities, phase space holes [13,14] and the ion dynamics [51,52]. The wavenumber interval holding the power spectra of the filaments is bounded in our simulations by statistical noise at high wavenumbers, which is related to thermal noise [53]. The large box excludes finite-box effects.
Other 2D PIC simulations, e.g., that in [29], have examined a plane that contains the initial plasma velocity vector. The study in [30] has investigated statistical properties of the filaments with the help of a particle-fluid hybrid code, showing that the number N f of filaments is decreasing in time t as N f ∼ t −0.9 , if only the electron dynamics is taken into account. In [31], a model for the interacting current filaments is specified, taking into account filament mergers. The equation of motion of the filaments is determined and is used to predict their collective dynamics. The evolution of the resulting magnetic fields and their correlation length λ B is followed in configuration space by a PIC simulation and it has increased like λ B ∼ t 0.7 and λ B ∼ t 1.1 for the nonrelativistic and the relativistic case, respectively.
Here, we examine the current distribution that is provided by the PIC simulations. We confirm that the current filaments merge by magnetic field reconnection through current annihilation [54] and that their size increases in time [42]. The filaments are sheaths, due to a weak magnetic pinch. In the considered nonrelativistic and relativistic cases, the wavenumber k M of the maximum of the power spectrum of the currents evolves as k M ∼ t −1 . The term power spectrum always refers to the spatial power spectrum. We test the developing current filaments for self-similarity [55,56]. Self-similarity implies the absence of any characteristic scale in the system. A function f(x) is scale free if f(x) = λ −D f(λx). A power spectrum P(k) ∼ k D fulfills this requirement. We find that the current distribution, which is developing out of the instability, is not necessarily self-similar. The current component parallel to the plasma flow velocity vector shows a power-law slope in its power spectrum only during the linear growth phase of the filamentation instability. This self-similarity could arise from a phase transition [57,58]. On the other hand, the current components perpendicular to the initial plasma flow velocity vector show a power-law slope in their power spectrum during the nonlinear evolution of the instability. Here, the self-similarity could be due to preferential attachment [56], since the filament perimeter influences the merger probability. The wavenumber range of the power-law distributions is not sufficient to demonstrate unambigously such a spectrum, but we think that it allows for such an interpretation. The self-similar current distribution parallel to the initial flow direction also induces a self-similar magnetic field. This would agree with the observed self-similarity of the perpendicular magnetic field in [27], if the ion filamentation instability would still be in its linear growth phase. This is probably the case, since the magnetic fields have been measured close to the initial collision boundary.
The filaments in our simulations can be considered to be a warm one-component plasma outside the intervals, in which the filaments overlap. This may imply that the self-similarity of the current filaments also holds for their velocity fields. Velocity fields in hydrodynamics that show power-law correlations evidence super-diffusive systems with short-range attractive and long-range repelling forces [59]. The filamentation instability may constitute a plasma system exposing such a behaviour.
The plasma thermalization by the filamentation instability yields, for the considered plasma parameters, an electron energy distribution that is dropping faster than exponential for increasing 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT energies. This suggests that the power-law distribution obtained for the electron energies in [28] may not be a consequence of the filamentation instability but of the electrostatic fields [60]. In contrast, the rapidly dropping electron distribution in [61] is in good agreement with our final electron distribution. The simulations in the [28,61] employ different beam density ratios and the one in [61] is close to that we consider here.
The present paper is structured as follows. In section 2, we discuss our simulation method and the chosen initial conditions. In section 3, we present our simulation results and we discuss their implications in section 4.

Code description and initial conditions
Our PIC simulation method [62] solves the Maxwell equations for the electric and magnetic fields and the relativistic Lorentz equation for the plasma particles in two spatial and three momentum dimensions.
The code fulfils ∇ · E = ρ/ 0 and ∇ · B = 0 to round-off precision. The electromagnetic fields are solved on a 2D grid with periodic boundary conditions and a constant cell size x in the x-and y-directions. Each computational particle with index cp, charge q cp and mass m cp moves on a continuous trajectory and represents a number of physical particles that is determined by the available computer resources. In this study, we consider only electrons. In what follows, we refer to a computational particle (electron) as CP. While the charge q cp and the mass m cp of a CP can be much larger than the elementary charge −e and mass m e of a physical electron, q cp /m cp = −e/m e . The electromagnetic fields are interpolated to the position of the respective CP in order to update its momentum The position of the CP is updated with the help of its velocity v cp and the simulation time step t . The current contribution q cp v cp of each CP is interpolated back on to the grid. The summed contributions of all CPs give the global current, with which the electromagnetic fields are updated.
We model two oppositely propagating electron species, each with the plasma frequency ω pe = 10 5 2πs −1 , that move at the speed modulus v b1 = 0.26 c in simulation 1 (index 1) and v b2 = c/ √ 2 in simulation 2 (index 2) along the z-direction. The interpenetration speeds are thus v c1 = 2 v b1 /(1 + v 2 b1 /c 2 ) = c/2 and v c2 ≈ 0.94 c, respectively. The v c1 is similar to the speed 0.6 c of the nonrelativistic case discussed in [31]. The speed v c2 is less than the relativistic case considered in [31] with its = 10 . In the rest frame of the respective electron beams, the thermal electron speeds v t 1,t 2 = (K B T 1,2 /m e ) 1/2 fulfil v b 1 /v t 1 = v b 2 /v t 2 = 5. The Boltzmann constant and temperatures are K B and T 1,2 . The minimum speed, at which the filamentation instability outruns the two-stream instability, is comparable to v c1 , provided that no guiding 5 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT magnetic field is present [36,38]. Speeds much higher than v c2 give rise to a rapid growth of the mixed mode instability [36,39]. We resolve the spatial grid by 1400 2 rectangular simulation cells, each with x = λ s /13.5, where λ s = c/ω p and ω p = √ 2ω pe are the electron skin depth and total plasma frequency. Each electron beam is resolved by 160 particles per cell. Both simulations cover the time T sim = 118/ω p , which they resolve by 4500 time steps t . During this time, ions would probably participate in the filamentation, introducing a second time and spatial scale. We deliberately neglect their contribution in order to separate the electron dynamics from the ion dynamics.

Simulation results
We must verify that the simulation results are robust with respect to the grid resolution. For this purpose, we perform a set of small-scale test simulations that resolve a rectangular domain with the size 19 λ s and that represent each electron beam by about N 1 = 10 7 CPs. One simulation resolves the domain by 256 2 cells, giving the cell size x = λ s /13.5, and it represents each electron beam by 160 CPs per cell, which is identical to that we employ in our main simulations. The evolution of the electric field energy in this simulation is compared to that of simulations that resolve the sub-domain by 64 2 , 128 2 and 768 2 simulation cells. We introduce the electric field energy as The volume of each cell is V and i and j denote the cell indices. The total electron kinetic energy is which is summed over both electron species, i.e. N t = 2 N 1 . Figure 1 shows The curves E E (t) disagree for the simulation with the lowest resolution, for that with the intermediate resolution and for both high resolution simulations. The match of the curves in both high-resolution simulations indicates that they are reliable. The grid resolution x = λ s /13.5, we select for the main simulations is comparable to the resolution x = λ s /10 used in [31].
In what follows we examine the current densities measured at each grid cell in the simulations. The current component C z is the component parallel to the initial flow velocity vector and it determines the magnetic fields in the simulation plane. The current components C x , C y in the simulation plane are connected to the magnetic field component parallel to the flow velocity vector. Since our simulation conditions are spatially isotropic, we expect no preferential direction of the currents in the simulation plane and we can introduce the perpendicular current C x + i C y , where i is the imaginary unit. Initially, the spatially uniform current distribution of C z implies that the currents due to both slabs annihilate each other. Statistical fluctuations in the currents should result in a single humped distribution. The filamentation instability channels the initial electron flow and a double-humped distribution will develop. This is confirmed by figure 2. The current distribution is symmetric to C z = 0 in both simulations.
The electron velocities are initialized with identical random number sequences in both simulations. Since the grid resolution and the number of particles per cell are identical, the same random numbers are used for each particle, which is then placed in the same cell in both simulations. In each respective simulation, however, the initialization of the velocities is random and no systematic spatial correlations are introduced. The scaling with the thermal speed and the Lorentz transformation has apparently not altered the initial (noise) current distributions out of which the filamentation instability develops, because the current distribution  The distribution log 10 P(C z , t) and for C z > 0 for the slow beam simulation (panel a) and for the fast beam simulation (panel b). Both panels show an initially single humped noise distribution. As the filamentation instability evolves, a peak at higher current densities C z develops. The distribution is almost symmetric to C z = 0 and, after the onset of the instability, P(C z , t) has two maxima.  are clearly visible. The sheaths are closed and encircle areas with a low current density. This topology indicates a different time evolution of the current channels than a merger of small circular current filaments to larger ones. That the small current filaments merge, however, is demonstrated by the figure 7, which displays snapshots of the mean electron momentum p z /m e c at four consecutive simulation times. The figure 7(a) shows the current distribution, when the fast beam simulation saturates. A circular current filament is visible in the region 2 < x ω p /c < 4 and 2 < y ω p /c < 4 that is attached to a smaller current filament at y ω p /c ≈ 5. Both filaments are separated by a narrow band with a reduced modulus of p z /m e c . The sequence figures 7(a)-(d) evidences the attachment of two further current filaments along a diagonal between the coordinates (xω p /c, yω p /c) that are (0, 4) and (6, 0) to form a current sheath.
To obtain information about the distribution of the electrons in terms of p z /m e c, we slice the displayed box in figure 7(a) along a line that is parallel to the y-axis. This slice is five cells wide to improve the particle statistics. The electron phase space distribution in figure 8 demonstrates that the electrons have not thermalized but that they form localized phase space islands. The modulus p z /m e c of the electrons is lower in the interior of the islands. The bulk distributions of the electrons of adjacent current filaments, with opposite flow directions, are spatially separated. However, electrons diffuse at y ω p /c ≈ 5 and 1.5, which yields a current annihilation.
We slice the figure 5(a) parallel to the x-axis, integrate the electron distribution over five simulation cells close to y ω p /c ≈ 6 and show the phase space distribution in figure 9. Electron bunches are present that overlap only over a limited space, each of which could be considered to be a warm single-fluid electron plasma. The reduction of p z /m e c in figure 5(b) in the interval x ω p /c > 4 and y ω p /c > 2 is due to a reduction of v z . The thermal spread yields electrons with negative and positive v z , explaining the density distribution in figures 5(c) and (d). Current annihilation in the interval 1 < x ω p /c < 2 in figure 9(b) reduces the average p z momentum magnitude in figure 5(b) compared to figure 4(b). A reduction of the mean speed magnitude | v z |, as well as current annihilation, are responsible for the decreasing C z in figure 2. The overlap regions between the current filaments facilitate magnetic field reconnection [54].
We can extract further statistical information from the current density distributions C z and C x + i C y by Fourier transforming them. Our periodic boundary conditions ensure the absence of aliasing in the wavenumber spectrum. The spatially isotropic current distribution and the high grid resolution and domain size, with N g = 1400 cells in the x-and y-directions, ensure a power spectrum that is azimuthally isotropic in wavenumber space. We obtain the power spectrum of the currents C z and C x + i C y , denoted by the index j, with the help of two steps. Firstly, the spatial Fourier transform where C j (x, y) is defined on the grid cells with the index m along the x-direction and n along the y-direction. The integers k x and k y range from −L + 1 to L, where L = N g /2 is the maximum   wavenumber set by the sampling theorem. We introduce the power spectrum C 2 j (k) for a given time t, which is a function of the integer scalark. We relatek to the integer value I k x ,k y of (k 2 x + k 2 y ) 1/2 . We integrate the power spectrum of the current distribution over the azimuthal angle with the help of a second sum Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 10. The power spectrum P(k, t) of the C z current density component in the fast beam simulation. The colour scale denotes the 10-logarithm of P(k, t).
Overplotted is the curve k ∼ t −1 .
Clearly, at low values ofk, this function will fluctuate due to jumps in the number of contributing cells. We define k = 2πk/N g x , which has no subscript to distinguish it from the k x , k y components in the simulations. With 2π/N g x ≈ 1.8 × 10 −4 m −1 , the conditionk 1 thus yields the requirement kc/ω p 0.06 to obtain a smooth spectrum. We calculate C 2 j (k, t) for the fast beam simulation. The time evolution of the power spectrum of C z is displayed in figure 10. The spectral distribution has a pronounced maximum at k M that is gradually shifting in time towards lower k. The overplotted curve k ∼ t −1 is accurately following this shift. The spectrum is sufficiently far away from the minimum and maximum resolved kc/ω p ≈ 0.06 and 40. A similar evolution towards lower k is revealed for C x + iC y by the figure 11. The overplotted curve is identical to that in figure 10, and the structures are thus correlated in their scale size. Figures 10 and 11 suggest a power-law drop of the power spectrum at large k. This is partially confirmed by the figure 12, which slices the figures 10 and 11 at the times ω p t = 24 and ω p t = 112. During early times, just after the filamentation instability has saturated, the power spectrum of C z shows a power-law slope in the interval 3 < kc/ω p < 15 in figure 12(a). The spectrum of the perpendicular current C x + i C y follows that of C z , but only at high wavenumbers does it show a power-law slope. At the late simulation time, the figure 12(b) suggests that only the power spectrum of the C x + i C y component follows a (harder) power-law. The lines of constant colour of the power spectrum at high k are parallel in figure 11 but not in figure 10, supporting this conclusion.
A similar evolution of the power spectra takes place in the slow beam simulation. We show the time evolution of the power spectrum of C z in the figure 13. The power spectrum is single peaked at the wavenumber k M of the spectral maximum, and it follows the overplotted curve k M ∼ t −1 . The curve k ∼ t −0.7 , reported for the nonrelativistic case in [31], does not constitute a good fit for the power maximum and for our considered plasma parameters. It could, however, fit lines of constant colour at higher k and lower power. The shift towards lower k of the power  Overplotted are the curves k ∼ t −1 and k ∼ t −0.7 . spectrum of C z in the slow beam simulation is followed by the corresponding power spectrum of C x + iC y , as the figure 14 reveals. The power spectrum apparently drops like a power-law for large increasing k. To clarify this, we slice the power spectra in the figures 13 and 14 at the times ω p t = 40 and 112 and we display the curves in figure 15. During early simulation times, when the filamentation instability has not yet saturated, the power spectrum of the C z current component follows a power law in the interval 2 < kc/ω p < 10. The power-law fit is, however, not as precise as in figure 12. The power spectrum of the perpendicular currents C x + i C y follows the curve corresponding to C z . At late simulation times in figures 15(b), the power spectrum of C z shows a hump, while the curve for C x + i C y follows a (harder) power law.
Our simulations show a self-similar current distribution C z parallel to the flow velocity vector, until the time when the filamentation instability saturates. This would probably result in a magnetic field component perpendicular to the initial flow velocity vector that is also self-similar [27]. After the nonlinear saturation, however, the power-law slope in the power spectrum of C z disappears and harder power-law slopes develop in the C x + i C y components instead. The origin of the perpendicular current is, presumably, the Lorentz force due to particles moving in the z-direction that crosses the magnetic field in the x, y plane. A self-similar current distribution of C x + i C y but not C z during the nonlinear stage of the filamentation instability implies, however, that the magnetic field B z component should be self-similar, but not the perpendicular B x + i B y component. We test this assumption in the figure 16. The figures 16(a) and (b) show the distribution of the parallel and perpendicular magnetic fields in the full simulation box of the fast beam simulation. The B z component shows circular and sheath structures, while the B x + i B y components are spaghetti-distributed. Figure 16(c) shows that the spaghetti distribution in figure  16(b) is azimuthally isotropic in k x , k y and therefore has no preferential direction.
To reveal a self-similarity in the magnetic field distributions in figures 16(a) and (b), we integrate the power spectrum, e.g., in figure 16(c), over the orientation angle in k-space as we did  for the currents. We compare the power spectra for the parallel and perpendicular magnetic field in the figure 16(d). The power spectrum of the B z component shows a power-law in the interval 0.4 < kc/ω p < 5 and we can probably associate the self-similarity of the C x + i C y current in figure 12(b) with a self-similarity of the B z magnetic field component. The absence of a powerlaw slope and the hump in the power spectrum of B x + i B y is consistent with the power spectrum of C z in figure 12(b).
The spaghetti distribution in figure 16(b) has magnetic x-points that result in magnetic reconnection. Figure 17 displays one of these x-points. Magnetic reconnection enables the filament merger and it is thus of key importance for the nonlinear evolution of the instability. The x-point in figure 17 is well-resolved. However, a coarser grid is likely to underresolve the x-point. This is probably the reason why the results of the test simulations with the lower resolutions were different in figure 1.
The figures 8 and 9 have shown electron acceleration in the overlap regions of the current filaments. The initial electron energy distributions are symmetric around v z = 0 and can be displayed as a function of the relativistic energy m e c 2 ( (v) − 1). The distributions are integrated over the box intervals considered by the figures 4 and 5. The time evolution of the electron energies in the fast beam simulation is shown in the figure 18 and that in the slow beam simulation is shown in figure 19. Initially, both simulations have a beamed electron distribution. The beam in the fast beam simulation has a narrower relative energy spread than the slow beam simulation. We have kept in both simulations equal ratios v b1 /v t1 = v b2 /v t1 in the reference frames of the respective beams. In the box frame, however, the Lorentz transformation compresses the velocity indicates. The electron energy distributions at the end of the slow and fast beam simulations decrease faster than exponential. This implies that the filamentation instability is probably not an efficient particle accelerator unless, possibly, the current sheaths would be very large.

Discussion
In this study we have examined two case studies related to the filamentation instability. We have modelled with PIC simulations the collision-less interaction of two infinitely long and equally dense electron plasmas. The 2D simulation box has been oriented orthogonally to the plasma flow velocity vector, as in the [31,35]. We have thus excluded the two-stream instability and the mixed mode instability [38,39,63], which would be competing with the filamentation instability in the simulation with the nonrelativistic and relativistic flow, respectively. The saturation of these waves results in electron holes, which couple to the filaments in 3D. We have also neglected the ion dynamics, which is important during the nonlinear evolution of the instability [27,51,52]. Our simulations are thus not necessarily realistic, but the restriction to two dimensions enabled us to resolve a large spatial domain at a high resolution.
We have confirmed previous findings that current filaments form in the plasma. During the linear growth phase of the filamentation instability, the current component aligned with the initial plasma flow velocity vector has shown a power spectrum at high wavenumbers that could be approximated by a power-law distribution. We may understand this self-similarity through an analogy. Each moving electron can be associated with a microcurrent. Initially, the system consists of two cool electron beams. In the simulation frame, the microcurrents have the same modulus. Half the currents point in the positive z-direction and the other half in the negative z-direction. The individual electrons couple through their magnetic fields. The filamentation instability cannot flip the direction of the microcurrents, but the electrons can move in the x, y plane. If we replace the current vectors with spin vectors, we would obtain a model that differs from the 2D Ising model [64] only in how the electrons react to collective forces. The filamentation instability may behave similarly to a ferromagnet close to its critical temperature. The phase transition could yield the observed spatial self-similarity [55,57,58,64]. The power spectrum displayed by the currents in the simulation plane has also resembled a power-law. After the nonlinear saturation, the flow-aligned current components have developed humped power spectra, breaking the self-similarity. The currents in the simulation plane have continued to show self-similarity at short spatial scales. Since larger filaments are more likely to merge with other filaments, this self-similarity may be due to preferential attachment [56].
The PIC simulation in [27] has modelled the relativistic collision of two electron-ion plasmas in 3D. A self-similar magnetic field distribution orthogonal to the initial plasma flow velocity vector has been found. This requires a self-similar current component along the flow velocity vector. The study in [27] agrees with our findings, if the magnetic field distribution has been measured during the linear growth phase of the filamentation instability. In [27], the ion filamentation dominates energetically the electron filamentation. The magnetic field distribution has been measured close to the initial collision boundary, implying a short convection time and Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 19. The electron energy distribution P(E kin ) in the slow beam simulation, where E kin = m e c 2 ( − 1) and E 0 is the relativistic kinetic energy associated with the initial beam speed v b1 . Panel (a) shows the time evolution of P(E kin ) and panel (b) shows the initial (black curve) and the final (blue curve) distibutions. thus wave growth time. The results in [31] and those obtained in this study probably agree, provided the 2D and 3D system lend themselves to a comparison.
During the nonlinear evolution of the instability, the self-similar current distribution in the simulation plane has resulted in a self-similar magnetic field component parallel to the flow velocity vector, while the magnetic field in the simulation plane has developed into a spaghettidistribution. We have found magnetic x-points in the latter distribution. The partial overlap of the sheaths violates the conservation of magnetic flux and triggers reconnection [54] and filament merger. The sheath formation and the development of a current system in the simulation plane are probably a thermal effect. The magnetic field has to overcome a high thermal pressure, in order to pinch the current tubes. Therefore, the filament surface is not minimized to give a circular perimeter and the current filaments are not fully separated in space. The magnetic fields generated by the flow-aligned current are thus not immersed in a vacuum but interact with electrons. This interaction probably results in the perpendicular currents. The analytic models specified in the [31,42] for the interaction between current filaments should be expanded to 3D current distributions, if they were to address hot plasmas. Our electron temperatures (1.5 and 10 keV) are comparable to the 100 keV temperatures for gamma ray burst ejecta [65] and the keV temperatures of microquasar plasmas [21] and thermal effects should thus be important in many astrophysical plasmas.
We have to emphasize that the simulations showed a power-law spectrum of the current distributions and of the parallel magnetic field only over a limited range of wavenumbers, typically about an order of magnitude wide. A power-law fit is thus not unambiguous. The power spectrum at high wavenumbers has been limited by simulation noise. The noise power is proportional to the mass of a computational electron [53] and it can thus be reduced significantly by adding more particles. Future studies can thus expand the power-law slope over significantly larger wavenumber intervals, to test our findings. A clear demonstration of a power-law current distribution is important. It would confirm that the filamentation instability indeed constitutes a self-similar system. Moreover, it has to be examined if self-similarity holds also for the velocity field of the filaments. This is suggested by our simulations that show that the filaments could be considered to be spatially separated beams with a well-defined mean speed. Self-similar velocity fields would imply that the filamentation instability is super-diffusive, giving rise to an attractive force on short distances and to a repelling force on long distances, as is known for hydrodynamic turbulence [59], with potentially observable consequences on all plasma scales. However, while self-similarity may hold also in 3D systems [27], a potentially super-diffusive filamentation instability would be competing with instabilities like the two-stream instability, which can be diffusive [7].
Future study has to assess how a coupled electron-ion dynamics affects the system. A promising first result in this respect is that the 3D PIC simulation in [27] has revealed a self-similar magnetic field for a system of colliding electron-ion plasmas that cascades to longer wavelengths. Our simulations support this observation and have shown that the wavenumber k M , at which the power spectrum of the current filaments and thus the magnetic field peaks, is shifted in time t as k M ∼ t −1 .
The intervals, in which current filaments with opposite mean speeds overlap, are electron acceleration sites. The electron energy distribution in both simulations has decreased faster than exponential. Our final electron distribution drops as fast as the one obtained in 3D PIC simulations of the filamentation instability if both colliding plasma slabs have about the same density [61]. A power-law electron energy distribution with a spectral index −2.7 has been found in the 3D PIC simulation of colliding plasma slabs with differing densities [28], which is identical to that resulting from an electrostatic two-stream instability [60]. This may be the consequence of the quasi-electrostatic mixed mode instabilities that dominate if the plasma slabs have strongly differing densities [38,39]. The final electron distribution may thus serve as an indicator of whether the plasma has been thermalized by electrostatic or electromagnetic instabilities. The electron acceleration has been weak in our simulations. However, the acceleration efficiency increases with the filament size. Large relativistic current filaments and sheaths may therefore be efficient particle accelerators [66]- [68 ].