The plasma filamentation instability in one dimension: nonlinear evolution

The plasma filamentation instability or beam-Weibel instability generates magnetic fields and accelerates particles in collisionless astrophysical plasma. This instability has been examined with multi-dimensional particle-in-cell (PIC) simulations, demonstrating the formation of current flux tubes. Such simulations could not model a statistically significant number of filaments. Here, we model with a PIC simulation the filamentation instability that is driven by nonrelativistic, cool electron beams in one spatial dimension at an unprecedented resolution. We show unambiguously that the gradient of the magnetic pressure which develops during the quasi-linear evolution of the filamentation instability, gives rise to an electrostatic field component. The interplay of the magnetic and electrostatic fields results in a wavenumber spectrum of the magnetic field that is a power-law, which has been reported previously for multi-dimensional PIC simulations. The magnetic field power spectrum decreases with the exponent −5.7 and that of the electrostatic field with −3.8, yielding a ratio of 3:2. The electromagnetic fields thermalize the electrons. The electrons develop a velocity distribution in the simulation direction that decreases exponentially at low speeds and faster at high speeds. The filamentation instability can thus not efficiently accelerate electrons to high energies. The filaments develop into a stationary final state. The probability distribution of the filament sizes is a Gumbel distribution. In astrophysical settings, this implies that the long exponential tail of this distribution may lead with a reasonable probability to large current filaments, if the filamentation instability develops in a large enough volume. The coherent magnetic fields of large filaments are required to explain the synchrotron emissions of gamma ray bursts.

The plasma filamentation instability or beam-Weibel instability generates magnetic fields and accelerates particles in collisionless astrophysical plasma. This instability has been examined with multi-dimensional particle-incell (PIC) simulations, demonstrating the formation of current flux tubes. Such simulations could not model a statistically significant number of filaments. Here, we model with a PIC simulation the filamentation instability that is driven by nonrelativistic, cool electron beams in one spatial dimension at an unprecedented resolution. We show unambiguously that the gradient of the magnetic pressure which develops during the quasi-linear evolution of the filamentation instability, gives rise to an electrostatic field component. The interplay of the magnetic and electrostatic fields results in a wavenumber spectrum of the magnetic field that is a power-law, which has been reported previously for multi-dimensional PIC simulations. The magnetic field power spectrum decreases with the exponent −5.7 and that of the electrostatic field with −3.8, yielding a ratio of 3:2. The electromagnetic fields thermalize the electrons. The electrons develop a velocity distribution in the simulation direction that decreases exponentially at low speeds and faster at high speeds. The filamentation instability can thus not efficiently accelerate electrons to high energies. The filaments develop into a stationary final state. The probability distribution of the filament sizes is a Gumbel distribution. In astrophysical settings, this implies that the long exponential tail of this distribution may lead with a reasonable probability to large current filaments, if the filamentation instability develops in a large enough volume. The coherent magnetic fields of large filaments are required to explain the synchrotron emissions of gamma ray bursts.

Introduction
Plasma instabilities, which are involved when nonthermal plasmas relax on timescales that are short compared to the diffusion times [1,2], have recently attracted attention as a means to generate strong magnetic fields in various astrophysical environments, e.g. in the jets of gamma ray bursts [3]- [5] and on cosmological scales [6]. In a simplified system, in which the nonthermal distribution is due to two cool and spatially homogeneous beams of electrons that interpenetrate with a relative flow velocity vector v b , we find three classes of linear instabilities. The two-stream instability results in electrostatic waves with a wavevector k v b [7,8]. The mixed modes are partially electromagnetic and their k is oblique to v b [9]. The filamentation instability leads to the aperiodic growth of electromagnetic waves with k ⊥ v b [10]. The filamentation instability is important, if the modulus v b of the flow velocity vector v b is not small compared to the speed of light c in vacuum, when both plasma beams have a comparable density [11] and if no ambient magnetic field is initially present [12]- [14]. The filamentation instability typically converts about 10% of the kinetic energy into magnetic field energy, as particle-in-cell (PIC) simulations show [15,16].
Recently, a PIC simulation of unmagnetized counter-propagating electron beams has been performed [17]. By modelling only the two spatial dimensions orthogonal to v b , all modes but the filamentation mode have been excluded. However, this restriction has enabled a simulation box size that resolved a large number of current filaments from which some statistical properties could be extracted. The wavenumber spectrum could be approximated by power-laws during the initial growth of the filamentation instability and during the nonlinear stage, when the filaments interact and merge [18]. While the power-law distribution of the wavenumber spectrum during the nonlinear stage of the filamentation instability could originate from a preferential attachment of the filaments [17,19], the origin of the power-law during the initial stage has not been obvious. In principle, it could be due to a phase change [17] or self-organization, as we find it in ferromagnetic and other related systems [20]- [22].
The filamentation modes in a PIC simulation grow from noise that does not show a powerlaw |k| spectrum [23] to their saturation at a growth rate, which is also not a power-law function of |k|. It is thus unclear, how the transition takes place from the initial, linear growth phase of the filamentation instability to a system that displays a power-law |k| spectrum. The identification of 3 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT this mechanism would be difficult in the two-dimensional (2D) simulation geometry employed in [17], because the linear growth and nonlinear filament merging processes occur simultaneously. In this work, we thus resort to a 1D spatial geometry, which permits us to extract the mechanism that is responsible for the power-law generation. It also allows us to examine the statistical properties of the current distribution when the filamentation instability is saturated, since the filament merging is practically inhibited. The large simulation box furthermore captures a statistically significant number of filaments.
Section 2 introduces the PIC simulation method and the initial conditions. Section 3 examines the time-evolution of the filamentation instability for modes with k ⊥ v b . We confirm previous observations that the filamentation instability develops an electrostatic component, once the magnetic fields grow to a large amplitude. The 1D simulation geometry unambiguously reveals the spatial correlation of the magnetic pressure gradient with the electrostatic field. The latter develops in wavenumber space at large k = |k|. The simulation exhibits that the k-spectrum of the magnetic field develops into a power-law function once the electrostatic component grows to a significant amplitude. Once the filamentation instability saturates, the evolution of the electron distribution is locked and the current forms a system of quasi-stationary domains. This final state is, however, stable only in one space dimension [15]- [17]. We analyse the size distribution and correlation of the domains, and we propose that the domain sizes follow a Gumbel distribution [24]. In section 4, we discuss the relevance of our findings with respect to models of interacting current filaments that are employed, for example, in the context of gamma ray burst jets [3,18,25]. The typically invoked interaction of the current filaments through magnetic fields may have to be expanded by the inclusion of electrostatic forces. This is the case even in the system considered here and in [17], in which linear electrostatic instabilities are excluded. The latter can lead to electrostatic waves with a k-spectrum that resembles white noise [26] and to phase space holes [27]- [29].

The plasma equations and the PIC method
We consider two counter-propagating identical electron beams. The beams move along the z-direction with the velocity modulus v b = |v b |. Thus, in the chosen reference frame, the total current vanishes. Particle collisions are neglected, which is justified in many astrophysical environments [4], and the system thus evolves under the influence of the Vlasov-Maxwell set of equations. With µ 0 and 0 being the vacuum permeability and permittivity, with the elementary charge e and electron mass m e , the governing equations are We solve equations (1)-(6) with a PIC simulation code that is based on the virtual particle method [30]. In contrast to electromagnetic Vlasov simulation codes [31]- [33] that solve equation (6), PIC codes resort to the method of characteristics. The continuous distribution f(x, v, t) is replaced by an ensemble of phase space elements, which we refer to as computational particles (CPs). The momentum p j and the position x j of the CP with the index j are evolved in time with the help of the Lorentz equation of motion dp

This is repeated for all
CPs. Each CP has a mass m c m e and charge q c /e 1, but q c /m c = e/m e . The microcurrent q j v j of the jth CP is interpolated onto the numerical grid. This is repeated for all CPs, which gives the global current J. The electromagnetic fields are updated by J with the help of the equations (1) and (2). The new fields are interpolated to the position of each CP, updating its momentum. This cycle is repeated for a given number of time steps.

The simulation set-up
The 1D simulation box is aligned with the x-direction, defined by the unit vector e x , and it uses periodic boundary conditions. We denote the spatial coordinate by the scalar x. Two electron beams are modelled that move along the z-direction. The charge is compensated by a background of positively charged immobile ions. This system corresponds to spatially homogeneous beams of infinite extent in all directions, which are cut perpendicularly to their flow velocity vector by the simulation box. This excludes all wave modes except the filamentation mode. The resolution of only one perpendicular direction implies, that the developing filamentation waves are planar and have a wavevector k e x , even though the growth rates are isotropic for all k ⊥ v b . However, the plane wave approximation is valid during the initial linear growth phase, as long as interactions between waves (weak turbulence) are negligible.
Each electron beam has a density n e , that corresponds to the plasma frequency ω e = (e 2 n e /m e 0 ) 1/2 = 2π × 10 5 s −1 . The total plasma frequency p = √ 2 ω e . Since the grid cell size in PIC simulations should be comparable to the Debye length λ D = v e / p , where the electron thermal speed v e = (k B T/m e ) 1/2 , we cannot model the cold plasma that is frequently invoked when filament interactions are examined [17,18,25]. We select v e = 5 × 10 6 m s −1 and v b /v e = 18 in our simulation and thermal effects should thus be small. The x-direction is resolved by N g = 3 × 10 4 simulation cells, each with the length x = 0.9λ D . Typically, the waves that are driven by the filamentation instability have wavelengths λ = 2π/k ∼ λ e , where λ e = c/ p is the electron skin depth. The box length is L = N g x = 444 λ e and the lowest wavenumber k = 2 π/L. The filamentation instability is thus well-resolved. Each beam is represented by N e = 6 × 10 6 particles, or 200 particles per cell. We let the system evolve in time for a simulation duration T S = 533/ p that is subdivided into a total of 6 × 10 4 time steps t . Initially, E = 0 and B = 0. The simulation noise due to the finite number of CPs introduces broadband thermal noise [23], which provides the initial amplitude for the wave growth.

Field evolution
We analyse the kinetic energy density E k,  Figure 1 shows E k,j , E E and E B during an early simulation time. The energy densities are normalized to E 0 = E k,1 + E k,2 at t = 0. The mean 0 E 2 z , which is the electric field energy density associated with the linear filamentation instability [9], remains two orders of magnitude below that of E E , while the mean values of 0 E 2 y and B 2 z /2µ 0 stay at noise levels (not shown here). At t p ≈ 5, electromagnetic filamentation modes are growing and E B (t), which grows at twice the rate of the wave amplitude, can be approximated by the exponential exp ( i t) during the time 5 < t p < 20 with i = 0.45 p . The peak growth rate of the cold filamentation instability is Max /ω e ≈ v b /c if both beams have equal densities and if ≈ 1 [9]. For our v b = 0.3 c, the growth rate is Max is an integral over the energy contributions of all wave modes, its growth rate is close to i . This is because the 6 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT fastest growing modes dominate the k-spectrum and thus E B . After t p ≈ 20, we observe in figure 1(b) a growth of E E . Electron beams flowing in the z-direction cannot drive an electrostatic two-stream instability in the x-direction. This, together with the observation from figure 1 that E E grows at twice the rate of E B , implies that the E x field is driven by the filamentation instability.
A simple 1D fluid model reveals the cause. It includes the magnetic pressure P b = B 2 y /2µ 0 , which acts along the x-direction. The pressure gradient is dP B /dx = µ −1 0 B y (dB y /dx). In the absence of ions, the charge density is related to the electron density n e (x). The linearized density variation δn e (x, t) evolves in time according to The linear filamentation instability driven by cold beams results in an aperiodic growth of where k 0 is a lower bound for the unstable wave spectrum [10]. A spatial variation of B y arises from the random simulation noise, which yields the initial amplitude B y (x, t 0 ) at the time t 0 = 0. The filamentation instability is characterized by a high, slowly changing growth rate over a broad k-interval. The superposition of many unstable waves allows for an independent growth of fluctuations of B y at different locations. Consider two neighbouring positions x 0 and x 1 = x 0 + dx that have an initial (noise) amplitude B y (x 0 , t 0 ) and B y (x 1 , t 0 ). These fields grow at the same rate and dB We find for the pressure gradient the growth B y dB y /dx ∼ exp ( i t). The pressure introduces an aperiodically growing charge density modulation δn e (x) ∼ exp ( i t) and thus an electric restoring force E x (x) ∼ exp ( i t). The electric field energy density grows as E E ∼ exp (2 i t), in line with figure 1(b). At t p ≈ 28, E B and E E saturate and remain approximately constant on the logarithmic scale until t p = 80. Eventually an almost steady state distribution is reached and after t p = 400 the field energies change by less than 10% (not shown). The electron beams lose about 10% of their initial energy to the growing electromagnetic fields, in line with previous results [15,17]. A Fourier transform over the x-direction in the form of is applied to E x and B y , the power spectrum F 2 (k k , t) = F(k k , t) · F * (k k , t) is computed and shown in figure 2. Initially, only B y grows at a k ≈ p /c, in line with figure 1. Eventually, the E x component grows in a k-interval that coincides with the high-k tail of the B y -spectrum. This overlap suggests that the E x field is connected to rapid (high-k) variations of B y , while the slow B y variations do not give rise to an E x . This can also be explained with the help of equation (7). The sharpest gradients in B y result also in rapid variations in the pressure force and, thus, in growing E x fields. Once E x has reached a high power, both wave spectra shift towards lower k. At all times we observe from figure 2 the correlation between E 2 x (k, t) and the high-k tail of B 2 y (k, t). Figure 2 demonstrates that a steady state is reached after t p ≈ 100, when the energy densities in figure 1 approach an equilibrium.
The real frequency of the waves growing due to the filamentation instability is zero (aperiodic growth). We thus expect that the wave spectrum of B y is dominated by stationary structures after the instability saturates. If E x is connected to the filamentation instability, we also expect that the electrostatic structures are time-independent. We apply the Fourier transform equation (8) over space and time. We transform E x and B y over the full simulation box and over the time interval 400 < t p < 533 and compute the power spectrum. The spectra are displayed in figure 3. Both field components have their peak power at = 0 and are thus time-stationary. The peak power in E 2 x (k, = 0) is two orders of magnitude smaller than that of B 2 y (k, = 0). The spectrum in figure 3(a) further reveals noise propagating on the electromagnetic wave branch at low k and | | ≈ p as well as noise from the Langmuir wave, which has a weak electromagnetic component in the presence of B y = 0. The frequency spread of the stationary structure in figure 3(a) at = 0 arises from the finite sampling time, while the broadband wave component in figure 3(b) is typical of thermal noise [23].
The power spectra E 2 x (k, t) and B 2 y (k, t) from figure 2 are each integrated over the time interval 400 < t p < 533, during which the fields are quasi-stationary, and the spectra are compared in figure 4. Both field components show spectra that are well-approximated by powerlaws over a wide range of wavenumbers. At high k, the structures become predominantly electrostatic. The power spectrum of B 2 y in figure 4 is steeper than that of P ⊥ (k) = B 2 y (k) + B 2 z (k) given in [17]. In [17], the 2D power spectrum has been integrated over the azimuth angle, which corresponds to a multiplication of the spectrum by k. The slope of B 2 y (k) above kc/ p = 2 is equal to 3/2 times that of E 2 x (k) in the same interval. Note, that we cannot obtain this steady state result from equation (7) alone. By equating the squared electrostatic restoring force and the squared magnetic pressure gradient in k-space, we would get the relation E 2 x (k) ∼ B 2 y (k)(kB y [k]) 2 = k 2 B 4 y (k) at high k.

Electron distribution
In what follows, we concentrate on the electron phase-space distribution f(x, v z ), since it is in this plane that the plasma structures form that develop in response to the filamentation instability. The distributions at four different simulation times are illustrated in figure 5. Only the beam with v z > 0 is shown. Since both beams are identical, the second beam has a similar distribution. At the time t p = 0 the beam is spatially homogeneous and has the speed modulus v b = 18 v e . During the initial growth phase 5 < t p < 20 of the instability, the electrons are redistributed along the x-direction, as the snapshot at t p = 18 shows. However, the phase space density modulation is weak and no significant E x is yet present.At the time t p = 28, when the filamentation instability in figure 1(a) has started to saturate, the electrons form phase space islands that are separated by depleted intervals. This space charge leads to the strong electrostatic field component in figure 1(b). The phase space islands still evolve in time. Eventually, a time-stationary phase-space distribution is reached at the simulation's end, shown in figure 5(d). The evolving phase-space distribution in a sub-interval of the simulation box is also shown in the supplementary movie. The movie demonstrates that the beam evolves from the initial, spatially homogeneous state to a quasi-stationary final state. Figure 6 confirms that the final E x and B y field components are directly related to the electron distribution f(x, v z ). The peak value B max of B y in the simulation box is significant, since it results in an electron gyrofrequency ω c = eB max /m e = p /5. Figure 6(a) reveals that E x = 0 when dB y /dx = 0 or when B y = 0, which results in a vanishing magnetic pressure gradient in equation (7). The k-spectrum of the E x -component must thus peak at larger k than that of B y , which is in line with the figures 2 and 4.
Since only E x and B y reach a significant amplitude, the E × B force points along the z-direction. The electrons at the equilibrium point x 0 p /c ≈ 2.5 in figure 6 experience no electromagnetic fields and move freely, except for their interaction with the weak E z component of the filamentation instability [9]. As we move away from x 0 to lower or to larger x, the moduli of E x and B y grow linearly. Since E x B y > 0, as we move away from x 0 , the v z -modulation is symmetric relative to x 0 and it is identical for both beams. The v × B force deflects the electrons of the beam with v b > 0 towards x 0 and those of the second beam are deflected away from x 0 . Each beam is thus pinched or diluted, depending on the v z -direction and on the sign of the magnetic field gradient. Figure 5 illustrates the structure formation within f(x, v z ). The v y -component is unaffected by the filamentation instability in the 1D simulation geometry employed here, since B x ≡ 0. While no clear structures develop in the f(x, v x )-plane (not shown here), the electrons are accelerated due to their large v z and the strong, spatially varying B y field. Figure 7 displays the final, box-averaged electron distributions. The electron v y distribution is the initial Maxwellian distribution. The v z distribution has been broadened due to the filamentation instability, but most electrons still move at speeds ≈ v b . The f(v x ) distribution shows an exponential slope up to a value v x ≈ v b , above which it is dropping faster than exponential. Figure 8(a) depicts the electron phase-space distribution at t p = 533. Both electron beams form similar filaments. The density maxima of each beam spatially coincide with the minima of the other beam. The current component J z calculated with equation (5) is displayed in figure 8(b) and it shows an ordered sequence of current domains with a practically constant value of J z . As we go from one domain to the next, the current changes its sign. The width of the current domains corresponds to the distance between two consecutive zero-crossings in figure 8(b).   This width is not constant. Thence, there is a distribution function for the probability of finding a domain with a certain width. However, with a simulation box length of L = 444 λ e , the total number of current domains is only 234 at t p = 533. This number is not sufficient to obtain a statistically significant distribution function.
We thus expand our simulation box length by a factor 15 and reduce the simulation time, while we keep all other plasma and simulation parameters unchanged. We measure J z at the time t p = 390. We exclude the zero-crossings with a separation of less than 0.3λ e , which removes the small scale (noise) fluctuations around J z ≈ 0. A trend in the size distribution can be extracted, by plotting the size of the (j + 1)th domain versus that of the jth domain. Figure 9 reveals a tendency of the filaments to cluster according to their size. Neighbouring domains preferably have a similar size. Most large domains with a size > 4λ e are found next to other large domains. However, this is not an exclusive rule because we find in figure 9 various domains j with a size above 4, which are next to domains (j + 1) with a size of around one.
The large simulation, which resolves ≈6.6 × 10 3 λ e , contains about 3500 current domains. The size of the current domains in figure 9 can reach values of up to 5.5λ e . The range of possible filament sizes 0.3λ e < S < 5.5λ e is subdivided into intervals with a constant width S = 0.0444λ e . The number of filaments that fall into a given size interval is counted and the resulting count rate (probability distribution) is displayed in figure 10. The count rates follow a curve with a single maximum that can be approximated by the Gumbel distribution [24], as for example discussed in [21]. A distribution P D = 210 × exp (−u − exp [u]) with u = 1.45 × (x/λ e − 1.65) is fitted, which reproduces the single maximum and the slopes of the measured distribution at low and large sizes.

Discussion
Here, we have examined the growth and saturation of the filamentation instability in one spatial dimension. This reduced dimensionality has firstly allowed us to suppress the filament merging [18], which competes with the wave growth. Secondly, in 1D denoted by x, we can separate an electrostatic E x -component from the spatially varying strong B y and weak E z electromagnetic components of the filamentation instability [9]. This has allowed for an in-depth analysis of the electrostatic E x component, which is developing during the quasi-linear growth phase of the filamentation instability. We have also found that the electron velocity distribution along the wavevector k decreases exponentially in a 1D geometry and that the final distribution of current domain sizes follows a Gumbel probability distribution [24].
The low plasma temperature (thermal pressure) has led to a high signal-to-noise ratio for E x [23] and to a dominance of the magnetic pressure over the thermal pressure, facilitating its investigation during the quasi-linear growth phase, when the electrons have not yet been heated. We could thus unambiguously identify the relation between the electrostatic and the magnetic energy densities as well as the k-spectrum and the frequency spectrum of the electrostatic component. Since the E x -component is driven by the spatial gradient of the aperiodically growing magnetic pressure P b = B 2 y (x, t)/2µ 0 ≈ B 2 y (x) exp ( i t)/2µ 0 , it goes to zero when either B y (x) or dB y (x)/dx vanish. Consequently the characteristic wavenumber of the E x -component is twice that of B y and the real part of its frequency vanishes. We have confirmed this here.
At the simulation's end, the k-interval that contained the strong E x field has coincided with the high-k tail of B y . In the overlap interval, the power spectra of both components could be approximated by power-law distributions. Previously only power-law distributed wavenumber spectra for the magnetic fields have been reported, without specifying their physical origin [16,17]. The power-law cannot be explained with the filamentation instability alone, which shows an almost constant growth rate at high k, because the initial amplitude distribution given by noise is not a power-law [23]. The overlap of a strong E 2 x (k) signal with the power-law distributed high-k tail of B 2 y (k) suggests that the E x component results in the power-law slope of B 2 y (k) at large k. The slopes are such that B 2 y (k) ∼ k −5.7 and E 2 x (k) ∼ k −3.8 , giving an exponent ratio of 3:2. For sufficiently high k, the quasi-linear filamentation instability thus becomes electrostatic. The functional dependence of E x and B y on k during the nonlinear, time-stationary equilibrium of the filamentation instability does, however, not simply arise from an equilibrium between the magnetic pressure gradient and the electrostatic restoring force, which led to the growth of the E x component during the quasi-linear phase. If we assume a drop E x (k) ∼ k −1.9 and thus E 2 We attribute this to the strong electron heating in the x direction, which gives rise to a significant thermal pressure contribution.
Once the final equilibrium is reached, the electrons show an exponentially decreasing velocity distribution f(|v x |) up to the initial beam velocity v b and an even faster drop above v b , in line with that in 2D simulations [17]. It confirms that the filamentation instability is not an efficient electron accelerator in one or two spatial dimensions. Thus, electrostatic instabilities [17,34] rather than the filamentation instability are probably responsible for the power-law distributed electron energy spectrum reported in [35]. This is because the different beam density ratios in [35] have favored electrostatic over electromagnetic instabilities [11].
The large simulation box has revealed the probability distribution for the current filament size (magnetic domains). The size distribution is not a Gaussian, which would have suggested a random and unconstrained current accumulation process. The distribution is best fitted to an extended Gumbel distribution F(x) ∼ exp {a(y − exp [y])} with y = 1.45(S j /λ e − 1.65) and importantly a = 1. The case a = 1 corresponds to the original Gumbel distribution for the maximum value in a random draw from a set of statistically independent values, each with an exponential distribution [24]. Chapman et al [21] further discusses the extended Gumbel distribution and its relation to a range of physical effects such as the intensity distributions of sunspots.
In our simulation, we have N e particles and a box averaged mean value J z = 0. Since both beams are initially cool, each particle contributes an initial microcurrentJ z ≈ v b q c z. Due to the initial random noise and the filamentation instability that result in a broad k-spectrum of growing waves, these particles will cluster together and form randomly distributed current filaments. Each of these filaments j will try to maximize its size S j , until all particles are distributed. Since each filament consists to a good approximation either of electrons with v ≈ +v b or v ≈ −v b the total current, which in analogy to the magnetic domain structure we call the magnetization M j , is an extremal value. When the time-stationary, final equilibrium is reached, each filament has a spatially constant macroscopic J z (x) and a magnetization M j ∼ S j J z . Each M j represents the maximum magnetization that can be reached locally through an accumulation of the microcurrents of the individual particles. The distribution of these maxima will then obey the Gumbel distribution with a = 1.
With respect to astrophysical plasmas, we find that the filamentation instability can result in flux tubes that interact through both, electrostatic and electromagnetic forces. A previous model [25] of flux tubes in the jets of gamma ray bursts, which result from the filamentation instability [3] and interact only through magnetic forces, may have to be expanded. The (Gumbel) distribution P(S j ) of the current filaments with size S j , which is apparently generated by the filamentation instability, has a long tail. If the volume, in which the filamentation takes place, is large compared to the electron skin depth λ e , we may find many large current flux tubes. Due to their large perimeter, these flux tubes will rapidly grow by coalescence with neighboring filaments. The latter are typically also large as the simulation in this work evidences. The comparatively small box volumes of the 3D PIC simulations [15,16,35] may underestimate the maximum size of the current filaments, on which a standard model for the gamma ray burst emission is based [3]. Future work will address with PIC and hybrid simulations [36,37] if ion beam filamentation also results in a Gumbel distribution. This would result in even larger flux tubes and associated coherent magnetic fields.