Both the pressure anisotropy-driven Weibel instability and the momentum anisotropy-driven current filamentation instability make a quasi-static magnetic field linearly grow. In some conditions, this growth couples with electrostatic perturbations, and an electrostatic field component growing twice as fast as the magnetic field was noticed since the early numerical simulations of these phenomena. We herein provide an interpretation of this process in terms of the electron density concentration induced by the differential rotation of current filaments around the maxima of the magnetic field. We then discuss how this effect, which is both of second order with respect to the amplitude of the electromagnetic Weibel mode and an ingredient of the linear instability itself, anisotropically couples with fluctuations of the distribution functions associated with the pressure tensor components. The analytical estimates are consistent with nonlinear kinetic simulations performed with both the semi-Lagrangian Vlasov code VLEM and with a reduced multi-stream model for the Vlasov–Maxwell system.

## I. INTRODUCTION

The validity of the linear approximation in nonlinear continuum systems is formally delimited by the negligibility of nonlinear coupling terms. Unfortunately, the quantitative importance of the latter in nonlinear numerical calculations may be not as much evident as mathematics would require. This issue concerns also the response of a marginally stable system to small amplitude perturbations, which one would expect to be described in the framework of a linear eigenmode analysis. Plasma physics is rich in examples of this type. They may be related, e.g., the role of particle trapping in the stability of electrostatic modes vs Landau damping^{1–11} or in connection with the role played by linear modes in the turbulence cascade (see, e.g., Refs. 12–14 and references therein). These examples have stimulated the debate—not solved, yet, it seems—about finding the quantitative boundaries between the notion of “linear” and “nonlinear” mode. They also exemplify the general interest in identifying and characterizing the extent at which the notion of nonlinear coupling enters as an ingredient in kinetic effects, such as particle trapping and wave–particle resonances. In this context, it becomes interesting to single out how collective effects of kinetic nature, for example, related to the notion of particle temperature of pressure, intervene in the coupling between linear fluctuations of electromagnetic (e.m.) and of electrostatic nature. This coupling may have implications also for the channels of energy conversion in kinetic plasma turbulence, in which the so-called pressure–strain interaction, macroscopically related to the anisotropic energy transfer induced by a velocity shear on the internal energy of the plasma,^{15–17} has been identified^{18–22} to play a key role. The coupling between magnetic, electrostatic, and pressure components may be relevant also to electrostatic turbulence met in tokamak plasmas, where the nonlinear coupling between pressure and electrostatic fluctuations can be an ingredient in the effects induced by the Reynolds stress tensor, which can affect the evolution of low frequency turbulence with second-order contributions (see, e.g., Refs. 23 and 24).

An interesting—and in some respects paradigmatic—example is provided by the comparison between the linear analysis of the anisotropy-driven Weibel instability^{25} and the nonlinear numerical integration of its evolution, during what is typically recognized as its “linear” stage. This is the subject of the present work. We consider the nonlinear kinetic evolution of a classical Weibel instability in a one-dimensional (1D) non-uniform geometry, in which electrostatic fluctuations are linearly decoupled from e.m. modes. We focus on the characterization of the electrostatic mode that is observed to grow,^{26} since the beginning of the linear stage, as the second harmonic of the unstable mode—we recall that the (“pure,” or transverse) Weibel mode corresponds to an exponential growth, with rate *γ _{B}* and null oscillation frequency, of the magnetic field component that is orthogonal to the plane in which both the wave-vector and the pressure (or temperature) anisotropy lay.

### A. Context: Former studies on this subject

The exponential growth of the electrostatic field with a rate $ 2 \gamma B$ has been observed and discussed in previous works,^{26–29} both on the anisotropy-driven Weibel instability and on the current filamentation instability^{30–35} (CFI), strictly related to it. The CFI can be indeed made correspond^{36,37} to a limit case of the Weibel instability, where the temperature anisotropy associated with a double-Maxwellian velocity distribution is replaced by two counter-propagating cold electron beams of initial densities *n*_{1} and *n*_{2}.

Taggart *et al.*^{26} were probably the first ones to identify this electrostatic growth as a second-order effect in the kinetic linear analysis of Weibel instabilities, by noting that it results from the nonlinear coupling of the electromagnetic mode with itself. This interpretation was supported by a comparison of analytical estimates of the kinetic terms responsible for the nonlinear coupling, with the numerical nonlinear integration of the instability. In Refs. 27 and 29, this coupling has been further explained as due to the nonlinear growth of the current density in the longitudinal direction, which induces the growth of the electric field via the displacement current in Ampère's law. In Ref. 28, the same phenomenon has been related to the coupling between the Weibel mode and the Langmuir waves, which had been already noted^{38} to develop in the numerical integration of the CFI. An exponential growth of the electrostatic field during the linear stage of the CFI, with a rate larger than *γ _{B}*, has been indeed observed and discussed in several works about the CFI.

^{30–35}A quantitative explanation of the rapid electrostatic mode growth—larger than

*γ*but generally different from $ 2 \gamma B$—was given in Ref. 30, in the framework of the linear theory of the CFI in a cold fluid model: it was therein suggested that

_{B}*γ*is a function of the asymmetry (i.e., the ratio $ n 1 / n 2$) of the two counter-propagating beams. Based on such analysis, however, this electrostatic growth disappears in a symmetric configuration (

_{B}*n*

_{1}=

*n*

_{2}), where an exact decoupling between electrostatic and Weibel-type modes is possible (as pointed out, e.g., also in Ref. 34): this effect is, therefore, different from the result discussed in Ref. 38, where the nonlinear electrostatic field growth, with rate $ 2 \gamma B$, has been highlighted in the cold fluid limit with

*n*

_{1}=

*n*

_{2}and has been attributed to the coupling between the Weibel mode and the Langmuir waves. More in general, the electrostatic growth with rate $ 2 \gamma B$ is recovered in the kinetic simulations of the CFI instability, which have been initialized with a small, yet finite temperature of two identical opposing electron beams (see, e.g., Fig. 1 of Ref. 31, Ref. 33, Fig. 23 and related comments in Ref. 34, and Fig. 2 and comments of Ref. 35). In these cases, again, the nonlinear nature of the electrostatic growth has been pointed out, together with the fact that it appears as a second-order effect, which becomes manifest at a “later time” with respect to the initial magnetic field growth. In Ref. 33, in particular, this effect was interpreted as due to the magnetic gradient force that accelerates electron currents, which, in turn, amplify the electrostatic field.

### B. Outline of the present article

In this work, we interpret the $ 2 \gamma B$ growth of the electrostatic field in purely e.m. Weibel mode as due to the second-order effect of charge separation, which is induced by the filamentation of the electron current density related to the spatial inhomogeneity of the magnetic field perturbation amplified by the Weibel mode. This electron density concentration, to which hereafter we refer for brevity as “magnetic lensing,” is due to the differential rotation in space of electron beams around the local maxima of the magnetic field: this tends to concentrate electrons halfway between the points of inversion of the magnetic field, by increasing the current density and by thus amplifying the magnetic field—mechanism that is at the basis of the CFI, indeed, to which Weibel instability can be related. This variation of the electron density, however, also induces a charge separation.

This mechanism is discussed via an analytical toy model in Sec. III, after having recalled in Sec. II some general features of both Weibel and CFI. Its implications are discussed in a kinetic non-relativistic framework, by considering in particular its effects on the anisotropic heating of the plasma. In Sec. IV, we discuss from an analytical point of view how the fluctuations of the electron pressure tensor components couple in an anisotropic way to the first-order magnetic fluctuations and to the second-order electrostatic fluctuations. This analysis is first performed in an extended electron fluid model where the pressure tensor equation is retained but the heat flux is neglected:^{39} this is a fitting approximation for the description of the linear evolution of the Weibel instability^{39–41} (Sec. IV A). We then revise (Sec. IV B) the key features of the reduced multi-stream kinetic model,^{42} specialized to describe the Weibel instability;^{37} we recall that this model is grounded on the notion of translational symmetry and, more specifically, on the invariance of the transverse canonical moment component associated with the presence of a cyclic variable. We then discuss (Sec. IV C) how the coupling between electron density and pressure fluctuations recognized in the extended fluid model can be identified in the reduced-kinetic multi-stream approach. In Sec. V, we support the previous analysis with three cases of nonlinear numerical simulations. To this end, we compare results obtained both with the “Eulerian” VLasov ElectroMagnetic solver VLEM,^{43} which integrates the Vlasov–Maxwell system with a semi-Lagrangian algorithm (Sec. V A), and with the multi-stream numerical model of Ref. 37 (Sec. V B). Conclusions follow (Sec. VI). Some features of the linear analysis of Weibel and CFI are recalled in Appendices A and B, respectively, and analytical details of the toy model for the “magnetic lensing” are provided in Appendix C.

## II. WEIBEL AND CURRENT FILAMENTATION INSTABILITY: GENERAL FEATURES

Let us first recall some essential features of both Weibel instability and CFI.

Ions are supposed to be unable to respond, over the timescales of interest, due to their larger inertia: they constitute a uniform static background with density *n*_{0}, with respect to which only electrons evolve.

*x*,

*y*)-plane, with $ P y , 0 > P x , 0$ and a perturbation with $ k = k e x$, the dispersion relation of Weibel instability

^{25}is obtained by linearizing this system around the neutral equilibrium distribution function for electrons,

*m*is the electron mass, and $ P i , 0 = n 0 k b T i$, for $ i = x , y , z$ are the non-zero components of an initially diagonal electron pressure tensor

*k*being Boltzmann constant. The pressure tensor component led us to introduce the (squared) thermal velocities $ c i 2 = P i , 0 / m$, where $ i = x , y , z$. From now on, labels 0 and 1 will be used for equilibrium and first-order perturbations, respectively.

_{B}Some mechanisms, which can lead to Weibel-unstable pressure configurations, have been discussed in Ref. 40. Key details of the linear kinetic theory of Weibel instability are recalled in Appendix A.

^{36}the (cold) CFI corresponds to the limit case of Weibel instability in which the temperature (pressure) component orthogonal to the wave vector goes to zero, and the “parallel” temperature (pressure) is given by the superposition of two counter-propagating cold electron beams. These can be labeled with $ \alpha = b , p$, which are here named “beam” (

*b*) and “plasma” (

*p*) following a common habit used for beam–plasma instabilities. In a statistical description, the only formal difference with respect to Weibel's case is in the initial distribution function,

*J*= 0, initially.

_{y}This configuration allows for a simple and intuitive interpretation of the magnetic field amplification in terms of Lenz law: the mechanism of CFI can be explained^{44} as due to the repulsion of the opposing current densities associated with the beams. For beams initially streaming along *y* and for perturbation along *x*, the magnetic field component $ B z , 1$ is this way amplified (cf. Fig. 1).

Some details of the linear theory of Fried's CFI are recalled in Appendix B, whereas in Appendix C, the essential features of the cold fluid description^{46} of CFI are summarized.

Some logical difference, however, exists between the kinetic modeling first proposed by Fried, based on a single Vlasov equation for the distribution of Eq. (3), and the typical macroscopical “picture,” which is usually given of CFI, when the two counter-propagating beams are identified as two *distinct* populations of electrons that separately evolve, although coupled by the Maxwell equations, according to their own dynamics and continuity equations. This difference leads to some discrepancies between the kinetic and fluid modeling (especially related to the evolution of higher fluid moments, when they are included^{39,47} so as to account for thermal effects in CFI), which have been preliminarily discussed in Ref. 41. As they are only of marginal interest for this work, we will address them in greater detail somewhere else. In the following, we do not need to discuss these subjects any further, since we will focus on the “original version” of the temperature anisotropy-driven instability by Weibel. References to CFI will be made, occasionally, only for illustrative purpose.

### A. Energy source of Weibel-type modes

The cold-CFI picture mimics the single particle electron dynamics. This allows one to recognize how the current filamentation mechanism of CFI, which is driven by the *momentum* anisotropy, is at the basis also of the magnetic field amplification in the Weibel instability, which is formally driven by the *pressure* anisotropy, instead. In a broader sense, both instabilities are thus “momentum driven,” although the source of energy in CFI is strictly “(ordered) kinetic,” whereas in Weibel's case, it is “thermal-kinetic,” i.e., it is the internal energy. This can be formally seen by looking at the total energy conservation of the system.

^{48}in the $ ( x , v )$ configuration space as $ f ( x , v , t )$, Vlasov equation reads

*total*energy density of the electrons, $ E = ( m / 2 ) \u222b v 2 f ( x , v , t ) \u2009 d 3 v$, fulfills

*K*= 0 in both Fried's kinetic version of the CFI and in Weibel instability, since $ u = 0$. Note that in the kinetic description of CFI, $ u = 0$, since the distribution function $ f ( x , v x , v y , 0 ) = f 0$ of Eq. (3) must fulfill $ n = n b + n p$ and $ n b v b + n p v p = 0$, in order for a null current density to be granted at the equilibrium (this is one of the differences with respect to a reduced-fluid or to a multi-beam approach, in which each beam is assumed to evolve according to a distinct Vlasov equation, instead). One sees from the dispersion relations

^{25,36}that the source of energy in both the kinetic CFI and the Weibel instability is in the difference between $ E T , x$ and $ E T , y$ (or $ E T , z$) of Eq. (12).

## III. CURRENT FILAMENTATION MECHANISM AND ELECTRON DENSITY CONCENTRATION

A quantitative characterization of the current filamentation mechanism shown in Fig. 1 is provided by the multi-stream model,^{42} which we will discuss later, by means of which a formal correspondence^{37,49} can be established between the linear Weibel instability and a configuration consisting of two co-propagating cold electron beams and a third cold electron population with zero fluid velocity.

In order to make useful estimates and to get a deeper insight of the filamentation mechanism, it is, however, possible to take a slightly different approach, simpler than the full kinetic treatment, in which the role of the single particle dynamics is made more explicit. Of course, this simplification must be paid in other means, namely, by making use of some heuristic estimate, which is based on the knowledge of the results of kinetic linear analysis (and which, as we will discuss at the end of this section, prevents us to recover the exact kinetic results).

To this purpose, let us first postulate a magnetic field perturbation in agreement with the Weibel instability case, so as to write $ B z ( x , t ) = B M \epsilon ( t ) \u2009 cos ( k x x )$. Here *B _{M}* is a reference value of the magnetic field amplitude—

*a posteriori*we will make it correspond to the maximum amplitude that the magnetic field attains at saturation—and $ cos ( k x x )$ is its spatial modulation corresponding to the wave-vector

*k*of the perturbation. In this writing, the spatial and temporal dependence of $ B z ( x , t )$ have been split: the function $ \epsilon ( t ) = \epsilon 0 e \gamma B t$ expresses the time dependence of the magnetic field amplitude that exponentially increases with growth rate

_{x}*γ*. In the comparison with Weibel modes, $ \epsilon ( t )$ can be read as a dimensionless parameter weighting the relative amplitude of the magnetic field perturbation at each time

_{B}*t*, until saturation is reached when $ \epsilon ( t ) \u223c O ( 1 )$.

In agreement with the results of the linear theory of Weibel-type modes, we can then assume that the time and space dependence can be separated for all scalar functions *F*, in the form $ F ( x , t ) = X F ( x ) T F ( t )$.

Maxwell's equations impose specific phase relations on the spatial and temporal components of the electric field and current density *E _{x}* and

*J*, which are detailed in Appendix C [Eqs. (C2) and (C3)]. It is then possible to establish a relation between Eqs. (C2) and (C3) and the function assumed for $ B z ( x , t )$ without resorting to the solution of the full Vlasov–Maxwell system but by looking instead at the previous set of equations from the point of view of the single particle dynamics. To this purpose, let us consider

_{x}*N*electrons, which at

*t*= 0 would be uniformly distributed over a volume $ V 0 = L x \xd7 L y \xd7 L z$ and would be at rest, if it were not for the motion related to the electron current density $ \delta J y$, related to an initial magnetic field perturbation $ \delta B ( x , t ) | t = 0 = ( 0 , \u2009 0 , \u2009 B z ( x , 0 ) )$ through Ampère's law. The number of electrons per unit volume

*V*

_{0}can change over time and depends on the region of space in the neighborhood of the coordinate value

*x*. Therefore, we write

*N*(

*x*,

*t*) and $ n ( x , t ) = N ( x , t ) / V 0$, with $ n 0 = n ( 0 , 0 )$. The symbol “ $ \u2026 \u0303$” we will use for the single electron velocity $ v \u0303 ( x , t )$ marks the difference with respect to the phase-space velocity coordinates $ v$, independent from

*x*, used before. Also note that, in this way, we are not self-consistently modeling the Weibel instability (the initial electron distribution lacks the momentum anisotropy that drives the instability), but we are only modeling the self-consistent dynamics of a population of electrons that is set into motion by an external, Weibel-like e.m.

*forcing*and back-reacts to it. Analytical expressions for $ v \u0303 x = \delta v \u0303 x , \u2009 v \u0303 y = \delta v \u0303 y , \u2009 n = n 0 + \delta n = ( N ( t = 0 ) + \delta N ) / V 0$, and $ J = \delta J , \u2009 B = \delta B , \u2009 E = \delta E$, where the quantities with

*δ*in front would be zero if $ B M = 0$, can be thus obtained by combining Maxwell's equations with Newton's equations and the continuity equation for a bunch of electrons described according to this toy-model.

These orderings—and the way they have been obtained (cf. Appendix C)—represent the main result, here. This helps us to shed light on some features, already known but somewhat “debated” in previous literature, of the nonlinear coupling between e.m. and electrostatic perturbations in Weibel modes. In particular, this analysis shows that the spatial modulation of the magnetic field induces the particle dynamics a concentration (or depletion) of electrons, in the regions around the local maxima (or minima) of the magnetic energy density $ | B ( x ) | 2$. For the chosen magnetic field perturbation, these regions correspond to the extrema of $ cos ( 2 k x x )$. This is a second-order effect with respect to the parameter $ \epsilon \u226a 1$ and to the spatial modulation $ cos ( k x x )$ of the magnetic field perturbation *B _{z}*, amplified by the Weibel instability. For future reference, we can call this effect “magnetic lensing.” It is a

*quadratic*—and thus “nonlinear”—effect of the current filamentation linearly correlated with the magnetic field perturbation [cf. Eq. (C1)], and it is one of the ingredients of the

*magnetic trapping*, which in their advanced nonlinear stage induces saturation of both the CFI and the Weibel instability (cf. Sec. IV C 1).

Finally, note that the analytical solutions found here for the perturbed quantities are accurate only for what concerns the spatial and time dependence, but not for what concerns multiplicative factors that depend on the plasma parameters. This is due to the non-self-consistent character of this modeling, in which the exponential growth of e.m. mode is only postulated. In other words, we have obtained analytical solutions that are compatible with an exponential growth of the magnetic field component *B _{z}* and provide a crude modeling of the mechanism depicted in Fig. 1, but they fail to account for the reservoir of energy driving the instability, not only for its saturation mechanism and but also for the quantitatively correct dispersion relation. This can be seen by comparing the growth rate of Weibel instability, which in the limit

^{25}$ T y \u226b T x$ reads $ \gamma W \u2243 k x c y / 1 + k x 2 d e 2$, with the result obtained by comparison of Eqs. (C1) and (C9) (cf. Appendix C), once the condition $ J y \u223c \u2212 n 0 e \delta v \u0303 y$ implied by $ \delta n \u223c \epsilon 2$ is used: the result $ \gamma B = \omega p e 1 + k x 2 d e 2$ so obtained not only lacks a cutoff wavelength—feature that is typical

^{39,46,50}of the cold fluid modeling of CFI—but it also displays the wrong spatial dispersion. In the expressions above, we have introduced the electron skin depth $ d e 2 \u2261 c 2 / \omega p e 2$ evaluated with respect to

*n*

_{0}.

## IV. INSIGHTS ON THE COUPLING OF THERMAL AND e.m. PERTURBATIONS FROM REDUCED MODELS

The “magnetic lensing” mechanism discussed above explains in dynamical terms how electrostatic fluctuations can nonlinearly couple to otherwise purely e.m. fluctuations in the evolution of Weibel instability. Since at these fast timescales electron density fluctuations relate both to charge separation effects and to pressure fluctuations, the spatial filamentation related to this electron density concentration can generally induce a coupling between electrostatic and thermal effects in non-uniformly magnetized plasmas.

The way this occurs is non-trivial and requires a kinetic treatment, which in turn leads to Vlasov or particle-in-cell (PIC) numerical approaches. These are powerful tools, but the richness of physics they describe has sometimes the drawback to make it difficult to single out the role played by the different processes at play. Reduced kinetic models, such as the multi-stream model,^{37,49,51,52} or extended fluid models^{39–41,47} that include the dynamics of the second-order anisotropic moment (i.e., the pressure tensor) of Vlasov equation, help in this: in addition to being self-consistent within well identifiable limits, they are particularly useful for the interpretation of the kinetic numerical results that we are going to discuss below, since they provide immediate insight in the dynamics of the lower order fluid moments—namely, density, fluid velocity, and kinetic pressure components—that are typically associated with macroscopically measurable quantities. The multi-stream model, moreover, provides a natural link between the kinetic and fluid description, as it has been preliminarily discussed in Ref. 41.

Information on the statistical distribution of particles is differently retained in the dynamics equations of each of these “reduced” models, but in spite of this, the energy equation always remains valid in the form of (11), provided its last right-hand side term is rewritten using the relevant heat flux and pressure tensor components. Note that the negligible effect, which the heat flux tensor has been shown to have on the fluid description of the linear dynamics of these e.m. instabilities,^{39} agrees with the fact that the right-hand side term of Eq. (11) does not contribute to the volume energy integral.

In the remainder of this section, for each type of reduced approach, we first recall some of the essential features of the model, which have been already detailed elsewhere. Then, we discuss some non-trivial consequences that can be drawn from each model about the description of Weibel-type modes and about the way that macroscopic quantities, such as density and the components of the fluid velocity and the pressure tensor, can mediate the nonlinear coupling between fluid and e.m. fluctuations.

### A. Fluid description of the Weibel instability

A linear description of the pressure-driven Weibel instability is made possible in an extended fluid description^{39,47} in which the full pressure tensor dynamics is retained. Reference 47 first considered the hydrodynamic limit, in which $ \xi 2 \u2261 k 2 c x 2 / \omega 2$ corrections are neglected [cf. Eqs. (A1)–(A3) and Appendix A]. In Ref. 39, it was, however, shown that a quite accurate description of the kinetic dispersion relation is possible also for finite values of $ \xi 2$, if the dynamics of the full pressure tensor $\Pi $ is retained while the heat flux gradient is completely neglected. This approximation was explicitly discussed in Refs. 39 and 41, where it was pointed out that the heat flux gradient brings relevant corrections only near saturation of the instability or, possibly, close to the cutoff wave-vector conditions, especially of the warm CFI. The efficiency of this fluid approach in describing the linear Weibel has been tested in the aforementioned works and also led to non-trivial predictions about the existence of roots of the kinetic dispersion relation of oblique Weibel-type modes.^{53}

^{54}This restriction is consistent with the use of the non-relativistic Vlasov equation for electrons when ions are at rest, which is adopted in the classical kinetic Weibel's approach (cf. Appendix A). This means that we must retain in the model the displacement current in Ampère's law [cf. Eq. (9)].

^{39}$ [ D ]$ with

Retaining instead the first-order expansion of (18) with respect to $ k 2 c x 2 / \omega 2$ yields the dispersion relation (A4), which coincides with the limit discussed by Basu in Ref. 47.

*D*= 0. Linearization of the fluid system (see Ref. 45, Sec. 3.3.4.1) yields the following equations, which are related to the velocity components $ u x , 1$ and $ u y , 1$:

_{yy}*x*,

*y*coordinates contribute to the electron fluid dynamics in the (

*x*,

*y*) plane, via the coupled equations,

*n*and $ u x , 1$, in agreement with the magnetic lensing mechanism discussed in Sec. III. The coupling between the in-plane diagonal components of $\Pi $ and the electrostatic fluctuations related to the gradients of $ u x , 1$ can be identified as due to electron density fluctuations, i.e., to charge density fluctuations, as it can be seen by noting that in this one-dimensional coordinate dependence, Eqs. (22) and (23) can be read as the linear form of the polytropic-type equations,

*ε*the “instantaneous” characteristic amplitude of the Weibel mode during its linear stage, we can consistently order the contributions of the fluid momentum equation as

_{yy}that can be deduced from its evolution equation,

The ordering $ \Pi x y \u223c \epsilon $ in Eqs. (28) and (29) has been chosen for consistency with Eq. (28), in which the last right-hand side term, proportional to *B _{z}*, is of the same order. As we will see, this non-trivial result is confirmed by full kinetic simulations (cf. Sec. V A 1). From Eqs. (27)–(29) follow indeed the non-isotropic way by which pressure fluctuations (i.e., the components of the pressure tensor) couple with electrostatic and magnetic fluctuations for the case of a growing, purely e.m. Weibel mode. Only the in-plane diagonal components of $\Pi $ are ordered like density and electrostatic fluctuations (i.e., they are nonlinear fluctuations), whereas the off-diagonal elements are ordered like the magnetic perturbation (i.e., they are linear fluctuations). This feature becomes manifest in the fluid description only when the degrees of freedom related to the evolution of the full pressure tensor components are retained. The fact that, depending on the direction of the spatial inhomogeneity or of the wave-vector, the components of the pressure tensor are anisotropically affected by the transfer of energy from the bulk fluid motion (e.g., in a shear-driven anisotropization mechanism

^{15}) or by linear mode oscillations (e.g., in the propagation of magneto-acoustic waves

^{55}), which have been already evidenced in extended fluid models describing Vlasov dynamics at ions scales. In this regard, it should be noted the consistency of the ordering $ \delta \Pi x x \u223c \delta \Pi y y \u223c ( \delta \Pi x y ) 2$ for the Weibel instability, with the fact that only the trace of the pressure tensor contributes to the internal energy $ E T , x + E T , y$ of the plasma: in diagonalizing the Π

_{ij}matrix for $ i , j = x , y$ so as to determine the local contributions to the internal energy, the off-diagonal corrections $ \Pi x y = \delta \Pi x y \u223c \epsilon $ contribute like the diagonal variations of order $ \epsilon 2$ since $ \Pi x x = P x , 0 + \delta \Pi x x \u223c 1 + \epsilon 2$ and $ \Pi y y = P y , 0 + \delta \Pi y y \u223c 1 + \epsilon 2$. This is analogous to the constraint on the weight of finite gyro-radius (FLR) corrections in the diagonal and off-diagonal terms of the dispersion matrix, which play a key role, for example, in the accounting of the correct dispersive kinetic effects of magneto-acoustic modes (see Refs. 56 and 57 and the more recent Ref. 55).

*K*. Therefore, $ K = K x + K y + K z$. We have also neglected the energy contributions that are not relevant to the geometry we have considered (i.e., $ \delta K z = \delta E T , z = 0$).

#### 1. Weibel instability and “generalized” Ohm's law

^{54}the same way non-polytropic closures are introduced in the so-called extended magnetohydrodynamic models. The evolution of

**in EMHD is as usual obtained by combining Faraday equation with the curl of (16)**

*B***, and from this Alfvén theorem, the frozen-in law and the conservation of magnetic helicity all follow (see, e.g., Ref. 58).**

*u*The violation of the “magnetic field topology” is allowed when the Lagrangian invariance of ** B** alone is broken. However, differently from magnetic reconnection, the Weibel-type instabilities of interest here do not require, strictly speaking, such violation to occur. Nevertheless, since in the geometry adopted in this work we require a (static) magnetic field component $ \delta B z$ to grow, we can reason similarly to the reconnection case, so as to look at the terms that allow $ \u2202 B / \u2202 t \u2260 0$ (cf. e.g., Refs. 59 and 60). Even in the absence of the pressure-related contribution, this can happen because of a non-zero rotational of the fluid velocity: this is the case of the CFI. Or, it can happen even if the rotational of the fluid velocity is null, provided the last term of Eq. (32) is non-zero: this is the case of both the Weibel instability and the Biermann's battery mechanism

^{61}(cf. Ref. 62). However, while Weibel instability requires a seed of the magnetic field perturbation, $ \delta B$, which is exponentially amplified (i.e., it is a linear instability), Biermann's mechanism can generate a magnetic field from zero that grows at a much slower rate (typically, almost linearly in time) than in Weibel's case, provided the equilibrium density and pressure tensor satisfy $ \u2207 n 0 \xd7 \u2207 \xb7 \Pi 0 \u2260 0$.

#### 2. Fluid vs single-particle dynamics

*d*Φ

_{Ω}(

**)/**

*B**dt*

_{ $ v \u0303$}≡ ∂Φ

_{Ω}(

**)/∂**

*B**t*+ $ v \u0303 \xb7 \u2207 \Phi \Omega $(

**), and which Faraday's equation implies to satisfy**

*B**cold*electrons), the magnetic flux must be conserved and the magnetic field amplification follows from the shrinking of the integration surface around the local extrema of the magnetic field, in turn related to the “magnetic lensing” mechanism previously discussed, under a global, fluid perspective the magnetic field grows

*because of the violation*of the magnetic flux conservation: a non-zero contribution of the right-hand side terms of (38) is needed if the surface Ω is chosen to cover the entire subdomain where the sign of

*B*is fixed (positive or negative).

_{z}This means that in Weibel instability, the growth of the magnetic field could be interpreted as due to a statistical dominance of electrons whose trajectories enclose surfaces, which, on average, shrink over time. This is compatible with the instability condition $ P y , 0 > P x , 0$, since in the rest frame of the electron plasma, the initial velocity of each particle is given by its thermal motion. The trajectory of a particle initially moving in the *x* direction, along which a magnetic field perturbation is sinusoidally modulated, will experience an alternation of deflections toward and away from the magnetic extrema, depending on the local sign of $ \delta B z ( x )$. For a spatially periodic modulation of $ \delta B z$, this alternation will give a null net contribution to the flux variation over a surface element enclosing the spatial fluctuations of the magnetic field. Instead, particles moving along *y* will undergo the change of trajectory sketched in Fig. 1 that leads to a global magnetic field increase. An oblique motion will encompass both effects, but the global growth of the magnetic field is granted by the statistical dominance of particles having $ v \u0303 y > v \u0303 x$, provided by the instability condition $ P y , 0 > P x , 0$. This allows a re-interpretation of the role of the pressure gradient term in Eq. (36) from the perspective of the single particle dynamics.

#### 3. Electrostatic growth due to longitudinal current density perturbations

The fluid model discussed in this section allows us to formally recover the interpretation provided in Refs. 27 and 29, where the growth of the electrostatic field was explained as needed to compensate for the current density generated along the wave-vector of the perturbation. We can show that thermal fluctuations play a key role in this.

*δt*a solenoidal electric component $ \delta E y ( x , \delta t ) \u223c [ \epsilon ( 0 ) ] \delta t$ follows from Ampère's equation (9) and $ \delta \Pi x y ( x , \delta t ) \u223c [ \epsilon ( 0 ) ] \delta t$ follows from Eq. (24). At the next step, these perturbations generate $ \delta u y ( x , 2 \delta t ) \u223c [ \epsilon ( 0 ) ] \delta t 2$ [cf. Eq. (21)] and $ \delta \Pi x x ( x , 2 \delta t ) \u221d [ \delta \Pi x y ( x , \delta t ) \delta B z ( x , \delta t ) ] \delta t \u2243 [ \delta \Pi x y ( x , \delta t ) \delta B z ( x , 0 ) ] \delta t \u223c [ \epsilon ( 0 ) ] 2 \delta t 2$ [cf. the

*xx*-component of Eq. (17)]. At the next step, due to Eq. (20) and using $ \delta B z ( x , 2 \delta t ) \u223c \u2202 x [ \delta E y ( x , \delta t ) ] \delta t \u223c [ \epsilon ( 0 ) ] \delta t 2$ [Eq. (8)], we obtain

*x*-component of Ampère's law requires then that an electrostatic component $ \delta E x | 3 \delta t = \u2212 4 \pi \delta J x | 3 \delta t$ exists, as already noted in Refs. 27 and 29. This formal reconstruction of the steps sketched in the cartoon of Fig. 1 elucidates the role that pressure fluctuations play in this collective process, in agreement with the discussion provided at the end of Sec. IV A 2 when a perturbation is given just on the magnetic field. Should the perturbation be given in a more self-consistent way in fulfillment of Ampère's equation, the essence of the argument would not change: only the ordering with respect to the power of

*δt*would change for some quantity, and the first perturbation on $ \delta u x$ would occur already at the first step as $ \delta u x ( x , \delta t ) \u221d [ \delta u y \delta B z ] | 0 \delta t \u223c [ \epsilon ( 0 ) ] 2 \delta t$.

### B. Multistream models for the Weibel instability

^{42}so as to provide the exact evolution of a sub-class of solutions of the Vlasov–Maxwell problem. It is, therefore, suited to describe purely transverse Weibel-type instabilities,

^{37,49}in which spatial dependence is assumed only along one direction (

*x*, in this case). Let us name $ P c \u22a5 = p \u22a5 \u2212 ( e / c ) A \u22a5 ( x , t )$ the invariant components of the electron canonical momentum

*perpendicular to the wave-vector*of the perturbation (in the Coulomb gauge $ \u2207 \xb7 A \u22a5 = 0$). The multi-stream model follows the evolution of each beam of particles having the same value of the $ P c \u22a5$ component, by associating each of them with its own Hamiltonian. Naming $ C N$ the numerical, constant value of the canonical momentum component of an electron of the “

*N*th” bunch, its Hamiltonian in the non-relativistic limit reads

In the 2D geometry configuration, we consider hereafter, the dynamics along the *z-axis* is neglected. We can thus take $ C N = C N e y$ and $ A \u22a5 = A y e y$.

Also note that, for purely transverse Weibel-type instabilities, which are fully decoupled from electrostatic two-stream-like instabilities, the distribution function $ f ( x , v x , v y , t )$ is spatially homogenous at equilibrium. Therefore, the dependence on *x* of each *F _{N}* of Eq. (46) is only due to the perturbed density.

Finally, note that assuming that each beam evolves with its own Vlasov equation (41) makes the multi-stream model alike to the fluid description of the CFI, in which each beam is let separately evolve according to its own momentum equation. Without such assumption, the set of equation (41) would be replaced by $ \u2202 ( \u2211 N f N ) / \u2202 t + [ \u2211 N H N , \u2211 N f N ] = 0$. This difference may lead to some discrepancies with respect to the full kinetic treatment, for example, if the number of beams is not large enough to grant a sufficiently dense sampling of the momentum coordinate in their phase-space so that the initially well-separated momenta evolve by “overlapping” with respect to that error. This, however, is a drawback shared by all reduced kinetic models. A well-known example is provided by the gyrokinetic model and reduced models derived from it, where alike particles are made separately evolve according to the value of their energy and/or of their magnetic moment. In the remainder of this work, we do not address further this point: we assume the validity of (41) and, in the following, we will interpret in this light the possible differences between results obtained from full Vlasov–Maxwell and reduced multi-stream simulations. We will see that having separated the evolution equation of each beam demands special care in the evaluation of quantities that are averaged over the whole electron population (i.e., over the sum of the beams). This especially concerns the calculation of higher order fluid velocity moments like the pressure tensor components.

#### 1. Case of “cold” multi-streams

*v*velocity component. This means we take

_{x}*N*fluid equations obtained from the zeroth and first-order moment of (42),

*B*. In the case of the Weibel instability, one typically neglects $ e A y ( x , 0 ) / c$ so that $ u y , N ( x , 0 ) = u y , N ( 0 ) = C y , N$.

_{z}*N*.

^{37}that at least three cold streams, one of which must have zero velocity, are necessary and sufficient to recover the linear dispersion relation of the transverse Weibel mode: the two streams with non-zero velocity are necessary to fix an initial $ P y , 0$ pressure component, whereas the stream with zero velocity is needed to provide a “bulk” electron density corresponding to the center of a Maxwellian distribution. Further correspondences between the two models, also in the nonlinear evolution of the instability, have been discussed in Refs. 41, 49, 51, 52, and 63.

#### 2. Case of “warm” multi-streams

Integrating the reduced Vlasov equation for a distribution (46) in which $ f N ( x , v x , t )$ differs from a delta of Dirac allows a more complete inclusion of thermal effects. To this end, an improvement from the point of view of both analytical simplicity and computational cost consists in further reducing the model for the evolution of the $ f N ( x , v x , t )$ distribution. The “waterbag” approach can be used for this. It consists in expressing $ f N ( x , v x , t )$ as a superposition of waterbag distributions.^{64,65} For Weibel-type instabilities, this approach has been first used for relativistic distributions^{66} and then for non-relativistic ones, too.^{49,63} Operationally speaking, we express $ f N ( x , v x , t )$ in terms of Heaviside step functions having value $ F l , N > 0$ for $ v x \u2208 [ \u2212 w N l , \u2212 ( x , t ) , + w N l , + ( x , t ) ]$ and zero otherwise, $ w N l , \xb1 ( x , 0 )$ being assigned at the initial time *t* = 0 for each waterbag *l*.

^{−1}. Owing to Liouville's theorem, the dynamics of the bag is completely determined by its contours,

^{64,65}that is, by

*w*corresponds to the fluid velocity component $ u x , N$ of Eqs. (48) and (49) of the cold stream case. Note that the closure condition (64) also defines the $ \Pi x x , N$ component in this modeling. Equation (64) can be directly computed from the first-order moment of Vlasov equation, by identifying $ m n N 3 / ( 12 n 0 2 A N 2 )$ as $ \Pi x x , N$ after direct comparison with the fluid equations. The argument can be generalized to the case of many bags, say for $ l = 1 , . . , M$, where Eq. (60) is substituted by

_{N}*N*th stream reads

_{xx}in the reference frame of the whole multi-waterbag distribution [i.e., using (70)] and from the evaluation of the second-order moment of Vlasov equation for the distribution (65), one deduces that

### C. Correspondence between full kinetic, extended-fluid, and reduced-kinetic descriptions

*f*(namely, initial density, average velocity, and “temperature”), this allows for an “exact” correspondence between the two linear descriptions of the instability.

*n*and $ C y , N 2$. Then, compatibly with the non-relativistic limit $ c 2 \u226b ( \omega / k ) 2$, Eq. (75) reduces to Weibel's result,

_{N}^{25}

#### 1. Magnetic lensing and magnetic trapping in the multi-stream approach

*N*characterized by the same value of the canonical momentum $ P y , N$, uniformly “sees” the effects of the electromagnetic field from a local point of view. In particular, given a spatial modulation of the magnetic field, particles with equal values of $ C y , N$ will be equally affected by the “magnetic lensing” mechanism described in Sec. III and will also equally experience a

*magnetic trapping*mechanism with a bounce frequency $ \omega b N ,$ expressible in terms of the peak value of the magnetic field (see details of the calculations in Ref. 52, and Sec. III for some definitions of quantities),

*E*are much weaker in amplitude, when compared to those on

_{x}*B*, their presence is crucial for mediating resonant wave–particle interactions in the Weibel instability. Electrostatic perturbations can become visible in the phase-space patterns of the distribution function, although it is

_{z}*a priori*difficult to recognize them, since they are generally masked by the more intense effects of the magnetic field, such as the magnetic trapping. However, the multi-stream approach allows us to neatly identify their role in the electron population associated with the “central” stream (

*N*= 0) with null canonical momentum, which corresponds to the value $ C y , N = 0 = 0$ and, therefore, to a zero magnetic bounce frequency (cf. Fig. 15 in Sec. V B 1).

## V. NONLINEAR ELECTROSTATIC COUPLING IN VLASOV–MAXWELL SIMULATIONS

We now discuss the results of nonlinear numerical simulations of the Weibel instability, in light of the previous analytical results. We first analyze two nonlinear simulations performed with the 1D-2V version of the semi-Lagrangian relativistic Vlasov–Maxwell solver VLEM.^{43} Then we will compare one of them with a simulation performed with the multi-stream numerical model of Ref. 63. Both the semi-Lagrangian code and the multi-stream code, thanks to their Eulerian phase-space structure, are initially noiseless, since they are based on a statistical approach where the graininess parameter $ g = 1 / ( n 0 \lambda D ) d$ is formally set to zero, *λ _{D}* being the electron Debye length and

*d*a parameter characterizing the dimensionality of the system. With this approach, it is then possible to follow first the emergence and then the growth of the various components of the electromagnetic field from the very beginning of the simulation, and with a remarkable accuracy. In particular, it is necessary to initialize the marginally stable equilibrium with a perturbation, in order to start the instability, due to the lack of numerical noise. This makes any other perturbation start with epsilon-machine amplitude.

### A. Full Vlasov–Maxwell simulations

Both simulations performed with the 1D-2V version of the solver VLEM^{43} have been initialized with the distribution function (1), perturbed with a magnetic fluctuation modulated along the coordinate *x*, only, so as to excite only the purely transverse Weibel instability.

(

*S1*) In the first simulation, the initial pressure components are $ P y , 0 = 20 P x , 0 = 20 \u2009 keV$. Since for this value of the pressure anisotropy the linear dispersion relation $ \epsilon y y = 0$ given by Eq. (A2) predicts a maximal growth rate $ \gamma M = 0.129 \omega p e$ for $ k M d e = 1.462$, we choose a spatial domain with extension $ L x = 2 \pi / k M$, so as to make certain that we destabilize such fastest growing mode. The magnetic field perturbation is given on the*z*-component as $ B z ( x , 0 ) = B 0 \u2009 cos ( k M x )$, with $ B 0 = 10 \u2212 5$. This is not an exact eigenmode of the system, which implies that after a transient time also some electrostatic fluctuations can be destabilized. However, the initial amplitude of these perturbations is much smaller than the starting $ \u223c 10 \u2212 5$ amplitude of the growing Weibel mode. The other numerical parameters are*N*= 512, $ N v x = N v y = 96$, and $ \Delta t = 0.005 \u2009 \omega p \u2212 1$. In this simulation, we focus on the correlation of density fluctuations and electrostatic fluctuations on the one side and of magnetic fluctuations on the other one, with the linear perturbations of the pressure tensor components (cf. Sec. IV A)._{x}(

*S2*) The second simulation, in which we focus instead on the phase-space dynamics, is initialized with similar physical parameters to (*S1*), $ P y , 0 = 25 P x , 0 = 25 \u2009 keV$, which correspond to a fastest Weibel mode with $ \gamma M = 0.142 \omega p e$ for $ k M d e = 1.40$. Again, we have fixed $ L x = 2 \pi / k M$, and the system is perturbed with $ B z ( x , 0 ) = B 0 \u2009 sin ( k M x )$, with $ B 0 = 10 \u2212 5$. The other numerical parameters are*N*= 256, $ N v x = N v y = 256$, and $ \Delta t = 0.005 \u2009 \omega p \u2212 1$. The difference with respect to the simulation case (_{x}*S1*) is in a slightly larger perpendicular temperature than before (so as to “enhance” the evolution of the Weibel mode, of which we are in this case interested also in the nonlinear stage) and in the larger resolution in the velocity space (somehow “compensated” by a smaller resolution in the space coordinate), here chosen so as to better follow the dynamics of the distribution function in the phase-space.

In both simulations, whenever not specified, times must be read as normalized to $ \omega p \u2212 1$, lengths to *d _{e}*, energies to

*mc*

^{2}, the components of the electric field to $ m \omega p c / e$, and the components of the magnetic field to $ m \omega p / e$

#### 1. Simulation case (*S1*): Electrostatic coupling via magnetic lensing and pressure tensor fluctuations

Figures 2–6 correspond to the simulation case (*S1*).

Figure 2 shows the time evolution of the total magnetic energy $ \u03f5 m , z = 1 / ( 8 \pi L x ) \u222b 0 L x B z 2 \u2009 d x$, the (inductive) electric energy component $ \u03f5 e , y = 1 / ( 8 \pi L x ) \u222b 0 L x E y 2 \u2009 d x$, and the electrostatic energy $ \u03f5 e , x = 1 / ( 8 \pi L x ) \u222b 0 L x E x 2 \u2009 d x$. Both the energy components related to $ B z 2$ and $ E y 2$ grow with a rate $ \u223c 0.26 \u2243 2 \gamma M$, in excellent agreement with the linear theory, whereas $ \u03f5 e , x$ grows with a rate, which is about $ \u2243 4 \gamma M$, in agreement with the fact it is generated by nonlinear coupling. While the components related to the Weibel mode start with comparable amplitudes, the perturbations related to the electrostatic fluctuations start with an amplitude, which is about 15 orders of magnitude smaller (i.e., $ B z 2 \u223c 10 30 E x 2$).

The transient time related to the fact we have initialized the system with a perturbation, which is not an exact eigenmode, is about $ \u223c 10 \u2009 t \omega p e$, if we look at the energy components related to the Weibel mode, but it is even shorter ( $ \u223c 5 \u2009 t \omega p e$) if we measure it with respect to the electrostatic component $ \u03f5 e , x = \u222b E x 2 \u2009 d x$. On the one hand, this highlights the role that the nonlinear coupling between electrostatic and e.m. fluctuations has in the building up of the correct eigenmodes, since we start from an approximated eigenmode expressed as a fluctuation on $ B z , 1$ alone. On the other hand, it should be noted that such an early growth of the electrostatic component due to the nonlinear coupling has been never observed in previous numerical kinetic simulations of both Weibel instability and CFI, where it has been usually noted,^{27–29,31,33–35} instead, only in the nonlinear stage of the instability: here it is made possible by the initial almost noiseless character of the Eulerian-type semi-Lagrangian Vlasov solver used. The small-amplitude temporal oscillations visible in the magnetic energy component during the transient interval have a frequency close to twice that of the e.m. wave propagating in the plasma, $ \omega e m = \omega p 2 + ( k M c ) 2 \u2243 \omega p$. All energy components indicate that saturation occurs around $ t \omega p \u2243 75$.

Figure 3 confirms the nonlinear origin of the electrostatic field growth. It shows the spatial profiles of the electron density *n*(*x*, *t*) and of the magnetic field $ B z ( x , t )$ at time $ t \omega p e = 0.5$, i.e., at an early time after initialization, and compares them with the corresponding profiles at *t* = 5: while the magnetic field is spatially modulated over the fundamental mode $ k M c / \omega p e = 2.017$, the electron density oscillates with $ 2 k M$. Owing to the “current filamentation/magnetic lensing” mechanism discussed in Sec. III and related to Fig. 1, we deduce from the results displayed in Fig. 3 that electrons are partially concentrated around regions where *B _{z}* = 0 and where the derivative of

*B*with respect to

_{z}*x*is positive or negative, depending on whether the sign of the

*v*electron velocity component is positive or negative, respectively. This concentration is initially weak but sufficient to create a local charge imbalance and, therefore, an electrostatic field that grows in time as the Weibel mode grows. This growth is at twice the $ \gamma B \u2243 \gamma M$ growth rate, as discussed in Sec. III. It must be noted that the variation of the $ u x ,$ fluid velocity component is itself a second-order fluctuation related to the nonlinear electrostatic coupling, that is, $ u ( x , t ) \u2212 u x , 0 \u2243 u x , 2$, with respect to linearized quantities (not shown here).

_{y}Figure 4 displays the time evolution of the fluctuations of the diagonal pressure tensor components, i.e., $ \delta \Pi x x ( x , t ) = \Pi x x ( x , t ) \u2212 P x , 0$ and $ \delta \Pi y y ( x , t ) = \Pi y y ( x , t ) \u2212 P y , 0$, evaluated in the position where their maxima are respectively localized (cf. Fig. 5). The components $ \Pi i j ( x , y )$ are computed from the kinetic simulations by integrating the second-order moment of the distribution function according to $ \Pi i j ( x , y ) = ( 1 / \u222b f d 3 v ) \u2009 \u222b ( v i \u2212 u i ) ( v j \u2212 u j ) f d 3 v$ and $ u i = ( 1 / \u222b f d 3 v ) \u2009 \u222b v i f d 3 v$. The logarithmic scale plots in Fig. 4 show that both the $ \Pi x x , 1$ and $ \Pi y y , 1$ grow exponentially with a growth rate $ 2 \gamma M$, which equals that of the electrostatic fluctuations. Coherently, the spatial profiles of $ \delta \Pi x x ( x , t )$ and $ \delta \Pi y y ( x , t )$, shown in Fig. 5 at $ t \omega p e = 0.5$, display the $ \u223c cos ( 2 k m x )$ spatial oscillation of the first-order harmonic of the fastest growing mode like the density fluctuations of Fig. 3.

The off-diagonal component $ \delta \Pi x y = \Pi x y ( x , t )$ (we recall that at equilibrium $ \Pi x y = 0$), instead, exponentially grows in time and is modulated in space as the linear Weibel mode, i.e., it has growth rate *γ _{M}*, and spatial oscillation corresponding to the wave vector

*k*(Fig. 6).

_{M}In summary, these numerical results confirm the orderings of Eq. (30).

#### 2. Simulation case (*S2*): Electrostatic coupling via magnetic lensing and phase-space dynamics

Figures 7–13 correspond to the simulation case (*S2*).

The saturation mechanism is linked to magnetic trapping, which appears in the form of a trapping structure on the $ k M = k 0$ mode (which is the fundamental mode *k*_{0}, here). However, also the $ 2 k M$ harmonic plays a key role in it. For a more detailed view of this, we show in Fig. 8 the time evolution of the components of the electromagnetic field evaluated at the center of the simulation box. The local electric field components $ e E x / m \omega p c$ (top) and $ e E y / m \omega p c$ (middle) and the magnetic field component $ e B z / m \omega p$ (bottom frame) are displayed therein. The dominant oscillation for all these quantities takes place at saturation at the magnetic bounce frequency *ω _{b}*. The first oscillation due to magnetic trapping, which is related to the bounce frequency of the $ k M = k 0$ mode, is altered by the second harmonic corresponding to the oscillation at about $ 2 \omega b$ of the $ 2 k M = 2 k 0$ mode. This double period appears in the various components of the electromagnetic field, but it is more evident in the

*E*component (top frame of Fig. 8).

_{x}All these results mirror those of Fig. 2 [case (*S1*)], so as the results shown in Fig. 9: in the latter, the spatial profiles of the normalized components $ e E x / m \omega p c$ and $ e B z / m \omega p$ are, respectively, shown on the left and right frames at different times. At time $ t \omega p = 7.5$ (top panels), still in the transient interval, the *E _{x}* fluctuation, coherently with its “nonlinear” nature, has increased from the machine error to an amplitude about the square of the

*B*component (initially excited with amplitude $ e B 0 / m \omega p \u223c 10 \u2212 5$): this can be seen also in the spatial profile, modulated for

_{z}*E*as the first harmonic of the dominant Weibel mode, visible as the fundamental mode in the spatial oscillation of the

_{x}*B*component. This correspondence persists until the beginning of the saturation phase, at time $ t \omega p = 82.5$ (middle panels), when the electrostatic harmonic perturbation begins to steepen. Later, at $ t \omega p = 90$, i.e., at saturation, the appearance of the harmonic 4

_{z}*k*in the

*E*profile indicates the beginning of a wave-breaking process, while the magnetic field is still dominated by the fundamental mode (bottom panels).

_{x}Despite the saturation stage seems to be dominated by magnetic trapping, and in spite of the strong decrease in the electrostatic energy component $ \u03f5 e , x$, the electrostatic activity is not completely negligible in the advanced nonlinear regime, as a certain level of electrostatic fluctuations persists after $ t \omega p \u2273 90$ (cf. Fig. 7, bottom frames). This dynamics seems to be related to the nonlinear “cascade” of the electrostatic field oscillations evidenced in bottom frames of Fig. 9. The role of these electrostatic fluctuations in the saturation stage and in the way they affect the global dynamics of the plasma can be ascertained by looking at the phase-space patterns of the distribution function over time.

In Fig. 10, we have represented the distribution function $ f ( p x , p y )$ at $ t \omega p = 82.5$, i.e., at the beginning of the saturation stage, at three different spatial locations: at $ x = L x / 4$ (upper panel), at the center of the box, $ x = L x / 2$ (middle panel), and at $ x = 3 L x / 4$ (lower panel). These points correspond to maximum, zero, and minimum values of $ B z ( x , t ) | t = 82.5$ (cf. center right panel of Fig. 9).

The deformation of the distribution function in Fig. 10 depends on the sense of rotation (at the cyclotron frequency Ω) imposed by the magnetic field. Note in this regard the perfect symmetry in the momentum space of *f*, when *B _{z}* = 0 (middle panel of Fig. 10).

In Fig. 11, the distribution function $ f ( x , p y )$, averaged over *p _{x}*, is shown at three times in the saturation stage. The

*Y*-pattern of these contour-plots, which had been already observed in Refs. 27 and 52, is due to the alternating sign, in space, of the

*y*-component of the electron velocity (condition required by the spatial modulation of the magnetic field, cf. Fig. 1). In addition to the characteristic

*Y*-pattern, a spatial modulation related to the 2

*k*mode is also present since the beginning of the saturation regime. It is visible in the top panel at $ t \omega p \u2243 79$, although its amplitude is much smaller than that modulated at

*k*. At later times, it becomes dominant by mixing with higher order harmonics, coherently with the behavior of the

*E*component shown in Fig. 9. This means that, despite the dominance of the magnetic trapping mechanism, wave–particle resonance due to electrostatic fluctuations induced by nonlinear coupling become progressively more important, as the evolution progresses in the nonlinear stage.

_{x}The effects of electrostatic fluctuations on the shape of the distribution function become more evident if we look at the contour-plots of $ ( 1 / L y ) \u222b f ( x , p x , p y ) \u2009 d p y$, shown in Fig. 12 at the same times of the frames of Fig. 11. The spatial modulation is in this case dominated by the 2*k* mode, coherently with the ordering $ u x \u2243 u x , 2 \u223c ( n ( x , t ) \u2212 n 0 ) \u223c E x$ previously established [cf. Eqs. (27) and (30)]. The appearance of the 4*k* harmonic is also visible, starting from $ t \omega p = 86$ (cf. also the bottom panel of Fig. 12, at time $ t \omega p = 90$, with the bottom left panel of Fig. 9).

A more detailed picture of the effects induced by the electrostatic field on the distribution function in the $ ( x , p x )$ phase space is obtained by projecting *f* on this plane for specific values of *p _{y}*. This allows us to distinguish the way electrostatic fluctuations differently affect the plasma dynamics, depending on whether one considers the bulk electron plasma or the regions of the distribution function corresponding to a very low particle density (e.g., the tail of the distribution function). Figure 13 shows the contour-plots of $ f ( x , p x )$ at time $ t \omega p = 82.5$ for various “populations” of particles having the same value of the canonical moment, $ C y , N = p y + ( e / c ) A y$ (or $ p y \u2243 C y , N$ if $ | ( e / c ) A y | \u226a | C y , N |$).

For a centered particle distribution, the value *p _{y}* = 0 (top panel of Fig. 13) selects the bulk of the electron plasma. This gives the main contribution to the electron density

*n*(

*x*,

*t*). The spatial modulation at the mode $ 2 k M$ indicates that here electrostatic effects dominate, in agreement with the discussion of Sec. III and with the orderings of Eq. (30). The phase-space dynamics in this region strongly differs from that shown in the other panels, where $ p y = m v y t h$ (center panel) or $ p y = 2 m v y t h$ (bottom panel) have been selected: the more increases the velocity of the electron population, the stronger magnetic effects become, and in the bottom panel the spatial modulation of the Weibel mode

*k*dominates. This implies that magnetic effects and magnetic trapping are expected to dominate in regions of lower particle density, since a larger particle velocity typically means a smaller particle density.

_{M}### B. Reduced kinetic simulations via the multi-stream model

We now compare the results of the numerical simulations of Sec. V A with a simulation performed with the multistream-waterbag code^{63} using the same set of physical parameters (and normalization for quantities) that has been used for the simulation case *S2*. Hereafter we refer to these results as to those of the simulation case (*S3*).

(

*S3*) The simulation has been initialized with an anisotropy in the pressure components equivalent to that of the case (*S2*), i.e., $ P y , 0 = 25 P x , 0 = 25 k e V$, and by perturbing the marginal equilibrium with a sinusoidal perturbation $ B z = B 0 \u2009 sin ( k M x )$ on a spatial domain with extension $ L x = 2 \pi / k M$, where $ k M d e = 2.017$ (corresponding to $ \gamma M = 0.129 \omega p e$) and $ B 0 = 10 \u2212 5$. The equivalent anisotropy has been realized by consistently choosing the parameters*n*for a total number of seven streams ( $ N = \u2212 3 , \u2212 2 , \u2026 , 3$) symmetrically distributed around $ C y 0 = m u y 0 = 0$ according to Eqs. (59) and (71), while a Gaussian distribution with normalized temperature $ k B T x , 0 = 1 \u2009 keV$ (i.e., $ P x , 0 = 1 \u2009 keV$ in normalized units) has been initially chosen for each beam. This means that we have taken $ f N ( x , v x , 0 ) \u221d \u2009 exp [ \u2212 m v x 2 / ( 2 k B T x , 0 ) ]$ in Eq. (46)._{N}For a more direct comparison with the numerical results of (

*S2*), five of the beams have been chosen in accordance with the values of*p*sampled in Fig. 12: $ u y 0 = 0 , \xb1 v y t h , \xb1 2 v y t h , \xb1 3 v y t h$. Thanks to the reduced computational cost of the multi-stream model, a finer phase-space mesh has been used with respect to the simulation (_{y}*S2*), namely, $ N x N v x = 256 \xd7 512$ points. Even in this simulation, the time step has been chosen to be $ \Delta t = 0.005 \u2009 \omega p \u2212 1$.

### 1. Simulation case (*S3*): Electrostatic coupling via magnetic lensing and phase-space dynamics

The results of this simulation globally display an excellent agreement with those of simulation case (*S2*) of Sec. V A 2. Some slight differences can be, however, recognized, which allow further insight in the nature of the nonlinear trapping mechanism and confirm the interpretation previously provided. Here below we focus on the discussion of these latter features.

The correspondence between the full Vlasov and the multi-stream simulation is here evidenced by Fig. 14, which shows, in normalized units, the time evolution of the mode *k _{M}* of the Fourier spectrum of the magnetic field amplitude $ | B z |$, in both linear (upper frame) and semi-logarithmic scale (bottom frame). The behavior, already discussed in relation to Fig. 7 of the simulation case

*S2*, is recovered: after a transient phase, dominated by the appearance of oscillations at the electromagnetic frequency $ \omega e m = \omega p 2 + k M 2 c 2 \u2243 1.72 \omega p$, the magnetic field

*B*grows exponentially with a rate $ \gamma M \u2243 0.140 \omega p$, until nonlinear saturation is reached at $ t \omega p \u2243 80$. Low frequency oscillations appear in the saturation regime, at a frequency close to the magnetic bounce frequency, i.e., at $ \omega b \u2243 0.21 \omega p$, in very good agreement with the results of the simulation (

_{z}*S2*).

For further comparisons between the full-Vlasov and multi-stream simulations of Weibel-type modes, we address the interested reader to look at previously published works (see, e.g., Refs. 37, 41, 43, 52, and 63). We now turn the attention to more specific features of the coupling between electrostatic and e.m. modes, which become even more evident in the multi-stream modeling.

Figure 15 must be compared with Fig. 13, of which it represents the same contour-plots (i.e., the contours of the distribution of electrons having *p _{y}* = 0, $ p y = m v y t h$, and $ p y = 2 m v y t h$ in the $ ( x , p x )$-plane), taken at equivalent times to those of the simulation case (

*S2*): the difference between the spatial modulation of the central stream (

*p*= 0, for which the electrostatic mode 2

_{y}*k*dominates) and that of the two other streams (for which the magnetic Weibel mode

*k*dominates) is in this case accentuated. In particular, in the multi-stream model, the mixing of electrostatic and of magnetic trapping effects is mitigated for the electron population with $ p y = m v y t h$ (compare the center frames of Figs. 13 and 15), for which magnetic trapping more evidently dominates. This is coherent with the arguments developed in Secs. V A and V B and with the fact that in the multi-stream modeling the dynamics of populations with different values of $ C y , N \u2243 p y$ are coupled one with each other more weakly than in the full Vlasov-Maxwell mode. In particular, the dynamics of the central beam is dominated by the $ 2 k M$ mode (Fig. 15, top), i.e., by electrostatic (second-order) effects, while the dynamics of the beams positioned at $ p y = p y , t h$ and $ p y = 2 p y , t h$ is more directly linked to magnetic trapping at the dominant $ k M = k 0$ mode.

Figure 16 displays the time evolution of the Fourier component $ 2 k M$ of the electrostatic field *E _{x}* (top frame) and the Fourier component

*k*of both the inductive electric field

_{M}*E*(middle frame) and

_{y}*A*vector potential (bottom frame). These results further support the interpretation given in Fig. 15 and agree very well with the results of simulation case

_{y}*S2*. The component $ | E \u0302 x ( 2 k M ) |$ reaches its maximum amplitude at a time similar to that of $ | E \u0302 y ( k M ) |$, and its growth is about $ 2 \gamma M$, i.e., twice the growth of the electromagnetic components $ | E \u0302 y ( k M ) |$ and $ | A \u0302 y ( k M ) |$. At saturation, the $ | E \u0302 y ( k M ) |$ component oscillates at about twice the bounce frequency

*ω*because of the fact that its average time value is zero. The apparent doubling of period follows from having evaluated the (Fourier transform) of the modulus of

_{b}*E*.

_{y}^{49}However, $ \u27e8 \Gamma \u27e9 \u2243 1$ in the non-relativistic limit considered in this simulation case. Using then the definitions (52) and (54), one verifies that at

*t*= 0 $ \u27e8 \Pi \u0303 x x \u27e9 L x \u2243 \u222b \Pi x x tot d x / L x$ [cf. Eq. (56)]. This provides an excellent approximation of $ \u222b \Pi x x tot d x / L x$ also at any subsequent time, as long as the initial symmetry of the beams is respected during their evolution. We thus see that at the beginning of the saturation phase, $ \u27e8 \Pi \u0303 x x \u27e9 L x$ behaves like the electric field

*E*, i.e., undergoes an amplification at time $ t \omega p = 60$, with a growth rate $ 2 \gamma B$ (cf. Sec. IV). Its saturation regime is, however, dominated by magnetic trapping, with a temporal modulation at low frequency close to the bouncing frequency

_{x}*ω*, the same of the saturation stage of the beams with initial non-zero velocity.

_{b}## VI. SUMMARY AND CONCLUSIONS

### A. Summary

In this work, we have revised the mechanism of linear and nonlinear coupling between a purely e.m., non-relativistic Weibel-type mode, and electrostatic fluctuations. From the theoretical point of view, we addressed the problem by making reference to the linear analysis of the anisotropy-driven Weibel instability using different reduced models, which we have revised and compared (Sec. IV) by extending, although in the non-relativistic limit only, the discussion earlier developed in Ref. 41. From a numerical point of view, we studied the nonlinear evolution of a Weibel instability using kinetic simulations (Sec. V) based on both the 1D-2V version of the semi-Lagrangian Vlasov–Maxwell solver VLEM^{43} (Sec. V A) and a numerical code implementing the multi-stream model^{37} (Sec. V B). The mechanism at the basis of the nonlinear coupling between e.m. and electrostatic fluctuations is here interpreted (Sec. III and Appendix C) as due to the spatial concentration of electron beams, which non-uniformly stream orthogonally to the wave-vector perturbation, depending on the sign and intensity of the spatial gradient of the local magnetic field. This electron density fluctuation over a uniform ion background induces space charge fluctuations. Due to the large timescale separation between the response of ions and electrons to electrostatic fluctuations, proportional to $ m i / m e$, this effect persists at least in the linear regime of the instability even when the ion dynamics—here neglected for simplicity—is self-consistently accounted for.

The main outcome of this analysis is the identification of a combination of two complementary but simultaneous features. On the one side, fluctuations of all the electron charge density, the electrostatic field, and the variation of the fluid electron velocity orthogonal to the initial beams (to which the time variation of electron density is related via fluid compressibility in the continuity equation) are of second order with respect to the e.m. Weibel mode perturbations; in this sense, these are “nonlinear” effects. This is consistent with what Taggart *et al.*^{26} and following authors^{27,29} already noted. Conversely, these second-order fluctuations represent an intrinsic feature of the *linear* evolution of the Weibel instability, since they correspond to the transversal separation of initially superposed beams and the concentration of the electron density current, both implied and required by the instability mechanism. This conclusion thus contributes to the debate present in the literature, about some delicate issues that have been already pointed out in previous works,^{1–11} in relation to the formal application of the linearization procedure for small amplitude modes: the Weibel instability provides an example in which electrostatic and thermal fluctuations, being both of compressive nature, intervene as effects of order $ \epsilon 2$ with respect to the amplitude *ε* of the linear e.m. Weibel mode. At the same time, a self-consistent description of the instability [cf. Eq. (29)] makes their evolution not separable from the growth of the linear mode, in spite/regardless of the applicability of the formal limit $ \epsilon 2 \u226a \epsilon $.

These features allow further insight in the interpretation of the mechanism of Weibel instability. In Sec. IV A 1, it has been discussed from the point of view of generalized Ohm's law in a fluid description. In this light, the instability condition $ P \u22a5 k > P | | k$ has been interpreted in Sec. IV A 2 in terms of the single particle dynamics, by comparing it to the fluid description.

Density and magnetic fluctuations are differently coupled with thermal fluctuations of the plasma, i.e., the components of the pressure tensor during the linear evolution of the Weibel instability (Sec. IV A). The analytical estimates agree with the fully kinetic results (Sec. V A): while the off-diagonal fluctuation (i.e., the $ \delta \Pi x y$ component) is linear, and thus proportional to the magnetic field amplitude, the consistency of the linear analysis requires $ \delta \Pi x x$ and $ \delta \Pi y y$ to be of second order with respect to the magnetic field amplitude (we recall that $ k = k e x$). An anisotropic response of thermal fluctuations to harmonic perturbations of coupled e.m./compressive type was already noted for other e.m. modes in collisionless plasmas.^{55}

As the Weibel mode enters its nonlinear stage, the differential rotation of the current filaments around the extrema of the magnetic field, which in the linear stage is responsible of the charge density fluctuations, gives rise to the well-known magnetic trapping structures (cf. Sec. IV C 1) that lead to the saturation of the e.m. instability. Numerical simulations run with a semi-Lagrangian multi-stream model (Sec. V B) allow us to single out the effect that electron density fluctuations, and thus electrostatic perturbations, have on the evolution of different “bunches” of electrons, depending on their average fluid velocity. This effect qualitatively agrees with the predictions of the toy-model of Sec. III. However, it is somewhat “enhanced” with respect to full Vlasov–Maxwell simulations (Sec. V A), in which it is also visible, but less evident. The reason is that in the multi-stream approach, streams with different values of the invariant momentum component are made evolve according to separate Vlasov equations coupled only by Maxwell equations (Sec. IV B), and their coupling is, therefore, weaker in than in full Vlasov–Maxwell simulations.

### B. Conclusions

The results discussed in this work provide a complementary interpretation to the occurrence of electrostatic fluctuations during the onset of linear Weibel-type instabilities, as they have been already discussed in the literature.^{26–35}

The analysis we have carried out may be of interest for the identification of energy coupling mechanisms during the turbulence cascade in collisionless plasmas, in which an anisotropic heating has been shown to play a key role.^{18–22} Of course, the results we have discussed would be just preliminary to this end, since they are limited to the dynamics of a single unstable mode and to a single propagation direction, conditions too reductive for a comparison with turbulence. Nevertheless, looking at the way small amplitude fluctuations balance each other in an extended fluid model including the full pressure tensor evolution, like we here did for Weibel's instability (cf. Sec. IV A), may be of some usefulness for the interpretation of kinetic simulations, especially at a local scale (e.g., near coherent structures like current and vorticity sheets, whose strong gradients locally break the spatial isotropy of the medium) or when the turbulent cascade results to be not isotropic with respect to the wave-vector orientation. We note in this regard that simulations of solar wind turbulence,^{18,69–71} where Weibel-type instabilities coupled with other beam–plasma instabilities may be in principle expected to also play a role^{16} close to the vorticity sheets separating vortices, have shown that electron pressure anisotropization and electrostatic fluctuations become progressively more important as the electron scales are approached.

Finally, it is interesting to extend the analysis here performed to the CFI instability, especially in relativistic regimes, which are relevant to both laser–plasma interactions and astrophysical shocks in novae. In this context, the coupling between electrostatic perturbations and thermal effects has been already discussed,^{34,72–78} although in somewhat different terms with respect to those considered in the present work. The extension of the present analysis to relativistic regimes requires, however, more effort, since the inclusion of relativistic effects makes the modeling more complex. On the one hand, a non-negligible ratio between the velocity of the beams and that of light is known^{46} to induce a linear coupling between electrostatic and e.m. modes, also for wave-vectors purely perpendicular to the beams; moreover, even for non-relativistic beams, an imbalance in the flux along the beams of the temperatures perpendicular to them produces the same effect^{39} (cf. also Refs. 74 and 76 for the relativistic case). On the other hand, the possible interpretation of the coupling between electrostatic and thermal fluctuations in fluid terms is complicated by the need to connect a covariant definition of a pressure tensor, possibly anisotropic also in the rest frame of reference, with the relativistic motion associated with beams of electrons streaming with different velocities. The relativistic multi-stream approach discussed in Refs. 37 and 49, and of which in this work we discussed the non-relativistic limit, seems to be a promising theoretical tool for the further extension of this kind of analysis to the relativistic anisotropy-driven Weibel instability^{66,79,80} and to the its time-resonant counterpart,^{40,68,81–83} and finally to relativistic oblique modes, in which a richer interplay between magnetic, electrostatic, and thermal effects has been pointed out.^{34,77,84–88}

## ACKNOWLEDGMENTS

We are grateful to two anonymous referees for their interesting comments and careful reading. In particular, we owe to one of them the suggestion that the electrostatic field discussed in Sec. III may be also interpreted as the one making electrons drift with respect to the self-generated magnetic field. Elaborating this argument led us to develop Sec. IV A 2, added during the peer-review. We thank the IDRIS computational center in Orsay, France for computer time allocation granted via GENCI under project no. A0110504028. This work has been partially carried out within the framework of the French Federation for Magnetic Fusion Studies (FR-FCM) and of the Eurofusion consortium (funding obtained through No. FR-FCM AAP 2023, “Evolution of current sheets in low-collision plasmas,” in particular are gratefully acknowledged). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**D. Del Sarto:** Conceptualization (lead); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (supporting); Validation (equal); Visualization (supporting); Writing – original draft (lead); Writing – review & editing (equal). **A. Ghizzo:** Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Visualization (lead); Writing – original draft (supporting); Writing – review & editing (equal). **M. Sarrat:** Conceptualization (supporting); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (supporting); Validation (equal); Visualization (lead); Writing – original draft (supporting); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: LINEAR KINETIC RESULTS ABOUT WEIBEL INSTABILITY

^{25}to the linear system $ [ \epsilon ] \xb7 E 1 = 0$, where the dispersion matrix $ [ \epsilon ]$ has a diagonal form with components

^{89}of argument $ \xi = \omega / ( 2 k c x )$. The condition $ \epsilon x x = 0$, associated with perturbations on $ E x , 1$, clearly describes electrostatic, Bohm–Gross type modes. The dispersion relation of two decoupled electromagnetic (e.m.) modes is instead described by the two equations $ \epsilon y y = 0$ and $ \epsilon z z = 0$. Perturbations on $ E y , 1$ are related to magnetic perturbations $ B 1 , z$, which are exponentially growing (Weibel instability) if $ P y , 0 > P x , 0$, whereas perturbations on $ E z , 1$ correspond to a Weibel instability that amplifies a magnetic field $ B 1 , y$ if $ P z , 0 > P x , 0$. The dispersion relation of the Weibel mode triggered by the pressure anisotropy in the (

*x*,

*y*) plane, is $ \epsilon y y = 0$, which clearly does not depend on the $ g ( v z ; P o , z )$ part of the distribution function (1). In the “hydrodynamic limit,” which corresponds to the condition of negligible thermal effect expressed by the parameter expansion $ k 2 c x 2 / \omega 2 \u226a 1$ (or equivalently to $ \xi \u226b 1$), the dispersion relation of Weibel modes (A2) related to a velocity anisotropy in the (

*x*,

*y*) plane becomes

*P*>

_{y}*P*is evident from the diagonal form of $ [ \epsilon ]$.

_{x}### APPENDIX B: LINEAR ANALYSIS AND ELECTROSTATIC COUPLING IN THE COLD CFI MODEL

^{36}of Weibel instability consisted in replacing the initial bi-Maxwellian distribution function with Eq. (3). This condition can be exactly described in a cold fluid model, which provides a quantitatively appropriate description of the CFI,

^{46}as long as particle kinetic heating processes are negligible. These become important for large values of the wave-vector, where the cold fluid model fails to provide a cutoff wavenumber. It is here interesting to recall the essential points of this modeling. Although the focus of this work is on the non-relativistic case, it is instructive to consider the fluid modeling of the cold relativistic CFI described by the system ( $ \alpha = a , b$),

*x*and wave-vector along

*y*.

It is important to notice that $ \omega 3 2 = 0$ regardless of the value of $ \Gamma \alpha , 0$ for symmetric beams, that is, when $ n b , 0 = n b , p$, and that $ \omega 3 2 \u2192 [ 4 \pi e 2 / ( m c 2 ) ] \u2211 \alpha ( 1 \u2212 \Gamma \alpha , 0 ) ( n \alpha , 0 u \alpha , 0 )$ as $ \Gamma \alpha , 0 \u2192 1$ for $ ( u \alpha , 0 / c ) 2 \u2192 0$, which means, again, that $ \omega 3 2 = 0$ in the non-relativistic limit where $ \Gamma \alpha , 0 = 1$.

The linear decoupling of the pure CFI, e.m. mode (i.e., perturbations on the $ E y , 1$ component) from the electrostatic mode (i.e., perturbations on the $ E x , 1$ component), is, therefore, exact in the limits where $ \omega 3 2 = 0$.

^{38}or, equivalently, to the growth of current densities parallel to the beams.

^{27,29}Both these mechanisms can be microscopically interpreted as due to the charge density perturbation caused by the magnetic lensing described in Sec. III.

### APPENDIX C: LINEAR ANALYSIS AND ELECTROSTATIC COUPLING IN THE COLD CFI MODEL

We here detail the algebra leading to the conclusion discussed in Sec. III.

Naming *B _{M}* the reference value of the magnetic field, which we choose to correspond to its maximum amplitude at saturation, the growing field related to the CFI reads $ B z ( x , t ) = B M \epsilon ( t ) \u2009 cos ( k x x )$, with $ \epsilon ( t ) = \epsilon 0 e \gamma B t$. Saturation corresponds to $ \epsilon ( t ) \u223c O ( 1 )$. All scalar functions

*F*(

*x*,

*t*) are assumed to be separable with respect to space and time, as $ F ( x , t ) = X F ( x ) T F ( t )$.

In a volume $ V 0 = L x \xd7 L y \xd7 L z$, in which ions form a fixed neutralizing background, we then consider *N* electrons that would be at rest at *t* = 0, if $ \epsilon | t = 0 = 0$, and for which $ n ( x , t ) = N ( x , t ) / V 0$ with $ n 0 = n ( x , 0 ) = N / V 0$. When, instead, $ \epsilon ( t ) \u2260 0$, the perturbation $ \delta B ( x , t ) | t = 0 = ( 0 , \u2009 0 , \u2009 B z ( x , 0 ) )$ implies through Ampére's law [second of Eq. (9)] the existence of $ \delta J = ( 0 , J y , 0 )$ given by Eq. (C1).

*x*, used previously), Newton equations of the single electron read

**, and thus between the e.m. fields in Newton's law and the particle population: writing $ v \u0303 x = \delta v \u0303 x , \u2009 v \u0303 y = \delta v \u0303 y , \u2009 n = n 0 + \delta n$, and $ J = \delta J , \u2009 B = \delta B , \u2009 E = \delta E$, where all quantities with**

*J**δ*in front would be zero if $ B M = 0$, we thus look for approximated solutions of Newton equations. These solutions should be valid as long as $ max { \epsilon , \epsilon \u0307 \Delta t , \epsilon \u0308 ( \Delta t ) 2} \u226a 1$ (i.e., for $ \Delta t \u226a min { \epsilon 0 / \gamma B , \u2009 \epsilon 0 / \gamma B}$). In terms of the separation of variables used before, any scalar function

*F*(

*x*,

*t*) has $ X F ( x ) \u223c O ( 1 )$ while $ T F ( t )$ can be ordered like a positive power of

*ε*.

**, statistically defined by Eq. (6), as $ \delta J \u2243 \u2212 e n 0 \delta v \u0303 + O ( \delta n \delta v \u0303 )$, which is consistent with the orderings (13) and (C8). This allows us to close Eq. (C10) by using the continuity equation of charge density, as usual obtained from the combination of the first of Eqs. (9) and (10). At the leading order (albeit unknown, yet), the continuity equation can be formally written as**

*J*## REFERENCES

*q*in a magnetic field associated with a vector potential $\nA\n(\nx\n,\nt\n)$ in the absence of collisions is established in a phase-space where the canonical conjugated variables are $\n(\nx\n,\nP\n)$. The passage from the canonical momentum $\nP\n=\np\n+\nq\nA\n(\nx\n,\nt\n)$ to the kinematic momentum

**that allows the immediate rewriting of Vlasov equation for $\n\nf\n\xaf\n(\nx\n,\nP\n,\nt\n)$ as a Vlasov equation for $\nf\n(\nx\n,\np\n,\nt\n)$ in the “phase-space” $\n(\nx\n,\np\n)$ (or for $\nf\n(\nx\n,\nv\n,\nt\n)$) in the $\n(\nx\n,\nv\n)$ configuration space for the non-relativistic limit $\np\n=\nm\nv$) is allowed by the equivalence of the infinitesimal volumes $\nd\nx\nd\nP\n=\nd\nx\nd\np$, which grants the validity of Liouville's theorem when performing the change of variables $\nP\n\u2192\np$ (cf., e.g., Ref. 45, §2.2.1.3).**

*p**The Vlasov Equation 1. History and General Properties*

*The Plasma Dispersion Function. The Hilbert Transform of the Gaussian*