# Nonlinear coupling of electromagnetic and electrostatic modes via density and pressure fluctuations: The case of Weibel instabilities

D. Del Sarto, A. Ghizzo, M. Sarrat

► To cite this version:

D. Del Sarto, A. Ghizzo, M. Sarrat. Nonlinear coupling of electromagnetic and electrostatic modes via density and pressure fluctuations: The case of Weibel instabilities. Physics of Plasmas, 2024, 31 (7), pp.072113. 10.1063/5.0207974. hal-04661965v2

HAL Id: hal-04661965

https://hal.science/hal-04661965v2

Submitted on 25 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Nonlinear coupling of electromagnetic and electrostatic modes via density and pressure fluctuations: The case of Weibel instabilities 🌐

D. Del Sarto 🧠 ; A. Ghizzo 🧠 ; M. Sarrat 🧠

<figure><img src="image_2.png" /></figure>

Check for updates

Phys. Plasmas 31, 072113 (2024)

https://doi.org/10.1063/5.0207974

<figure><img src="image_3.png" /></figure>

View
Online

<figure><img src="image_4.png" /></figure>

Export
Citation

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>
# Nonlinear coupling of electromagnetic and electrostatic modes via density and pressure fluctuations: The case of Weibel instabilities F

Cite as: Phys. Plasmas **31**, 072113 (2024); doi: 10.1063/5.0207974

Submitted: 12 March 2024 · Accepted: 29 June 2024 ·

Published Online: 25 July 2024

<figure><img src="image_8.png" /></figure>

View Online

<figure><img src="image_9.png" /></figure>

Export Citation

<figure><img src="image_10.png" /></figure>

CrossMark

D. Del Sarto,ª ID A. Ghizzo,ᵇ ID and M. Sarratᶜ ID

## AFFILIATIONS

Institut Jean Lamour, UMR 7198 CNRS – Université de Lorraine, F-54000 Nancy, France

## ABSTRACT

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.

© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution-NonCommercial 4.0 International (CC BY-NC) license (https://creativecommons.org/licenses/by-nc/4.0/). https://doi.org/10.1063/5.0207974

## 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¹⁻¹¹ 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,¹⁵⁻¹⁷ has been identified¹⁸⁻²² 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²⁵ 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

<footnote>ª) Author to whom correspondence should be addressed: daniele.del-sarto@univ-lorraine.fr</footnote>

<footnote>ᵇ) Electronic mail: alain.ghizzo@univ-lorraine.fr</footnote>

<footnote>ᶜ) Electronic mail: sarrat.mathieu@gmail.com</footnote>
nonlinear kinetic evolution of a classical Weibel instability in a one-
dimensional (1D) non-uniform geometry, in which electrostatic fluctu-
ations are linearly decoupled from e.m. modes. We focus on the char-
acterization of the electrostatic mode that is observed to grow,²⁶ since
the beginning of the linear stage, as the second harmonic of the unsta-
ble mode—we recall that the (“pure,” or transverse) Weibel mode cor-
responds to an exponential growth, with rate γ<sub>B</sub> 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,²⁶⁻²⁹ both on the
anisotropy-driven Weibel instability and on the current filamentation
instability³⁰⁻³⁵ (CFI), strictly related to it. The CFI can be indeed made
correspond³⁶,³⁷ to a limit case of the Weibel instability, where the tem-
perature anisotropy associated with a double-Maxwellian velocity dis-
tribution is replaced by two counter-propagating cold electron beams
of initial densities $n_1$ and $n_2$.

Taggart et al.²⁶ were probably the first ones to identify this elec-
trostatic growth as a second-order effect in the kinetic linear analysis
of Weibel instabilities, by noting that it results from the nonlinear cou-
pling 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³⁸ to develop in the numerical integration of the CFI. An expo-
nential growth of the electrostatic field during the linear stage of the
CFI, with a rate larger than $\gamma_B$, has been indeed observed and discussed
in several works about the CFI.³⁰⁻³⁵ A quantitative explanation of the
rapid electrostatic mode growth—larger than $\gamma_B$ 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 $\gamma_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 ($n_1 = n_2$), where an
exact decoupling between electrostatic and Weibel-type modes is pos-
sible (as pointed out, e.g., also in Ref. 34): this effect is, therefore, differ-
ent 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 simu-
lations 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 nonlin-
ear 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:³⁹ this is a fitting approximation for the description of the linear evolution of the Weibel instability³⁹⁻⁴¹ (Sec. IV A). We then revise (Sec. IV B) the key features of the reduced multi-stream kinetic model,⁴² specialized to describe the Weibel instability;³⁷ 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,⁴³ 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.

In the case of Weibel instability, the electron dynamics is
described by the Vlasov equation, coupled to the Maxwell equations.
Assuming an initial pressure anisotropy in the (x, y)-plane, with
$P_{v,0} > P_{x,0}$ and a perturbation with $k = ke_x$, the dispersion relation of Weibel instability²⁵ is obtained by linearizing this system around the neutral equilibrium distribution function for electrons,

$$f_0 = \frac{mn_0}{2\pi} \left( \frac{1}{P_{x,0}P_{y,0}} \right)^{\frac{1}{2}} g(v_z; P_{z,0}) e^{-n_0 \left( \frac{m v_z^2}{2P_{x,0}} + \frac{m v_y^2}{2P_{y,0}} \right)}, \quad (1)$$

where $g(v_z; P_{z,0}) = (mn_0)^{1/2}/(2\pi P_{z,0})^{1/2} e^{-m v_z^2/(2P_{z,0})}$, $m$ is the electron mass, and $P_{i,0} = n_0 k_b T_i$, for $i = x, y, z$. From now on, labels 0 and 1 will be used for equilibrium and first-order perturbations, respectively.

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.

CFI can be seen as a simplified version of Weibel instability: as shown by Fried,³⁶ 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,

$$\mathbf{\Pi}_0 = \begin{pmatrix} P_{x,0} & 0 & 0 \\ 0 & P_{y,0} & 0 \\ 0 & 0 & P_{z,0} \end{pmatrix}, \quad (2)$$

chosen so that $J_y = 0$, initially.

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⁴⁴ 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⁴⁶ 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³⁹,⁴⁷ 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.

Re-writing the non-relativistic phase-space electron distribution function $\tilde{f}(x, P, t)$, where $P \equiv p - eA(x, t)$ is the electron canonical

<figure><img src="image_11.png" /><figcaption>FIG. 1. Mechanism of magnetic field amplification in Weibel instability sketched from top to bottom. Top frame: Two counter-propagating cold electron beams (along $y$) are perturbed by a magnetic field $\delta B = B_{z,1}e_z$ modulated with respect to $x$ (red signs). Center frame: The magnetic field deflects electrons via Lorentz force, which concentrates electrons in regions where the modulus of the resulting magnetic field is weaker and empties regions where it is stronger. A current density fluctuation $\delta J = \delta J_y e_y$ (blue arrows), spatially modulated by the inhomogeneity of the magnetic field profile, is this way locally generated. Bottom frame: The current density induces a magnetic field (blue signs) which reinforces, in turn, the initial magnetic perturbation. The process goes on at the expenses of the initial electron velocity along $y$, which is reduced until saturation by the $\delta J_z$ fluctuation so generated (red arrows). Adapted from Ref. 45, reproduced with permission from the author.</figcaption></figure>
momentum,¹⁸ in the $(x, v)$ configuration space as $f(x, v, t)$, Vlasov equation reads

$$ \frac{\partial f}{\partial t} + v_x \frac{\partial f}{\partial x} - \frac{e}{m} \left( E_x + \frac{v_y}{c} B_z \right) \frac{\partial f}{\partial v_x} = 0. \quad (4) $$

Assuming ions at rest, Eq. (4) implies that the equation for the total energy density of the electrons, $\mathcal{E} = (m/2) \int v^2 f(x, v, t) d^3 v$, fulfills

$$ \frac{\partial \mathcal{E}}{\partial t} + \frac{\partial}{\partial x} \int \frac{mv^2}{2} v_x f(x, v, t) d^3 v = \mathbf{E} \cdot \mathbf{J}, \quad (5) $$

where

$$ \mathbf{J}(x, t) \equiv -e \int v f(x, v, t) d^3 v. \quad (6) $$

This can be combined with the Poynting theorem, here expressible in the form

$$ \mathbf{E} \cdot \mathbf{J} = - \frac{\partial}{\partial t} \left( \frac{E^2}{8\pi} + \frac{B^2}{8\pi} \right) - \frac{c}{4\pi} \frac{\partial}{\partial x} (E_y B_z), \quad (7) $$

due to the form that Faraday's and Ampère's laws take in this geometry,

$$ \frac{\partial B_z}{\partial t} = -c \frac{\partial E_y}{\partial x}, \quad (8) $$

$$ 4\pi J_x + \frac{\partial E_x}{\partial t} = 0, \quad 4\pi J_y + \frac{\partial E_y}{\partial t} = -c \frac{\partial B_z}{\partial x}. \quad (9) $$

The particle distribution function indirectly contributes to the electrostatic field component via the definition of spatial density $n(x, t) \equiv \int f(x, v, t) d^3 v$ in the Gauss equation, where we assume an (hydrogen) ion density $n_0 = n(0, 0)$,

$$ \frac{\partial E_x}{\partial x} \equiv 4\pi e(n_0 - n). \quad (10) $$

From Eqs. (5) and (7) follows the equation of the total energy density in the usual continuity form,

$$ \frac{\partial}{\partial t} \left( \mathcal{E} + \underbrace{\frac{E_x^2}{8\pi}}_{\equiv \epsilon_{e,x}} + \underbrace{\frac{E_y^2}{8\pi}}_{\equiv \epsilon_{e,y}} + \underbrace{\frac{B_z^2}{8\pi}}_{\equiv \epsilon_{m,z}} \right) = -\frac{\partial}{\partial x} \left( \frac{c}{4\pi} E_y B_z + \int \frac{mv^2}{2} v_x f(x, v_x, v_y, t) d^3 v \right). \quad (11) $$

The vanishing of the right-hand side term after integration over the spatial domain shows that the e.m. energy components can grow only at the expenses of the $\mathcal{E}$ energy component. The latter consists both of the bulk kinetic energy $K \equiv m \int (\mathbf{u}^2/2) d^3 v$, where $\mathbf{u} \equiv (1/n) \int v f d^3 v$, and of the thermal energy components $\mathcal{E}_{T,i} \equiv (m/2) \int (v_i - u_i)^2 f d^3 v$ for $i = x, y, z$, according to

$$ \mathcal{E} = K + \mathcal{E}_T = K + \mathcal{E}_{T,x} + \mathcal{E}_{T,y} + \mathcal{E}_{T,z}. \quad (12) $$

Introducing the pressure tensor components $\Pi_{ij} \equiv m \int (v_i - u_i) \times (v_j - u_j) f d^3 v$, one sees that $\mathcal{E}_T = \text{tr}(\mathbf{\Pi})/2$. Therefore, $K=0$ in both Fried's kinetic version of the CFI and in Weibel instability, since $\mathbf{u} = 0$. Note that in the kinetic description of CFI, $\mathbf{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²⁵,³⁶ that the source of energy in both the kinetic CFI and the Weibel instability is in the difference between $\mathcal{E}_{T,x}$ and $\mathcal{E}_{T,y}$ (or $\mathcal{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,⁴² which we will discuss later, by means of which a formal correspondence³⁷,⁴⁹ 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 \dot{e}(t) \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_x$ of the perturbation. In this writing, the spatial and temporal dependence of $B_z(x, t)$ have been split: the function $\varepsilon(t) = \varepsilon_0 e^{i\gamma t}$ expresses the time dependence of the magnetic field amplitude that exponentially increases with growth rate $\gamma_B$. In the comparison with Weibel modes, $\varepsilon(t)$ can be read as a dimensionless parameter weighting the relative amplitude of the magnetic field perturbation at each time $t$, until saturation is reached when $\varepsilon(t) \sim 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_x$, 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 $N$ electrons, which at $t=0$ would be uniformly distributed over a volume $V_0 = L_x \times L_y \times 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, 0, 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
"..." we will use for the single electron velocity $\tilde{v}(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 $\tilde{v}_x = \delta\tilde{v}_x$, $\tilde{v}_y = \delta\tilde{v}_y$, $n = n_0 + \delta n = (N(t=0) + \delta N)/V_0$, and $J = \delta J$, $B = \delta B$, $E = \delta E$, where the quantities with $\delta$ 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.

The calculations, detailed in Appendix C, are valid for small enough values of all among $\epsilon, \dot{\epsilon}, \ddot{\epsilon} \ll 1$, which corresponds to time intervals $\Delta t \ll \min\{\epsilon_0/\gamma_B, \sqrt{\epsilon_0/\gamma_B}\})$. They yield the orderings,

$$ \delta B_z \sim \delta E_y \sim \delta J_y \sim \epsilon, \quad \delta \tilde{v}_y \sim \epsilon, \quad (13) $$

$$ \delta \tilde{v}_x \sim \delta J_x \sim \delta E_x \sim \delta n \sim \epsilon^2. \quad (14) $$

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(2k_x x)$. This is a second-order effect with respect to the parameter $\epsilon \ll 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 \gg T_x$ reads $\gamma_W \simeq k_x c_y / \sqrt{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 \sim -n_0 \epsilon \delta \tilde{v}_y$ implied by $\delta n \sim \epsilon^2$ is used: the result $\gamma_B = \omega_{pe} \sqrt{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 \equiv c^2 / \omega_{pe}^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 \equiv k_x^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.⁵³

In this model, the fluid equations to be coupled with Maxwell equations (8)–(10) read

$$ \frac{\partial n}{\partial t} + \nabla \cdot (n\mathbf{u}) = 0, \qquad (15) $$

$$ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{e}{m} \left( \mathbf{E} + \frac{\mathbf{u}}{c} \times \mathbf{B} \right) - \frac{\nabla \cdot \mathbf{\Pi}}{nm}, \qquad (16) $$

$$ \frac{\partial \mathbf{\Pi}}{\partial t} + \nabla \cdot (\mathbf{u}\mathbf{\Pi}) + \nabla \mathbf{u} \cdot \mathbf{\Pi} + (\nabla \mathbf{u} \cdot \mathbf{\Pi})^T = - \underbrace{\nabla \cdot \mathbf{Q}}_{=0 \text{ (by Hp.)}} - \frac{e}{mc} (\mathbf{\Pi} \times \mathbf{B} + (\mathbf{\Pi} \times \mathbf{B})^T). \qquad (17) $$

The validity of the model is restricted to non-relativistic electron fluid velocities and temperatures, but relativistic phase-space velocities are admitted, thanks to the fact that electron density fluctuations are retained. These violate quasi-neutrality because of (10), differently from what happens in the incompressible electron-magnetohydrodynamic (EMHD) limit.⁵⁴ 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)].

Linearization around a uniform equilibrium with zero fluid velocity and an initial diagonal pressure tensor given by (2) leads to a diagonal dispersion matrix³⁹ [D] with

$$ D_{yy} = \frac{\omega^2}{c^2} - k^2 - \frac{1}{d_e^2} \frac{\omega^2 + k^2(c_y^2 - c_x^2)}{\omega^2 - k^2 c_x^2}. \qquad (18) $$

This dispersion relation gives also the correct value of the cutoff wave-number of Weibel instability, as it is evaluated from the kinetic Vlasov-Maxwell model, even for values $k^2c_x^2/\omega^2 \ge 10^2$ (cf. Fig. 1 of Ref. 39).

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

It is interesting, for the interpretation of numerical results, to highlight the role played by the components of Eqs. (15)–(17) in the linear analysis, and in particular, in the contribution to the Weibel mode given by $D_{yy} = 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}$:

$$ \frac{\partial n_1}{\partial t} + n_0 \frac{\partial u_{x,1}}{\partial x} = 0, \qquad (19) $$

$$ \frac{\partial u_{x,1}}{\partial t} = -\frac{e}{m} E_{x,1} - \frac{1}{mn_0} \frac{\partial \Pi_{xx,1}}{\partial x}, \qquad (20) $$

$$ \frac{\partial u_{y,1}}{\partial t} = -\frac{e}{m} E_{y,1} - \frac{1}{mn_0} \frac{\partial \Pi_{xy,1}}{\partial x}. \qquad (21) $$

Only the components of $\mathbf{\Pi}$ related to the x, y coordinates contribute to the electron fluid dynamics in the (x, y) plane, via the coupled equations,

$$ \frac{\partial \Pi_{xx,1}}{\partial t} + 3P_{x,0} \frac{\partial u_{x,1}}{\partial x} = 0, \qquad (22) $$

$$ \frac{\partial \Pi_{yy,1}}{\partial t} + P_{y,0} \frac{\partial u_{x,1}}{\partial x} = 0, \qquad (23) $$

$$ \frac{\partial \Pi_{xy,1}}{\partial t} + P_{x,0} \frac{\partial u_{y,1}}{\partial x} = -\frac{e}{m} (P_{y,0} - P_{x,0}) B_{z,1}. \qquad (24) $$

It is worth insisting on the importance of the last equation (and, by extension, to the importance of the off-diagonal component of the pressure tensor) in order to correctly describe the linear mechanism of Weibel instability. The right-hand-side term of Eq. (24) implies the coupling between the initial pressure anisotropy and the magnetic perturbation $B_{z,1}$.

The current filamentation mechanism is driven by the spatial gradient of the $u_{y,1}$ velocity component in Eq. (24). The component $u_{y,1}$ is linearly coupled to the amplification of the magnetic field $B_{z,1}$ and, nonlinearly, to the growth of electrostatic fluctuations, associated with $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 $\mathbf{\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,

$$ \frac{d}{dt}(\Pi_{xx}n^{-3}) = 0 \quad \Rightarrow \quad \Pi_{xx,1} n_0^{-3} + P_{x,0} n_1^{-3} = cnst_1, \qquad (25) $$

$$ \frac{d}{dt}(\Pi_{yy}n^{-1}) = 0 \quad \Rightarrow \quad \Pi_{yy,1} n_0^{-1} + P_{y,0} n_1^{-1} = cnst_2. \qquad (26) $$

The single particle dynamics outlined in Sec. III can be thus recovered in this extended fluid analysis and yields non-trivial insight about the non-isotropic way by which the components of the tensor pressure are involved in the process: naming $\varepsilon$ the “instantaneous” characteristic amplitude of the Weibel mode during its linear stage, we can consistently order the contributions of the fluid momentum equation as

$$ \underbrace{\frac{\partial u_x}{\partial t}}_{\sim \varepsilon^2} + \underbrace{u_x \frac{\partial u_x}{\partial x}}_{\sim \varepsilon^2 \times \varepsilon^2} = -\frac{e}{m} \left[ \underbrace{E_x}_{\sim \varepsilon^2} - \underbrace{u_y}_{\sim \varepsilon \times \varepsilon} \underbrace{B_z}_{\sim O(1)} \right] - \frac{1}{mn} \underbrace{\frac{\partial \Pi_{xx}}{\partial x}}_{\sim \varepsilon^2}, \qquad (27) $$

$$ \underbrace{\frac{\partial u_y}{\partial t}}_{\sim \varepsilon} + \underbrace{u_x \frac{\partial u_y}{\partial x}}_{\sim \varepsilon^2 \times \varepsilon} = -\frac{e}{m} \underbrace{E_y}_{\sim \varepsilon} - \frac{1}{mn} \underbrace{\frac{\partial \Pi_{xy}}{\partial x}}_{\sim \varepsilon}, \qquad (28) $$

These conditions for the purely e.m. Weibel mode are generally valid also in a Vlasov-Maxwell kinetic description, since they follow from the first three velocity moments of the Vlasov equation. The same is true (as long as the heat flux divergence can be neglected in the modeling of Weibel instability, cf. Ref. 39), for the ordering on $\Pi_{yy}$ that can be deduced from its evolution equation,

$$ \underbrace{\frac{\partial \Pi_{yy}}{\partial t}}_{\sim \varepsilon^2} + \underbrace{\Pi_{yy} \frac{\partial u_x}{\partial x}}_{\sim \varepsilon^2 \times \varepsilon^2} + 2 \underbrace{\Pi_{xy} \frac{\partial u_y}{\partial x}}_{\sim \varepsilon \times \varepsilon} = -2 \underbrace{\frac{e}{mc} \Pi_{xy} B_z}_{\sim \varepsilon \times \varepsilon}. \qquad (29) $$

The ordering $\Pi_{xy} \sim \varepsilon$ 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
<figure><img src="image_12.png" /><figcaption>FIG. 2. Time evolution of the electrostatic energy $\epsilon_{e,x}$ (blue), of the transverse (inductive) electric energy $\epsilon_{e,y}$ (red) and of the magnetic energy $\epsilon_{m,z}$ (yellow) in logarithmic scale. The linear growth of the Weibel instability (curves of $\epsilon_{e,y}$ and $\epsilon_{m,z}$) starts around $t\omega_p \sim 20$, just after a transient regime where the longitudinal electric field is also excited out of the numerical noise related to the machine rounding error. The growth rate of the electrostatic energy grows twice as fast, i.e., as the mode $2k_M$. Simulation case (S1).</figcaption></figure>

result is confirmed by full kinetic simulations (cf. Sec. VA1). 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¹⁵) or by linear mode oscillations (e.g., in the propagation of magneto-acoustic waves⁵⁵), 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_{xx} \sim \delta\Pi_{yy} \sim (\delta\Pi_{xy})^2$ for the Weibel instability, with the fact that only the trace of the pressure tensor contributes to the internal energy $\mathcal{E}_{T,x} + \mathcal{E}_{T,y}$ of the plasma: in diagonalizing the $\Pi_{ij}$ matrix for $i, j = x, y$ so as to determine the local contributions to the internal energy, the off-diagonal corrections $\Pi_{xy} = \delta\Pi_{xy} \sim \epsilon$ contribute like the diagonal variations of order $\epsilon^2$ since $\Pi_{xx} = P_{x,0} + \delta\Pi_{xx} \sim 1 + \epsilon^2$ and $\Pi_{yy} = P_{y,0} + \delta\Pi_{yy} \sim 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).

We can summarize the result of this analysis with the orderings,

$$ \begin{aligned} n &\sim O(1), & u_y &\sim B_z \sim E_y \sim \delta\Pi_{xy} \sim \epsilon, \\ u_x &\sim \delta n \sim E_x \sim \delta\Pi_{xx} \sim \delta\Pi_{yy} \sim \epsilon^2. \end{aligned} \quad (30) $$

With the exception of the pressure tensor components, the orderings are the same as Eqs. (13) and (14).

As a consequence, the conservation of the volume integral of the total energy [(11) and (12)] during the evolution of a Weibel instability in a Vlasov plasma implies the following orderings among the energy contributions:

$$ \begin{aligned} & \int dx \left[ \underbrace{\delta K_y}_{\sim \epsilon^2} + \underbrace{\delta \mathcal{E}_{T,x}}_{\sim \epsilon^2} + \underbrace{\delta \mathcal{E}_{T,y}}_{\sim \epsilon^2} + \underbrace{\delta \epsilon_{e,y}}_{\sim \epsilon^2} + \underbrace{\delta \epsilon_{m,z}}_{\sim \epsilon^2} \right] \\ &= - \int dx \left[ \underbrace{\delta K_x}_{\sim \epsilon^4} + \underbrace{\delta \epsilon_{e,x}}_{\sim \epsilon^4} \right]. \end{aligned} \quad (31) $$

We have used the notation introduced in Eq. (12), by separating the contributions of the three components of the fluid velocity in the kinetic energy $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 \mathcal{E}_{T,z} = 0$).

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

It can be of interest to look at the analysis above from the perspective of generalized Ohm’s law. The model of Eqs. (15)–(17), combined with Maxwell equations, can be seen as an extension of the (incompressible) collisionless EMHD description,⁵⁴ the same way non-polytropic closures are introduced in the so-called extended magnetohydrodynamic models. The evolution of $B$ in EMHD is as usual obtained by combining Faraday equation with the curl of (16)

$$ \frac{\partial}{\partial t} \{ \mathbf{B} - \frac{m}{ec} \mathbf{\nabla} \times \mathbf{u} \} = \mathbf{\nabla} \times \{ \mathbf{u} \times [ \mathbf{B} - \frac{m}{ec} \mathbf{\nabla} \times \mathbf{u} ] \} + \frac{c}{e} \mathbf{\nabla} \times \left( \frac{\mathbf{\nabla} \cdot \mathbf{\Pi}}{n} \right), \quad (32) $$

where

$$ \mathbf{u} = - \frac{c}{4\pi ne} \left( \mathbf{\nabla} \times \mathbf{B} - \frac{1}{c^2} \frac{\partial E}{\partial t} \right), \quad (33) $$

as stated by Ampère’s law combined with $J = -en_e \mathbf{u}$ (immobile ions). If the last pressure-related term at the right-hand side is zero, Eq. (32) states the Lagrangian invariance of $\{\mathbf{B} - (m/\epsilon c)\mathbf{\nabla} \times \mathbf{u}\}$ with respect to the electron fluid $\mathbf{u}$, and from this Alfvén theorem, the frozen-in law and the conservation of magnetic helicity all follow (see, e.g., Ref. 58).

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 $\partial \mathbf{B} / \partial t \neq 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⁶¹ (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 $\nabla n_0 \times \nabla \cdot \mathbf{\Pi}_0 \neq 0$.

Let us now specialize Eq. (32) for the CFI or Weibel instabilities, for which $\partial/\partial t \to \gamma_B$ and $\nabla \to -i(k, 0, 0)$. We first note that

$$ \frac{m}{ec} \nabla \times \mathbf{u} = \frac{n_0 d_e^2}{n} \left[ \frac{1}{c^2} \frac{\partial^2}{\partial t^2} - \nabla^2 \right] \mathbf{B} - \frac{mc}{e} \left( \frac{\nabla n}{n} \times \frac{\mathbf{J}}{en} \right), \quad (34) $$

meaning that in Eq. (34) the contribution of the displacement current (i.e., the $\sim (1/c^2)\partial^2/\partial t^2$ term) is negligible as long as $\gamma_B^2/\omega_{pe}^2 \ll (kd_e)^2$, which is typically the case also for the linear stage of Weibel-type modes. Writing each scalar component as $f = f_0 + f_1$ and using $\delta\mathbf{B} = B_{z,1}e_z$, in the limit in which the displacement current effects are neglected, we obtain

$$ \frac{\partial}{\partial t} \left( B_{z,1} - \frac{m}{ec} \frac{\partial u_{y,1}}{\partial x} \right) = \frac{1}{n_0} \frac{\partial^2 \Pi_{xy,1}}{\partial x^2}, \quad (35) $$

that is,

$$ B_{z,1} = -i \frac{m}{ec} k u_{y,1} - \frac{k^2 \Pi_{xy,1}}{\gamma_B n_0}. \quad (36) $$

Once again, we recover the ordering of the first line of Eq. (30), $E_{y,1}$ being trivially ordered like $B_{z,1}$ via Faraday's equation.

## 2. Fluid vs single-particle dynamics

The interpretation of Weibel instability provided in terms of the generalized Ohm's law can be compared with that provided by the single particle dynamics of Sec. III. As already pointed out, electrons therein act as test particles, which do not affect the imposed, exponentially growing magnetic field. Their dynamics, however, although not self-consistently described, account for the induced polarization of the medium as far as charge separation effects are concerned. From the point of view of single particle dynamics, instead, each electron moving in a non-rectilinear trajectory generates a magnetic field as described by Lenz law, that is, by the surface integral of Faraday's equation. Therefore, the single particle motion occurs by preserving the magnetic flux through the plane on which the electron is moving since there is no electromotive force (consider, e.g., the motion in vacuum). From a formal point of view, this is seen by integrating Faraday's equation over a closed surface evolving in time, whose contour is the trajectory of the electron itself. The deformation in time of the surface element $\Omega$ whose contour is dragged by a velocity field $\tilde{v}$ is then described by the integral $\oint \mathbf{B} \cdot (\tilde{v} \times dl)$, which plays the role of the convective term in the material derivative of the magnetic flux defined as $d\Phi_{\Omega}(\mathbf{B})/dt_{\tilde{v}} = \partial\Phi_{\Omega}(\mathbf{B})/\partial t + \tilde{v} \cdot \nabla\Phi_{\Omega}(\mathbf{B})$, and which Faraday's equation implies to satisfy

$$ \underbrace{\int_{\Omega} \frac{\partial \mathbf{B}}{\partial t} \cdot d\mathbf{S} + \oint_{\partial \Omega} \mathbf{B} \cdot (\tilde{v} \times dl)}_{= \frac{d\Phi_{\Omega}(\mathbf{B})}{dt_{\tilde{v}}}} = \oint_{\partial \Omega} [\mathbf{B} \cdot (\tilde{v} \times dl) - c\mathbf{E} \cdot d\mathbf{l}]. \quad (37) $$

Therefore, the conservation of the magnetic flux across a surface whose border is delimited by the trajectory of a bunch of cold electrons moving with equal velocity $\tilde{v}$ implies that an electrostatic field is generated, which compensates the Lorentz force, i.e., $\delta E = -(\tilde{v}/c) \times \delta \mathbf{B}$ on the right-hand side of Eq. (37). In this sense, the electrostatic field $\delta E$ appears as a second-order effect in Weibel instability so as to allow the

bunch of electrons move in the self-consistently generated magnetic field. In particular, this electric field has components $\delta E_x = -\tilde{v}_y \delta B_z/c$ and $\delta E_y = \tilde{v}_x \delta B_z/c$, whose orderings are consistent with those of Eqs. (13) and (14). The actual motion of an electron in the plasma will be, however, determined also by the contribution to the magnetic field fluctuation coming from the other particles. This is provided by the fluid description, which accounts for the polarization effect on the plasma due to the collective particle response to the e.m. perturbations.

A final conclusion can be drawn from the comparison of the magnetic flux conservation valid for the bunch of cold, test-particle electrons moving with velocity $\tilde{v}$ [see Eq. (37)] with the evolution of the magnetic flux that can be deduced by surface integration of Ohm's law of Eq. (32),

$$ \frac{d\Phi_{\Omega}(\mathbf{B})}{dt_u} = \frac{m}{ce} \frac{d}{dt_u} \left[ \oint_{\partial\Omega} \mathbf{u} \cdot \mathbf{l} \right] + \frac{c}{e} \oint_{\partial\Omega} \frac{\nabla \cdot \mathbf{\Pi}}{n} \cdot d\mathbf{l}. \quad (38) $$

While for each isolated electron (or a bunch of 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 $\Omega$ is chosen to cover the entire subdomain where the sign of $B_z$ is fixed (positive or negative).

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 $\tilde{v}_y > \tilde{v}_x$, provided by the instability condition $P_{y,0} > P_{x,0}$. This allows a reinterpretation 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.

We consider the response of the nonlinear system to an initial magnetic perturbation $\delta B_z(x, 0) = B_M \epsilon(0) \cos(kx)$ introduced in a Weibel-unstable equilibrium (as we will later do in nonlinear kinetic simulations, cf. Sec. V). To this end, let us perform a finite difference integration for small time steps in the interval $[t_0, t_0 + \delta t]$, e.g.,
$\partial F(x, t)/\partial t = G(x, t_0) \Rightarrow F(x, t_0 + \delta t) \simeq F(x, t_0) + G(x, t_0)\delta t$. At the first step, $\delta t$ a solenoidal electric component $\delta E_y(x, \delta t) \sim [\epsilon(0)]\delta t$ follows from Ampère's equation (9) and $\delta\Pi_{xy}(x, \delta t) \sim [\epsilon(0)]\delta t$ follows from Eq. (24). At the next step, these perturbations generate $\delta u_y(x, 2\delta t) \sim [\epsilon(0)]\delta t^2$ [cf. Eq. (21)] and $\delta\Pi_{xx}(x, 2\delta t) \propto [\delta\Pi_{xy}(x, \delta t)\delta B_z(x, \delta t)]\delta t \simeq [\delta\Pi_{xy}(x, \delta t)\delta B_z(x, 0)]\delta t \sim [\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) \sim \partial_x[\delta E_y(x, \delta t)]\delta t \sim [\epsilon(0)]\delta t^2$ [Eq. (8)], we obtain

$$ \delta u_x(x, 3\delta t) = \left[ \frac{e\delta u_y}{mc} \delta B_z - \frac{1}{mn} \frac{\partial \delta \Pi_{xx}}{\partial x} \right] \Big|_{2\delta t} \delta t \\ \sim [\epsilon(0)]^2 \delta t^5 + [\epsilon(0)]^2 \delta t^3 \sim [\epsilon(0)]^2 \delta t^3, \quad (39) $$

which is of course consistent with the orderings (30). Due to the immobile ion hypothesis, this corresponds to a current density $\delta J_x|_{3\delta t} = -e(nu_x)|_{3\delta t} \simeq -en_0\delta u_x|_{3\delta t}$. The x-component of Ampère's law requires then that an electrostatic component $\delta E_x|_{3\delta t} = -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 $\delta 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) \propto [\delta u_y\delta B_z]|_0\delta t \sim [\epsilon(0)]^2\delta t$.

## B. Multistream models for the Weibel instability

The multistream model is a kinetic description that uses the translational invariance of some components of the electron canonical momentum,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 $\mathbf{P}_{c\perp} = \mathbf{p}_{\perp} - (e/c)\mathbf{A}_{\perp}(x, t)$ the invariant components of the electron canonical momentum perpendicular to the wave-vector of the perturbation (in the Coulomb gauge $\mathbf{V} \cdot \mathbf{A}_{\perp} = 0$). The multistream model follows the evolution of each beam of particles having the same value of the $\mathbf{P}_{c\perp}$ 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 "Nth" bunch, its Hamiltonian in the non-relativistic limit reads

$$ H_N(x, v, t) = \frac{mv_x^2}{2} + \frac{(C_N + (e/c)\mathbf{A}_{\perp}(x, t))^2}{2m} - e\phi(x, t). \quad (40) $$

Assuming that the "longitudinal" distribution function of each beam, $f_N(x, v_x, t)$, evolves according to a separate Vlasov equation in the form

$$ \frac{\partial f_N}{\partial t} + [H_N, f_N] = 0, \quad (41) $$

we write

$$ \frac{\partial f_N}{\partial t} + v_x \frac{\partial f_N}{\partial x} - \left( \frac{e}{m} E_x + \frac{\partial}{\partial x} \frac{(C_N + (e/c)\mathbf{A}_{\perp})^2}{2m^2} \right) \frac{\partial f_N}{\partial v_x} = 0. \quad (42) $$

The Gauss and Ampère's equations, which provide the only coupling between the beams, in this case take the form

$$ \frac{\partial E_x}{\partial x} = -4\pi e \left( n_0 - \sum_N n_N(x, t) \right), \quad (43) $$

$$ \frac{\partial^2 \mathbf{A}_{\perp}}{\partial t^2} - c^2 \frac{\partial^2 \mathbf{A}_{\perp}}{\partial x^2} = -4\pi c \sum_N J_{\perp,N}(x, t), \quad (44) $$

where

$$ \begin{gathered} n_N(x, t) \equiv \int f_n(x, t) \, d^3\mathbf{v}, \\ J_{\perp,N}(x, t) \equiv -e \int \mathbf{v}_{\perp} f_n(x, \mathbf{v}, t) \, d^3\mathbf{v}. \end{gathered} \quad (45) $$

The "complete" distribution of each beam reads

$$ F_N(x, v, t) \equiv f_N(x, v_x, t) \delta \left( v_{\perp} - \left[ \frac{C_N}{m} + \frac{eA_{\perp}(x, t)}{mc} \right] \right). \quad (46) $$

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_{\perp} = 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 $\partial(\sum_N f_N)/\partial t + [\sum_N H_N, \sum_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

A further possible simplification consists of assuming each stream to be "cold" with respect to the $v_x$ velocity component. This means we take

$$ f_N(x, v_x, t) = n_N(x, t) \delta(v_x - u_{x,N}(x, t)) \quad (47) $$

at any time during the plasma evolution. This automatically fixes a closure on the second-order moment of the reduced Vlasov equation.
Equation (42) combined with (47) is indeed equivalent to the N fluid equations obtained from the zeroth and first-order moment of (42),

$$ \frac{\partial n_N}{\partial t} + \frac{\partial}{\partial x}(n_N u_{x,N}) = 0, \qquad (48) $$

$$ \frac{\partial u_{x,N}}{\partial t} + u_{x,N} \frac{\partial u_{x,N}}{\partial x} = - \frac{e}{m} E_x - \frac{\partial}{\partial x} \left( \frac{C_{y,N} + (e/c)A_y}{2m} \right)^2. \qquad (49) $$

Using

$$ C_{y,N} = m u_{y,N}(x, t) - \frac{e}{c} A_y(x, t) \qquad (50) $$

in principle allows one to fix $C_{y,N}$ from $u_{y,N}(x, 0) = u_{y,N}^{(0)}$ and from $A_y(x, 0)$ that gives the initial perturbation on $B_z$. In the case of the Weibel instability, one typically neglects $eA_y(x, 0)/c$ so that $u_{y,N}(x, 0) = u_{y,N}^{(0)} = C_{y,N}$.

Regardless of the value of $u_{y,N}(x, 0) = u_{y,N}^{(0)}$,

$$ \Pi_{xx,N}(x, t) = \Pi_{xy,N}(x, t) = \Pi_{yy,N}(x, t) = 0, \quad \forall N, \forall t, \qquad (51) $$

where we have defined

$$ \Pi_{ij,N}(x, t) \equiv mn_N \langle (v_i v_j)_N - \langle v_i \rangle_N \langle v_j \rangle_N \rangle \qquad (52) $$

$$ = \left( m \iint v_i v_j F_N dv_x dv_y \right) - mn_N u_{i,N} u_{j,N} \qquad (53) $$

by introducing the ensemble average $\langle \dots \rangle_N$ over the distribution labeled with N.

To make comparison with Weibel instability, we consider the global pressure tensor evaluated for the superposition of all streams in the reference frame of the whole electron distribution,

$$ \Pi_{ij}^{tot}(x, t) \equiv m \iint \sum_{N} (v_i - \bar{u}_i)(v_j - \bar{u}_j) F_N dv_x dv_y, \qquad (54) $$

where

$$ \bar{u}_i \equiv \frac{1}{\sum_N n_N} \iint \sum_N v_i F_N dv_x dv_y = \frac{\sum_N n_N u_{i,N}}{\sum_N n_N}. \qquad (55) $$

Using $\langle (v_i - u_i)(v_j - u_j) \rangle = \langle v_i v_j \rangle - u_i u_j$, where $\langle \dots \rangle$ stands for the integral over the velocities $dv_x dv_y$, weighted by the relevant distribution function, we write

$$ \Pi_{ij}^{tot} = \sum_{N} \Pi_{ij,N} + m \sum_{N} n_N u_{i,N} u_{j,N} - m \sum_{N} n_N \bar{u}_i \bar{u}_j. \qquad (56) $$

In spite of Eq. (51), the total pressure components can be non-zero due to the last two contributions related to the fluid velocity. In the case of Weibel instability with initially cold beams, one obtains indeed, using (51) and the last of Eq. (55),

$$ \Pi_{ij}^{tot} = m \left[ \sum_{N} n_N u_{i,N} u_{j,N} - \frac{\left( \sum_{N} n_N u_{i,N} \right) \left( \sum_{N} n_N u_{j,N} \right)}{\sum_{N} n_N} \right]. \qquad (57) $$

Assuming the validity of (41), the case of cold streams with $u_{x,N}(x, 0) = u_{x,N}^{(0)} = 0 \forall N$ implies $P_{y,0} > 0$ and $P_{x,0} = 0$, identifiable as the $P_{y,0} \gg P_{x,0}$ limit of Weibel instability. A correspondence between the parameters of the kinetic distribution (1) and of the beams of Eq. (47) can be established thanks to (57). If we pose $A_y(x, 0) = 0$ in Eq. (50), we find

$$ n_0 = \sum_N n_N, \qquad (58) $$

$$ P_{y,0} = \frac{1}{m} \left[ \sum_N n_N C_{y,N}^2 - \frac{1}{n_0} \left( \sum_N n_N C_{y,N} \right)^2 \right]. \qquad (59) $$

It was shown³⁷ 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.⁶⁴,⁶⁵ For Weibel-type instabilities, this approach has been first used for relativistic distributions⁶⁶ and then for non-relativistic ones, too.⁴⁹,⁶³ Operationally speaking, we express $f_N(x, v_x, t)$ in terms of Heaviside step functions having value $\mathcal{F}_{f,N} > 0$ for $v_x \in [-w_N^{\pm}(x, t), +w_N^{\pm}(x, t)]$ and zero otherwise, $w_N^{\pm}(x, 0)$ being assigned at the initial time $t = 0$ for each waterbag $l$.

For the case of a single waterbag, we write

$$ \left\{ \begin{array}{l} f_N(x, v_x, t) = n_0 A_N [\mathcal{H}(v_x - w_N^{-}) - \mathcal{H}(v_x - w_N^{\pm})], \\ n_N(x, t) \equiv \int_{-\infty}^{+\infty} f_N(x, v_x, t) dv_x = n_0 A_N (w_N^{\pm} - w_N^{-}), \end{array} \right. \qquad (60) $$

where $A_N$ is the “amplitude” of the bag, in units of (velocity)$^{-1}$. Owing to Liouville’s theorem, the dynamics of the bag is completely determined by its contours,⁶⁴,⁶⁵ that is, by

$$ \frac{\partial w_N^{\pm}}{\partial t} + w_N^{\pm} \frac{\partial w_N^{\pm}}{\partial x} = - \frac{e}{m} E_x - \frac{1}{2m} \frac{\partial}{\partial x} (C_{y,N} + eA_y)^2. \qquad (61) $$

The evolution of the fluid moments of Eq. (60) is thus determined. Note that $w_N^{\pm}(x, t)$ and $w_N^{\pm}(x, t)$ separately evolve according to Eq. (61), whereas they appear combined in the fluid moments’ evolution. As first shown in Ref. 65 (cf. also Ref. 67, Secs. 2.4–2.5), Eq. (61) corresponds indeed to the combination of the two lowest order moment equations, closed with a polytropic-like condition with index 3 [cf. Eq. (25)],

$$ \frac{\partial n_N}{\partial t} + \frac{\partial}{\partial x} (n_N w_N) = 0, \quad w_N \equiv \frac{w_N^{\pm} + w_N^{-}}{2}, \qquad (62) $$

$$ \frac{\partial w_N}{\partial t} + w_N \frac{\partial w_N}{\partial x} = - \frac{1}{mn_N} \frac{\partial \Pi_{xx,N}}{\partial x} - \frac{e}{m} E_x - \frac{1}{2m} \frac{\partial}{\partial x} (C_{y,N} + eA_y)^2, \qquad (63) $$
$$ \Pi_{xx,N}^{(l)} = \frac{mn_N}{12} (w_N^+ - w_N^-)^2 = \frac{mn_{I,N}^3}{12n_0^2 A_{I,N}^2}. \quad (64) $$

Here, $w_N$ 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_{xx,N}$ component in this modeling. Equation (64) can be directly computed from the first-order moment of Vlasov equation, by identifying $mn_N^3/(12n_0^2 A_{I,N}^2)$ as $\Pi_{xx,N}$ after direct comparison with the fluid equations. The argument can be generalized to the case of many bags, say for $l = 1, \dots, M$, where Eq. (60) is substituted by

$$ \left\{ \begin{aligned} f_N &= \sum_{l=1}^{M} n_0 A_{l,N} \left[ \mathcal{H}(v_x - w_N^{l-}) - \mathcal{H}(v_x - w_N^{l+}) \right], \\ n_{l,N}(x, t) &= n_0 A_{l,N} (w_N^{l+} - w_N^{l-}), \\ n_N(x, t) &= \sum_{l=1}^{M} n_{l,N} = n_0 \sum_{l=1}^{M} A_{l,N} (w_N^{l+} - w_N^{l-}). \end{aligned} \right. \quad (65) $$

In this case, defining

$$ w_{l,N} \equiv \frac{w_N^{l,+} + w_N^{l,-}}{2}, \quad (66) $$

Eqs. (61)–(64) correspond to

$$ \begin{aligned} & \frac{\partial}{\partial t} \sum_{l}^{M} n_{l,N} + \frac{\partial}{\partial x} \left( \sum_{l}^{M} n_{l,N} w_{l,N} \right) = 0, && (67) \\ & \frac{\partial}{\partial t} \sum_{l}^{M} n_{l,N} w_{l,N} + \frac{\partial}{\partial x} \sum_{l}^{M} n_{l,N} w_{l,N}^2 \\ & = -\frac{\partial}{\partial x} \sum_{l}^{M} \frac{\Pi_{xx,N}^{(l)}}{m} - \sum_{l}^{M} n_l \left( \frac{e}{m} E_x + \frac{1}{2m} \frac{\partial}{\partial x} (C_{y,N} + eA_y)^2 \right), && (68) \\ & \Pi_{xx,N}^{(l)} \equiv mn_{l,N} (\langle v_x v_x \rangle_{l,N} - \langle v_x \rangle_{l,N}^2) \\ & = \frac{mn_{l,N}}{12} (w_N^{l,+} - w_N^{l,-})^2 = \frac{mn_{l,N}^3}{12n_0^2 A_{l,N}^2}. && (69) \end{aligned} $$

$$ \Pi_{xx,N}^{(l)} \equiv mn_{l,N} (\langle v_x v_x \rangle_{l,N} - \langle v_x \rangle_{l,N}^2) = \frac{mn_{l,N}}{12} (w_N^{l,+} - w_N^{l,-})^2 = \frac{mn_{l,N}^3}{12n_0^2 A_{l,N}^2}. \quad (69) $$

It is interesting to combine the results of the waterbag formalism with the multi-stream model. We use Eqs. (54)–(56) to evaluate the components of the pressure tensor, which, in the case of “cold streams,” were initially null. Then, the fluid velocity component $u_{x,N}$ of the whole water-bag distribution for the Nth stream reads

$$ u_{x,N} = \frac{\sum_{l=1}^{M} n_{l,N} w_{l,N}}{\sum_{l}^{M} n_{l,N}} = \frac{\sum_{l}^{M} A_{l,N} \left[ (w_N^{l,+})^2 - (w_N^{l,-})^2 \right]}{2 \sum_{l}^{M} A_{l,N} (w_N^{l,+} - w_N^{l,-})} . \quad (70) $$

From the definition of $\Pi_{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

$$ \begin{aligned} \Pi_{xx,N} &\equiv mn_N (\langle v_x v_x \rangle_N - u_{x,N} u_{x,N}) \\ &= \sum_{l=1}^{M} \left( 4\Pi_{xx,N}^{(l)} + mn_{l,N} w_N^{l,+} w_N^{l,-} \right) - mn_N u_{x,N}^2. \end{aligned} \quad (71) $$

This is the only difference with respect to the cold stream limit of Eq. (51): an analogous procedure to evaluate the other pressure components still yields $\Pi_{xy,N}(x, t) = \Pi_{yy,N}(x, t) = 0\sqrt{N}$. These conditions can be combined with Eq. (56) to evaluate the total pressure components $\Pi_{ij}^{(tot)}$ in the multi-stream/multi-waterbag model.

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

We recognize a correspondence between the components $\Pi_{ij}^{(tot)}$ obtained by integrating (i.e., averaging) over the $N = 1, \dots, N$ streams and the quantities $\{n_N, C_N; f_N(x, p_x, t)\}$ associated with each of them,

$$ \begin{aligned} \Pi_{xx}^{(tot)} &\leftrightarrow \{n_N, C_{y,N}; f_N\}, \\ \Pi_{yy}^{(tot)} &\leftrightarrow \{n_N, C_{y,N}^2; f_N\}, \\ \Pi_{xy}^{(tot)} &\leftrightarrow \{n_N, C_{y,N}; f_N\}. \end{aligned} $$

$\Pi_{yy}^{(tot)}$ displays a quadratic dependence on the coefficients $C_{y,N}$, making it insensitive of their sign, whereas $\Pi_{xx}^{(tot)}$ and $\Pi_{xy}^{(tot)}$ do depend on the sign of each $C_{y,N}$. The initial distribution defining $\Pi_{xx}^{(tot)}|_{t=0} = P_{x,0}$ and $\Pi_{yy}^{(tot)}|_{t=0} = P_{y,0}$ can be described using a sufficiently small number of “streams” by appropriately choosing the parameters $\{n_N, C_{y,N}\}$. The total number of streams can be conveniently chosen to be $2N + 1$, i.e., in pairs having non-zero velocities plus one stream with zero velocity (e.g., $N = -N, \dots, N$). Given whatsoever initial distribution function $f(x, p_x, p_y, t)$, it is then possible to determine the values of the parameters of the streams so that an equivalence in terms of the moments is established between the full kinetic and the reduced kinetic model. Since the linear growth of the Weibel instability depends just on the first three moments of $f$ (namely, initial density, average velocity, and “temperature”), this allows for an “exact” correspondence between the two linear descriptions of the instability.

Naming $u_{x,N}^{(0)} = u_{x,N}(x, 0)$, the elements of the full kinetic dispersion matrix [(A1) and (A2)] can be thus recovered,

$$ \varepsilon_{xx} = 1 - \frac{\omega_p^2}{n_0} \sum_N \frac{n_N}{\omega^2 - (ku_{x,N}^{(0)})^2}, \quad (72) $$

$$ \varepsilon_{yy} = \omega^2 - k^2 c^2 - \omega_p^2 + \sum_N \frac{n_N C_{y,N}^2}{n_0 m^2} \left( \frac{1}{d_e^2} - \frac{k^2 \omega_p^2}{\omega^2 - (ku_{x,N}^{(0)})^2} \right), \quad (73) $$

$$ \varepsilon_{xy} = -ik \frac{\omega_p^2}{n_0} \sum_N \frac{n_N (C_{y,N}/m)}{\omega^2 - (ku_{x,N}^{(0)})^2}. \quad (74) $$

The resonant denominators are related to the presence of a longitudinal temperature (cf. Refs. 40 and 68), thanks to the waterbag description. When the condition $\sum_N n_N C_{y,N} = 0$ is satisfied, e.g., when $2N + 1$ symmetric streams are taken into account for describing the distribution function in the perpendicular direction, then $\varepsilon_{xy} = 0$. This means that the longitudinal and transverse modes are decoupled and we recover the dispersion relation of the Weibel instability from $\varepsilon_{yy} = 0$,

$$ \omega^2 = k^2 c^2 + \omega_p^2 - \sum_{N=-N}^{+N} \frac{n_N C_{y,N}^2}{n_0 m^2} \left( \frac{1}{d_e^2} - \frac{k^2 \omega_p^2}{\omega^2 - (ku_{x,N}^{(0)})^2} \right). \quad (75) $$
The condition $\epsilon_{xx} = 0$ gives instead the dispersion relation of the longitudinal, electrostatic modes: when $u_{x,N}^{(0)} = 0$, the dispersion relation of Langmuir waves, $\omega = \pm \omega_p$, is recovered. The thermal velocity $(v_x^{th})^2$ in the waterbag approach depends on the velocities $(u_{x,N}^{(0)})^2$ [cf. Eq. (71)], whereas, in the multistream approach, $(v_y^{th})^2 = \sum_N (n_N/n_o)(C_{y,N}/m)^2$. The limit $P_{y,0} \gg P_{x,0}$ can be thus made correspond to $\Pi_{yy,N}^{(0)} \gg \Pi_{xx,N}^{(0)}$. This is compatible with the limit $u_{x,N}^{(0)} \to 0$, which only requires conditions on $n_N$ and $C_{y,N}^2$. Then, compatibly with the non-relativistic limit $c^2 \gg (\omega/k)^2$, Eq. (75) reduces to Weibel's result,²⁵

$$ \omega^2 \approx k^2 c^2 + \omega_p^2 \left( 1 + \frac{(kv_y^{th})^2}{\omega^2} \right). \qquad (76) $$

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

Each set of particles, belonging to the stream 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_{bN}$, 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),

$$ \omega_{bN} = \omega_p \sqrt{\frac{k_c e B_M C_{y,N}}{\omega_p m \omega_p mc}} = \omega_p \sqrt{kd_e \Omega \frac{C_{y,N}}{mc}}. \qquad (77) $$

The differential rotation of current filaments around the extremum values of the magnetic field, of which we highlighted the implications for the electron density concentration during the linear stage of the instability, is also the key element of the well-known magnetic trapping mechanism that occurs at saturation. Introducing in our discussion the notion of magnetic trapping is useful for the interpretation of the numerical results, which we are going to present in Sec. V, in relation to the nonlinear coupling between electromagnetic and electrostatic perturbations. As we will see, although perturbations on $E_x$ are much weaker in amplitude, when compared to those on $B_z$, 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 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.⁴³ 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, $\lambda_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⁴³ 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} = 20P_{x,0} = 20$ keV. Since for this value of the pressure anisotropy the linear dispersion relation $\epsilon_{yy} = 0$ given by Eq. (A2) predicts a maximal growth rate $\gamma_M = 0.129\omega_{pe}$ 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 \cos(k_M x)$, with $B_0 = 10^{-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 $\sim 10^{-5}$ amplitude of the growing Weibel mode. The other numerical parameters are $N_x = 512$, $N_{v_x} = N_{v_y} = 96$, and $\Delta t = 0.005 \omega_p^{-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).

(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} = 25P_{x,0} = 25$ keV, which correspond to a fastest Weibel mode with $\gamma_M = 0.142\omega_{pe}$ 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 \sin(k_M x)$, with $B_0 = 10^{-5}$. The other numerical parameters are $N_x = 256$, $N_{v_x} = N_{v_y} = 256$, and $\Delta t = 0.005 \omega_p^{-1}$. The difference with respect to the simulation case (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^{-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$.
<figure><img src="image_13.png" /><figcaption>FIG. 3. Spatial profile of the electron density $n(x, t)$ (blue) and of the magnetic field $B_z(x, t)$ (orange) at time $t\omega_{pe} = 0.5$, i.e., at an early time after initialization (left frame), and at a later time (right frame). The blue curve in the left frame is affected by numerical fluctuations related to the small value of the density perturbation at this early time of the instability, in which the dominant Weibel eigenmode is still in its early growth stage (cf. Fig. 2). Simulation case (S1).</figcaption></figure>

## 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 $\epsilon_{m,z} = 1/(8\pi L_x) \int_0^{L_x} B_z^2 dx$, the (inductive) electric energy component $\epsilon_{e,y} = 1/(8\pi L_x) \int_0^{L_x} E_y^2 dx$, and the electrostatic energy $\epsilon_{e,x} = 1/(8\pi L_x) \int_0^{L_x} E_x^2 dx$. Both the energy components related to $B_z^2$ and $E_y^2$ grow with a rate $\sim 0.26 \approx 2\gamma_M$, in excellent agreement with the linear theory, whereas $\epsilon_{e,x}$ grows with a rate, which is about $\approx 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 \sim 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 $\sim 10 t\omega_{pe}$, if we look at the energy components related to the Weibel mode, but it is even shorter ($\sim 5 t\omega_{pe}$) if we measure it with respect to the electrostatic component $\epsilon_{e,x} = \int E_x^2 dx$. 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,²⁷⁻²⁹,³¹,³³⁻³⁵ 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_{em} = \sqrt{\omega_p^2 + (k_M c)^2} \approx \omega_p$. All energy components indicate that saturation occurs around $t\omega_p \approx 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_{pe} = 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

<figure><img src="image_14.png" /><figcaption>FIG. 4. Time evolution of the fluctuations of the diagonal pressure tensor components, i.e., $\delta\Pi_{xx}(x, t) = \Pi_{xx}(x, t) - P_{x,0}$ (upper frame) and $\delta\Pi_{yy}(x, t) = \Pi_{yy}(x, t) - P_{y,0}$ (lower frame), evaluated in the position where their maxima are respectively localized (cf. Fig. 5): in the same plot the results are presented both in linear scale (orange curve) and in semi-logarithmic scale (blue curve). The latter allows a direct comparison with the growth rate $\gamma_M = 0.129$ of the magnetic field $B_z$ (black dashed line): both $\delta\Pi_{xx}(x, t)$ and $\delta\Pi_{yy}(x, t)$ grow as the density fluctuations, at twice the growth rate of the dominant Weibel mode (value). Simulation case (S1).</figcaption></figure>

fundamental mode $k_{Mc}/\omega_{pe} = 2.017$, the electron density oscillates with $2k_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_z$ with respect to $x$ is positive or negative, depending on whether the sign of the $v_y$ 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 \approx \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) - u_{x,0} \approx u_{x,2}$, with respect to linearized quantities (not shown here).

Figure 4 displays the time evolution of the fluctuations of the diagonal pressure tensor components, i.e., $\delta\Pi_{xx}(x, t) = \Pi_{xx}(x, t) - P_{x,0}$ and $\delta\Pi_{yy}(x, t) = \Pi_{yy}(x, t) - P_{y,0}$, evaluated in the position where their maxima are respectively localized (cf. Fig. 5). The components $\Pi_{ij}(x, y)$ are computed from the kinetic simulations by integrating the second-order moment of the distribution function according to $\Pi_{ij}(x, y) = (1/\int fd^3v) \int (v_i - u_i)(v_j - u_j)fd^3v$ and $u_i = (1/\int fd^3v) \int v_i fd^3v$. The logarithmic scale plots in Fig. 4 show that both the $\Pi_{xx,1}$ and $\Pi_{yy,1}$ grow exponentially with a growth rate $2\gamma_M$, which equals that of the electrostatic fluctuations. Coherently, the spatial profiles of $\delta\Pi_{xx}(x, t)$ and $\delta\Pi_{yy}(x, t)$, shown in Fig. 5 at $t\omega_{pe} = 0.5$, display the $\sim \cos(2k_Mx)$ 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_{xy} = \Pi_{xy}(x, t)$ (we recall that at equilibrium $\Pi_{xy} = 0$), instead, exponentially grows in time and is
<figure><img src="image_15.png" /><figcaption>FIG. 5. Spatial profiles of the fluctuations of the diagonal components of the pressure tensor $\delta\Pi_{xx} = \Pi_{xx,1}$ (blue line in left frame) and of $\delta\Pi_{yy} = \Pi_{yy,1}$ (blue line in right frame) at an early time of the simulation. Superposed, the profile of $B_z(x, t)$ (orange line) at the same time. Simulation case (S1).</figcaption></figure>

modulated in space as the linear Weibel mode, i.e., it has growth rate $\gamma_M$, and spatial oscillation corresponding to the wave vector $k_M$ (Fig. 6).

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).

Figure 7 displays in the bottom frame the time evolution of the magnetic energy $\epsilon_{m,z}$ in semi-logarithmic scale. From this, the numerical value of the growth rate of the dominant Weibel mode is estimated to be close to $0.141\omega_p$. This is in very good agreement with the maximal growth rate obtained from the dispersion relation, $\gamma_M \simeq 0.142\omega_p$ (top frame of the same figure). Saturation occurs at time $t\omega_p \simeq 80$ and a transient stage is visible for $t\omega_p \lesssim 20$, with time oscillations of the energy at about twice the frequency of e.m. waves propagating in the plasma, $\omega_{em} = \sqrt{\omega_p^2 + k_0^2 c^2} \simeq 1.72\omega_p$. The abrupt decrease in the electrostatic field energy $\epsilon_{e,x}$ at saturation indicates that the dominant mechanism is magnetic trapping, whose characteristic bounce frequency can be numerically estimated from the temporal oscillations of $\epsilon_{m,z}$ for $t\omega_p > 80$ to be $\omega_b \simeq 0.20\omega_p$. This value is again in good agreement, although somewhat larger than, with the theoretical value estimated from Eq. (77) and from the numerical results,

$$ \frac{\omega_b}{\omega_p} = \sqrt{\frac{k_0 c v_y^{th}}{\omega_p c} \frac{eB_M}{m\omega_p}} = \sqrt{1.4 \times 0.20 \times 0.10} \approx 0.167. \quad (78) $$

<figure><img src="image_16.png" /><figcaption>FIG. 6. Time evolution of the off-diagonal component $\delta\Pi_{xy} = \Pi_{xy}(x, t)$ (we recall that at equilibrium $\Pi_{xy} = 0$). This component evolves in time and is modulated in space as the linear Weibel mode, i.e., with a growth rate similar to that of WI (al left panel) and with a spatial oscillation corresponding to the fundamental wave vector (at right). This result confirms the ordering $\delta\Pi_{xy} \sim \delta B_z \sim \epsilon$ of Sec. IVA simulation case (S1).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 7. Top frame: growth rate $\gamma/\omega_p$ of the Weibel modes as a function of $kd_e$, numerically evaluated from Eq. (A2) for $P_{y,0}/P_{x,0} = 25$. Bottom frame: Time evolution of the magnetic field energy $\epsilon_{m,z}$ in a semi-log scale: the slope of the magnetic energy increase is $\sim 0.282\omega_p$, about twice the peak value of growth rate shown in top frame ($\gamma_M \simeq 0.142\omega_p$) (cf. with left frames of Fig. 8, for a linear time representation). Simulation case (S2).</figcaption></figure>

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 $2k_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 $eE_x/m\omega_p c$ (top) and $eE_y/m\omega_p c$ (middle) and the magnetic field component $eB_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 $\omega_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 $2k_M = 2k_0$ mode. This double period appears in the various components of the electromagnetic field, but it is more evident in the $E_x$ component (top frame of Fig. 8).

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 $eE_x/m\omega_p c$ and $eB_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_z$ component (initially excited with amplitude $eB_0/m\omega_p \sim 10^{-5}$): this can be seen also in the spatial
<figure><img src="image_18.png" /><figcaption>**FIG. 8.** Left frames: time evolution of the electromagnetic field components involved in the growth of the Weibel mode, measured at the center of the simulation box as a function of time. The two electric field components $eE_x/m\omega_p c$ (top frame), and $eE_y/m\omega_p c$ (middle frame), and the magnetic field component $eB_z/m\omega_p$ (bottom frame) are shown. Right frames: time evolution of the energy components corresponding to the field components in the left frames. The energy associated with the solenoidal electric component, $\epsilon_{e,x}$ (top frame), to the electrostatic component $\epsilon_{e,y}$ (middle frame) and to the magnetic field $e_{m,z}$ (bottom frame, cf. with bottom frame of Fig. 7 for its semi-log representation) is shown. In the saturation regime (for $t\omega_p \ge 80$), the components oscillate mainly at frequency $\omega_b$. However, the $2k_M$ mode is also present and induces a disturbance with a pulsation close to $2\omega_b$. Simulation case (S2).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**FIG. 9.** Spatial profile, in normalized units, of $E_x(x)$ (left panels), and of $B_z(x)$ (right panels) at different times of the evolution (from top to bottom). Simulation case (S2).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**FIG. 10.** Projection of the electron distribution function on the $(p_x, p_y) = m(v_x, v_y)$ plane at time $t\omega_p = 82.5$ at three different locations in space: at $x = L_x/4$, corresponding to a local maximum positive value of $B_z$ (upper panel); at the center of the box (middle panel), corresponding to a null value of $B_z$; and finally at $x = 3L_x/4$ (bottom panel), corresponding to a maximum negative value of $B_z$, i.e., to a local minimum. The sense of rotation of the distribution function is induced by the local sign of the magnetic field. For a zero magnetic field, the distribution remains perfectly symmetric in the momentum space. Simulation case (S2).</figcaption></figure>

profile, modulated for $E_x$ as the first harmonic of the dominant Weibel mode, visible as the fundamental mode in the spatial oscillation of the $B_z$ 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 $4k$ in the $E_x$ profile indicates the beginning of a wave-breaking process, while the magnetic field is still dominated by the fundamental mode (bottom panels).

Despite the saturation stage seems to be dominated by magnetic trapping, and in spite of the strong decrease in the electrostatic energy component $\epsilon_{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 \ge 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
<figure><img src="image_21.png" /><figcaption>**FIG. 11.** Contour-plot of the distribution function in the $(x, p_y)$ phase-space averaged along $p_x$ and represented at three different times of the saturation stage. In addition to the typical shape of the distribution function (the classical Y shape), there is a modulation on the second harmonic, clearly visible in the upper panel. This can be interpreted as due to the presence of electrostatic fluctuations (cf. Secs. IVA and VA1). Simulation case (S2).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**FIG. 12.** Contour-plot of $f(x, p_x)$ averaged over $p_y$ at different times of the saturation stage: the $2k_M$, dominant at $t\omega_p \simeq 82$ is later coupled to higher order harmonics (e.g., the mode $4k_M$). Simulation case S2.</figcaption></figure>

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 = 3L_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 $\Omega$) 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 $2k$ mode is also present since the beginning of the saturation regime. It is visible in the top panel at $t\omega_p \simeq 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_x$ 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.

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) \int f(x, p_x, p_y) dp_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 $2k$ mode, coherently with the ordering $u_x \simeq u_{x,2} \sim (n(x, t) - n_0) \sim E_x$ previously established [cf. Eqs. (27) and (30)]. The appearance of the $4k$ 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
<figure><img src="image_23.png" /><figcaption>FIG. 13. Contour-plot of $f(x, p_x, p_y, t)$ at different values of $p_y$ and at time $t\omega_p \simeq 82.5$ (cf. Fig. 12). In the center and bottom panels, the notation $p_{y,th} = mv_y^{th}$ has been used. Fixing a value of $p_y$ corresponds to select a population of electrons having the same canonical momentum $p_y \simeq C_{y,N}$, if $|(e/c)A_y| \ll |C_{y,N}|$. Each “population” of electrons experiences its own trapping mechanism which has a hybrid character, ranging from a process dominated by electrostatic fluctuations in which the $2k_M$ mode dominates (top panel) to a process where the e.m. trapping related to $k_M$ mode dominates (bottom panel). The middle panel shows an intermediate case where both magnetic and electrostatic fluctuations seem to have comparable weight. Simulation case (S2).</figcaption></figure>

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 \simeq C_{y,N}$ if $|(e/c)A_y| \ll |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 $2k_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 = mv_y^{th}$ (center panel) or $p_y = 2mv_y^{th}$ (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_M$ 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.

## B. Reduced kinetic simulations via the multi-stream model

We now compare the results of the numerical simulations of Sec. VA 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} = 25P_{x,0} = 25\text{keV}$, and by perturbing the marginal equilibrium with a sinusoidal perturbation $B_z = B_0 \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_{pe}$) and $B_0 = 10^{-5}$. The equivalent anisotropy has been realized by consistently choosing the parameters $n_N$ for a total number of seven streams ($N = -3, -2, \dots, 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\text{keV}$ (i.e., $P_{x,0} = 1\text{keV}$ in normalized units) has been initially chosen for each beam. This means that we have taken $f_N(x, v_x, 0) \propto \exp[-mv_x^2/(2k_B T_{x,0})]$ in Eq. (46).

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_y$ sampled in Fig. 12: $u_y^0 = 0, \pm v_y^{th}, \pm 2v_y^{th}, \pm 3v_y^{th}$. Thanks to the reduced computational cost of the multi-stream model, a finer phase-space mesh has been used with respect to the simulation (S2), namely, $N_x N_v = 256 \times 512$ points. Even in this simulation, the time step has been chosen to be $\Delta t = 0.005 \omega_p^{-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. VA 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_{em} = \sqrt{\omega_p^2 + k_M^2 c^2} \simeq 1.72\omega_p$, the magnetic field $B_z$ grows exponentially with a rate $\gamma_M \simeq 0.140\omega_p$, until nonlinear saturation is reached at $t\omega_p \simeq 80$. Low
<figure><img src="image_24.png" /><figcaption>**FIG. 14.** Illustration of the Fourier component $k_0$ of the electromagnetic field $B_z$, as a function of time (top). The same component is plotted in logarithmic scale at the bottom. The physical characteristics are perfectly found here: after a transient phase, dominated by the appearance of high frequency oscillations, the magnetic field $B_z$ undergoes an exponential growth, with a growth rate close to $\eta_{num} = 0.140\omega_p$, followed by the nonlinear saturation regime. The simulation was performed using the multistream code.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**FIG. 15.** Contour-plot of $f(x, p_x, p_y, t)$ at different values of $p_y$ and at time $t\omega_p \simeq 82.5$. These results, obtained with the multi-stream model, must be compared to those of Fig. 13, which they mirror. Simulation case (S3).</figcaption></figure>

frequency oscillations appear in the saturation regime, at a frequency close to the magnetic bounce frequency, i.e., at $\omega_b \simeq 0.21\omega_p$, in very good agreement with the results of the simulation (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 = mv_y^{th}$, and $p_y = 2mv_y^{th}$ 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_y = 0$, for which the electrostatic mode $2k$ 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 = mv_y^{th}$ (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. VA and

VB and with the fact that in the multi-stream modeling the dynamics of populations with different values of $C_{y,N} \simeq 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 $2k_M$ mode (Fig. 15, top), i.e., by electrostatic (second-order) effects, while the dynamics of the beams positioned at $p_y = p_{y,th}$ and $p_y = 2p_{y,th}$ 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 $2k_M$ of the electrostatic field $E_x$ (top frame) and the Fourier component $k_M$ of both the inductive electric field $E_y$ (middle frame) and $A_y$ vector potential (bottom frame). These results further support the interpretation given in Fig. 15 and agree very well with the results of simulation case S2. The component $|\hat{E}_x(2k_M)|$ reaches its maximum amplitude at a time similar to that of $|\hat{E}_y(k_M)|$, and its growth is about $2\gamma_M$, i.e., twice the growth of the electromagnetic components $|\hat{E}_y(k_M)|$ and $|\hat{A}_y(k_M)|$. At saturation, the $|\hat{E}_y(k_M)|$ component oscillates at about twice the bounce frequency $\omega_b$ 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 $E_y$.
<figure><img src="image_26.png" /><figcaption>FIG. 16. Illustration of the temporal dynamics of the plasma. The main Fourier modes have been plotted: of the electrostatic field $E_x$ (first harmonic $2k_0$ at the top), of the electric component $E_y$ ($k_0$ mode in the middle), and finally of the $A_y$ vector potential component ($k_0$ electromagnetic mode at the bottom), vs time. The simulation was performed using the multistream code for the simulation case (S3).</figcaption></figure>

Finally, in Fig. 17, we show the time evolution of

$$ \langle \tilde{\Pi}_{xx} \rangle_{L_x} \equiv \sum_{N=-3}^{+3} \int_{0}^{L_x} \frac{dx}{L_x} \int_{-\infty}^{+\infty} mv_x^2 f_N(x, v_x, t) dv_x. \quad (79) $$

This quantity is determined from the global distribution function data provided by the multi-stream code. Since the code integrates the relativistic Vlasov-Maxwell system, the second integral in Eq. (79) should be rather read by replacing $mv_x^2 dv_x \to p_x^2 / (\langle \Gamma \rangle m^2) dp_x$ with $\langle \Gamma \rangle$ relativistic Lorentz factor averaged over the beams.49 However, $\langle \Gamma \rangle \simeq 1$ in the non-relativistic limit considered in this simulation case. Using then the definitions (52) and (54), one verifies that at $t=0 \langle \tilde{\Pi}_{xx} \rangle_{L_x}$

<figure><img src="image_27.png" /><figcaption>FIG. 17. Time evolution of the quantity $\langle \tilde{\Pi}_{xx} \rangle_{L_x}$ of Eq. (79) At the beginning of the saturation phase, $\langle \tilde{\Pi}_{xx} \rangle_{L_x} \simeq \int \Pi_{xx}^{tot} dx/L_x$ behaves like the electric field $E_x$, i.e., undergoes an amplification at time $t\omega_p = 60$, with a growth rate $2\gamma_M$. The saturation regime remains dominated by magnetic trapping. The simulation was performed using the multi-stream code.</figcaption></figure>

$\simeq \int \Pi_{xx}^{tot} dx/L_x$ [cf. Eq. (56)]. This provides an excellent approximation of $\int \Pi_{xx}^{tot} dx/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, $\langle \tilde{\Pi}_{xx} \rangle_{L_x}$ behaves like the electric field $E_x$, 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 $\omega_b$, the same of the saturation stage of the beams with initial non-zero velocity.

# VI. SUMMARY AND CONCLUSIONS

## A. Summary

In this work, we have revised the mechanism of linear and non-linear 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 VLEM43 (Sec. VA) and a numerical code implementing the multi-stream model37 (Sec. VB). 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.²⁶ and following authors²⁷,²⁹ 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,¹⁻¹¹ 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 $\epsilon$ 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 \ll \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_{\perp k} > P_{\parallel 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_{xy}$ component) is linear, and thus proportional to the magnetic field amplitude, the consistency of the linear analysis requires $\delta\Pi_{xx}$ and $\delta\Pi_{yy}$ to be of second order with respect to the magnetic field amplitude (we recall that $k = ke_x$). An anisotropic response of thermal fluctuations to harmonic perturbations of coupled e.m./compressible type was already noted for other e.m. modes in collisionless plasmas.⁵⁵

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.²⁶⁻³⁵

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.¹⁸⁻²² 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,¹⁸,⁶⁹⁻⁷¹ where Weibel-type instabilities coupled with other beam-plasma instabilities may be in principle expected to also play a role¹⁶ 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,³⁴,⁷²⁻⁷⁸ 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⁴⁶ 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³⁹ (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⁶⁶,⁷⁹,⁸⁰ and to the its time-resonant counterpart,⁴⁰,⁶⁸,⁸¹⁻⁸³ and finally to relativistic oblique modes, in which a richer interplay between magnetic, electrostatic, and thermal effects has been pointed out.³⁴,⁷⁷,⁸⁴⁻⁸⁸

## 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

Linearization of Eqs. (4) and (8)–(10) leads<sup>25</sup> to the linear system $[\mathbf{e}] \cdot \mathbf{E}_1 = \mathbf{0}$, where the dispersion matrix $[\mathbf{e}]$ has a diagonal form with components

$$ \varepsilon_{xx} = \frac{\omega^2}{c^2} \left[ 1 + \omega_p^2 \frac{1 + \xi Z(\xi)}{k^2 P_{x,0}} \right], \quad (A1) $$

$$ \varepsilon_{yy} = \frac{\omega^2}{c^2} - k^2 - \frac{1}{d_e^2} \left[ 1 - (1 + \xi Z(\xi)) \frac{P_{y,0}}{P_{x,0}} \right], \quad (A2) $$

$$ \varepsilon_{zz} = \frac{\omega^2}{c^2} - k^2 - \frac{1}{d_e^2} \left[ 1 - (1 + \xi Z(\xi)) \frac{P_{z,0}}{P_{x,0}} \right]. \quad (A3) $$

Here $Z(\xi)$ stands for the plasma dispersion function<sup>89</sup> of argument $\xi = \omega/(\sqrt{2}kc_x)$. The condition $\varepsilon_{xx} = 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 $\varepsilon_{yy} = 0$ and $\varepsilon_{zz} = 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 $\varepsilon_{yy} = 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^2c_x^2/\omega^2 \ll 1$ (or equivalently to $\xi \gg 1$), the dispersion relation of Weibel modes (A2) related to a velocity anisotropy in the $(x, y)$ plane becomes

$$ \frac{\omega^2}{c^2} - k^2 - \frac{1}{d_e^2} \left[ 1 + \frac{k^2 c_y^2}{\omega^2} \right] = 0. \quad (A4) $$

The purely e.m. nature of Weibel instability for $\mathbf{k} = k_x \mathbf{e}_x$ and $P_y > P_x$ is evident from the diagonal form of $[\mathbf{e}]$.

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

Fried’s “reduced” model<sup>36</sup> 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,<sup>46</sup> 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)$,

$$ \frac{\partial n_{\alpha}}{\partial t} + \mathbf{\nabla} \cdot (n_{\alpha} \mathbf{u}_{\alpha}) = 0, \quad (B1) $$

$$ \frac{\partial \mathbf{p}_{\alpha}}{\partial t} + (\mathbf{u}_{\alpha} \cdot \mathbf{\nabla}) \mathbf{p}_{\alpha} = -e \left( \mathbf{E} + \frac{\mathbf{u}_{\alpha}}{c} \times \mathbf{B} \right), \quad (B2) $$

$$ \mathbf{p}_{\alpha} = m \Gamma_{\alpha} \mathbf{u}_{\alpha}, \quad \Gamma_{\alpha} = \sqrt{1 + \frac{p_{\alpha}^2}{mc^2}} = \frac{c}{\sqrt{c^2 - u_{\alpha}^2}}, \quad (B3) $$

$$ \mathbf{\nabla} \cdot \mathbf{E} = 4\pi \left( n_0 - \sum_{\alpha} n_{\alpha} \right), \quad \mathbf{\nabla} \cdot \mathbf{B} = 0, \quad (B4) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = -c \mathbf{\nabla} \times \mathbf{E}, \quad \frac{\partial \mathbf{E}}{\partial t} = c \mathbf{\nabla} \times \mathbf{B} - 4\pi c \left( e \sum_{\alpha} n_{\alpha} \mathbf{u}_{\alpha} \right). \quad (B5) $$

Linearization around an equilibrium with $\sum_{\alpha} n_{x,0} = n_0$ and $\sum_{\alpha} n_{x,0} \mathbf{u}_{x,0} = 0$ with $\mathbf{u}_{x,0} = \mathbf{u}_{x,0} \mathbf{e}_y$ for perturbations of the form $\delta f \sim (e^{i(kx-\omega t)} + e^{-i(kx-\omega t)})$, where indices 0 and 1 stand for equilibrium and first-order perturbations, leads again to $[\mathbf{e}] \cdot \mathbf{E}_1 = \mathbf{0}$, where the dispersion matrix $[\mathbf{e}]$ reads
$$ \begin{pmatrix} \omega^2 - \omega_2^2 - k^2 d_e^2 \left(1 + \frac{\omega_4^2}{\omega^2}\right) & kd_e \frac{\omega_3^2}{\omega} & 0 \\ kd_e \frac{\omega_3^2}{\omega} & \omega^2 - \omega_2^2 & 0 \\ 0 & 0 & \omega^2 - \omega_1^2 - k^2 d_e^2 \end{pmatrix}. \quad (B6) $$

Here, like in Refs. 46 and 50, we have introduced the characteristic frequencies,

$$ \begin{aligned} \omega_1^2 &= \frac{4\pi e^2}{m} \sum_{\alpha} \frac{n_{\alpha,0}}{\Gamma_{\alpha,0}}, & \omega_2^2 &= \frac{4\pi e^2}{m} \sum_{\alpha} \frac{n_{\alpha,0}}{\Gamma_{\alpha,0}^3}, \\ \omega_3^2 &= \frac{4\pi e^2}{m} \sum_{\alpha} \frac{n_{\alpha,0}}{\Gamma_{\alpha,0}} \frac{u_{\alpha,0}}{c}. \end{aligned} \quad (B7) $$

Note, however, the different choice of geometry with respect to the aforementioned articles, in which beams had been taken along 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 \to [4\pi e^2/(mc^2)] \sum_{\alpha} (1 - \Gamma_{\alpha,0})(n_{\alpha,0}u_{\alpha,0})$ as $\Gamma_{\alpha,0} \to 1$ for $(u_{\alpha,0}/c)^2 \to 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$.

One can then look at the nonlinear coupling when $\omega_3^2 = 0$, which in the most general case occurs in the non-relativistic limit, by retaining second-order terms (index 2) in the linearization of Eqs. (B1)–(B5) for $\Gamma_{\alpha,0} = 1$. In this non-relativistic limit, $\omega_1^2 = \omega_2^2 = \omega_p^2$. Differentiating once more Ampère's equation with respect to time and eliminating the magnetic field via Faraday's equation, one obtains the relation

$$ \frac{\partial^2 E}{\partial t^2} + \nabla(\nabla \cdot E) - \nabla^2 E = \sum_{\alpha} \left(\frac{\partial n_{\alpha}}{\partial t}\right) u_{\alpha} + \sum_{\alpha} n_{\alpha} \frac{\partial u_{\alpha}}{\partial t}. \quad (B8) $$

This, combined with momentum and continuity equations, leads to an equation for the $E_{x,2}$ component that can be written in the form

$$ \frac{\partial^2 E_{x,2}}{\partial t^2} - E_{x,2} = - \sum_{\alpha} n_{\alpha,0} u_{x,2,1} \frac{\partial}{\partial x} u_{x,2,1}. \quad (B9) $$

From Eqs. (B8) and (B9), one sees that the nonlinear coupling between e.m. and electrostatic perturbations in CFI can be indeed related to density fluctuations associated with Langmuir modes$^{28}$ 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 \dot{\varepsilon}(t) \cos(k_x x)$, with $\dot{\varepsilon}(t) = \varepsilon_0 e^{i\omega t}$. Saturation corresponds to $\varepsilon(t) \sim 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)$.

Direct integration of (8) yields $E_y(x, t) = -[B_M/(k_x c)] \dot{\varepsilon} \sin(k_x x) + E_y(0, t)$, which, combined with the second of Eq. (9) and with the boundary condition $E_y(0, t) = 0$, that we choose, gives

$$ J_y(x, t) = n_0 e k_x d_e^2 \Omega_c \left( \frac{\ddot{\varepsilon}}{k_x^2 c^2} + \varepsilon \right) \sin(k_x x), \quad (C1) $$

where $\Omega_c = eB_M/(mc)$ is a reference electron cyclotron frequency. Note that if we mean $B_z(x, 0)$ to correspond to a Weibel-like perturbation, $J_y(x, 0)$ given by Eq. (C1) is a perturbed quantity itself and, therefore, is not zero, initially. Defining then $\delta n(x, t)$ such that $n = n_0 + \delta n$, the first of Eqs. (9) and (10) imply, for example,

$$ X_{J_x}(x) = X_{E_x}(x), \quad 4\pi T_{J_x}(t) = -\dot{T}_{E_x}(t), \quad (C2) $$

$$ X'_{E_x}(x) = -X_{\delta n}(x), \quad T_{E_x}(t) = 4\pi e T_{\delta n}(t), \quad (C3) $$

where the prime and upper dot express derivation with respect to the space and the time variable, respectively.

In a volume $V_0 = L_x \times L_y \times L_z$, in which ions form a fixed neutralizing background, we then consider N electrons that would be at rest at $t=0$, if $\varepsilon|_{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, $\varepsilon(t) \neq 0$, the perturbation $\delta B(x, t)|_{t=0} = (0, 0, 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).

Using the symbol “...” for the single electron velocity $\tilde{v}(x, t)$ (note the difference with respect to the phase-space velocity coordinates $v$, independent from $x$, used previously), Newton equations of the single electron read

$$ \frac{\partial \tilde{v}_x}{\partial t} = \Omega_c \left\{ -\frac{cE_x}{B_M} - \varepsilon(t) \cos(k_x x) \tilde{v}_y \right\}, \quad (C4) $$

$$ \frac{\partial \tilde{v}_y}{\partial t} = \Omega_c \left\{ -\frac{cE_y}{B_M} + \varepsilon(t) \cos(k_x x) \tilde{v}_x \right\} \quad (C5) $$

and give a total energy equation in the form,

$$ \tilde{v}_x^2 + \tilde{v}_y^2 = -2 \frac{c \Omega_c}{B_M} \int_0^t (E_x \tilde{v}_x + E_y \tilde{v}_y) dt. \quad (C6) $$

In the expressions (C4)–(C6), the electric field components are to be treated as external, i.e., the effect of the particle dynamics on their evolution is neglected, since the so-far missing information on the electron distribution prevents us to relate the right-hand side of Eqs. (C6) and (7). However, the feedback that the electron motion has on the e.m. fields via Eqs. (8)–(10) suggests us to look for perturbative solutions of Eqs. (C4) and (C5), in which connection is made between $\tilde{v}$ and $J$, and thus between the e.m. fields in Newton's law and the particle population: writing $\tilde{v}_x = \delta \tilde{v}_x$, $\tilde{v}_y = \delta \tilde{v}_y$, $n = n_0 + \delta n$, and $J = \delta J$, $B = \delta B$, $E = \delta E$, where all quantities with $\delta$ 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, \dot{\epsilon}\Delta t, \ddot{\epsilon}(\Delta t)^2\} \ll 1$ (i.e., for $\Delta t \ll \min\{\epsilon_0/\gamma_B, \sqrt{\epsilon_0}/\gamma_B\}$). In terms of the separation of variables used before, any scalar function $F(x, t)$ has $X_F(x) \sim O(1)$ while $T_F(t)$ can be ordered like a positive power of $\epsilon$.

Taking $n_0 \sim O(1)$, one sees from Eqs. (C1), (C4), and (C5) that

$$\delta B_z \sim \delta E_y \sim \delta J_y \sim \epsilon, \quad \delta \bar{v}_y \sim \epsilon, \qquad (C7)$$

while Eqs. (C2)–(C4) combined with (13) imply

$$\delta \bar{v}_x \sim \delta J_x \sim \delta E_x \sim \delta n. \qquad (C8)$$

Direct integration of (C5) performed while neglecting $\sim \epsilon^2$ terms and using previous results on $E_y(x, t)$ yields

$$\delta \bar{v}_y(x, t) \simeq + \frac{\Omega_c}{k_x} \epsilon(t) \sin(k_x x) + O(\epsilon^2), \qquad (C9)$$

while Eq. (C4) becomes

$$\frac{\partial}{\partial t} \delta \bar{v}_x = -\Omega_c \frac{\delta E_x}{B_M} - \frac{\Omega_c^2}{2k_x} \epsilon^2(t) \sin(2k_x x) + O(\epsilon^3). \qquad (C10)$$

The simplistic assumptions done on the initial electron distribution allow us to approximate the perturbed current density $J$, statistically defined by Eq. (6), as $\delta J \simeq -en_0\delta\bar{v} + O(\delta n\delta\bar{v})$, 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

$$\frac{\partial}{\partial t} \delta n \simeq -n_0 \frac{\partial}{\partial x} \delta \bar{v}_x. \qquad (C11)$$

This implies

$$X_{\delta n}(x) = -X'_{\nu_x}(x), \quad \dot{T}_{\delta n}(t) = n_0 T_{\nu_x}(t), \qquad (C12)$$

which, combined with Eqs. (C2) and (C3) and with the spatial derivative of Eq. (C10), leads, after little algebra, to two conditions on $X_{\delta n}$ and $T_{\delta n}$,

$$X_{\delta n}(x) = -\cos(2k_x x), \qquad (C13)$$

$$\ddot{T}_{\delta n}(t) + \omega_{pe}^2 T_{\delta n}(t) = \Omega_c^2 \epsilon_0^2 e^{\gamma_B t}. \qquad (C14)$$

Using the trivial exponential solution admitted by Eq. (C14), one finally obtains

$$\delta n(x, t) = - \frac{\Omega_c^2}{\gamma_B^2 + \omega_{pe}^2} n_0 \epsilon^2(t) \cos(2k_x x), \qquad (C15)$$

from which the ordering

$$\delta \bar{v}_x \sim \delta J_x \sim \delta E_x \sim \delta n \sim \epsilon^2 \qquad (C16)$$

follows, and from which the analytical expressions of $\delta \bar{v}_x$, $\delta J_x$ and $\delta E_x$ can be deduced by using Eqs. (C2), (C3), and (C12).

REFERENCES

¹H. Schamel, “Hole equilibria in Vlasov-Poisson systems: A challenge to wave theories of ideal plasmas,” Phys. Plasmas 7, 4831 (2000).

²B. Afeyan, K. Won, V. Savchenko, T. W. Johnston, A. Ghizzo, and P. Bertrand, “Kinetic electrostatic electron nonlinear (KEEN) waves and their interaction driven by the ponderomotive force of crossing laser beams,” arXiv:1210.8105 (2003).

³M. Albrecht-Marc, A. Ghizzo, T. Réveillé, P. Betrand, and T. W. Johnston, “Semi-Lagrangian Vlasov-Maxwell simulation of self-sustained kinetic electron nonlinear electrostatic waves in relativistic laser-plasma interaction,” Transport Theory Stat. Phys. 34, 127 (2005).

⁴F. Valentini, T. M. O’Neil, and D. H. E. Dubin, “Excitation of nonlinear electron acoustic waves,” Phys. Plasmas 13, 052303 (2006).

⁵T. W. Johnston, Y. Tyshetskiy, A. Ghizzo, and P. Bertrand, “Persistent sub-plasma frequency kinetic electrostatic electron nonlinear waves,” Phys. Plasmas 16, 042105 (2009).

⁶H. Schamel, “Analytic BGK modes and their modulational instability,” J. Plasma. Phys. 13, 139 (1975).

⁷H. Schamel, “Cnoidal electron hole propagation: Trapping, the forgotten non-linearity in plasma and fluid dynamics,” Phys. Plasmas 19, 020501 (2012).

⁸F. Valentini, D. Perrone, F. Califano, P. Veltri, F. Pegoraro, P. Veltri, P. J. Morrison, and T. M. O’Neil, “Undamped electrostatic plasma waves,” Phys. Plasmas 19, 092103 (2012).

⁹H. Schamel, “Comment on ‘Undamped electrostatic plasma waves’ [Phys. Plasmas 19, 092103 (2012)],” Phys. Plasmas 19, 092013 (2012).

¹⁰F. Valentini, D. Perrone, F. Califano, P. Veltri, F. Pegoraro, P. Veltri, P. J. Morrison, and T. M. O’Neil, “Response to ‘Undamped electrostatic plasma waves’” [Phys. Plasmas 20, 034701 (2013)],” Phys. Plasmas 20, 034702 (2013).

¹¹A. Ghizzo and D. Del Sarto, “Nonlinear nature of kinetic undamped waves induced by electrostatic turbulence in stimulated Raman scattering,” Eur. Phys. J. D 68, 275 (2014).

¹²A. Schekochihin, A. A. Cowley, S. C. Dorland, G. W. Hammett, G. G. Howes, E. Quataert, and T. Tatusno, “Astrophysical gyrokinetics: Kinetic and fluid theory of turbulent cascades in magnetized weakly collisional plasmas,” Astrophys. J. Suppl. Ser. 182, 310 (2009).

¹³E. Papini, A. Cicone, L. Franci, M. Piersanti, S. Landi, P. Hellinger, and A. Verdini, “Spacetime Hall-MHD turbulence at sub-ion scales: Structures or waves?,” Astrophys. J. Lett. 917, L12 (2021).

¹⁴A. Mallet, S. Dorfman, M. Abler, T. Bowen, and C. H. K. Chen, “Nonlinear dynamics of small-scale Alfvén waves,” arXiv:2303.10192 (2023).

¹⁵D. Del Sarto, F. Pegoraro, and F. Califano, “Pressure anisotropy and small spatial scales induced by velocity shear,” Phys. Rev. E 93, 05303 (2016).

¹⁶D. Del Sarto and F. Pegoraro, “Shear-induced pressure anisotropization and correlation with fluid vorticity in a low collisionality plasma,” Monthly Notices R. Astron. Soc. 475, 181 (2018).

¹⁷P. A. Cassak and M. H. Barbhuiya, “Pressure-strain interaction. I. On compression, deformation, and implications for Pi-D,” Phys. Plasmas 29, 122306 (2022).

¹⁸Y. Yang, W. H. Matthaeus, T. N. Parashar, P. Wu, M. Wan, Y. Shi, S. Chen, V. Roytershteyn, and W. Daughton, “Energy transfer channels and turbulence cascade in Vlasov-Maxwell turbulence,” Phys. Rev. E 95, 061201 (2017).

¹⁹Y. Yang, W. H. Matthaeus, T. N. Parashar, C. C. Haggerty, V. Roytershteyn, W. Daughton, M. Wan, Y. Shi, and S. Chen, “Energy transfer, pressure tensor, and heating of kinetic plasma,” Phys. Plasmas 24, 072306 (2017).

²⁰W. H. Matthaeus, Y. Yang, M. Wan, T. N. Parashar, R. Bandyopadhyay, A. Chasapis, O. Pezzi, and F. Valentini, “Pathways to dissipation in weakly collisional plasmas,” Astrophys. J. 891, 101 (2020).

²¹P. A. Cassak, M. H. Barbhuiya, H. Liang, and M. R. Argall, “Quantifying energy conversion in higher-order phase space density moments in plasmas,” Phys. Rev. Lett. 130, 085201 (2023).

²²P. Hellinger, V. Montagud-Camps, L. Franci, L. Matteini, E. Papini, A. Verdini, and S. Landi, “Ion-scale transition of plasma turbulence: Pressure-strain effect,” Astrophys. J. 930, 48 (2022).

²³A. I. Smolyakov, P. H. Diamond, and M. V. Medvedev, “Role of ion diamagnetic effects in the generation of large scale flows in toroidal ion temperature gradient mode turbulence,” Phys. Plasmas 7, 3987 (2000).
<footnote n="24">²⁴A. Ghizzo and F. Palermo, “Shear-flow trapped-ion-mode interaction revisited. I. Influence of low-frequency zonal flow on ion-temperature-gradient driven turbulence,” Phys. Plasmas 22, 082303 (2015).</footnote>

<footnote n="25">²⁵E. S. Weibel, “Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution,” Phys. Rev. Lett. 2, 83 (1959).</footnote>

<footnote n="26">²⁶K. A. Taggart, B. B. Godrey, C. E. J. Rhoades, and H. C. Ives, “Second-order effects in the Weibel instability,” Phys. Rev. Lett. 29, 1729 (1972).</footnote>

<footnote n="27">²⁷A. Stockem, M. E. Dieckmann, and R. Schlickeiser, “Pic simulations of the thermal anisotropy-driven Weibel instability: Field growth and phase space evolution upon saturation,” Plasma Phys. Controlled Fusion 51, 075041 (2009).</footnote>

<footnote n="28">²⁸L. Palodhi, F. Califano, and F. Pegoraro, “Nonlinear kinetic development of the Weibel instability and the generation of electrostatic coherent structures,” Plasma Phys. Controlled Fusion 51, 125006 (2009).</footnote>

<footnote n="29">²⁹M. E. Innocenti, M. Lazar, S. Markidis, G. Lapenta, and S. Poedts, “Electron streams formation and secondary two stream instability onset in the post-saturation regime of the classical Weibel instability,” Phys. Plasmas 18, 052140 (2011).</footnote>

<footnote n="30">³⁰F. Califano, F. Pegoraro, S. V. Bulanov, and A. Mangeney, “Kinetic saturation of the Weibel instability in a collisionless plasma,” Phys. Rev. E 57, 7048 (1998).</footnote>

<footnote n="31">³¹F. Califano, T. Cecchi, and C. Chiuderi, “Nonlinear kinetic regime of the Weibel instability,” Phys. Plasmas 9, 451 (2002).</footnote>

<footnote n="32">³²M. Honda, J. Meyer-ter Vehn, and A. Pukhov, “Two-dimensional particle-in-cell simulation for magnetized transport of ultra-high relativistic currents in plasma,” Phys. Plasmas 7, 1302 (2000).</footnote>

<footnote n="33">³³M. E. Dieckmann, “The filamentation instability driven by warm electron beams: Statistics and electric field generation,” Plasma Phys. Controlled Fusion 51, 124042 (2009).</footnote>

<footnote n="34">³⁴A. Bret, L. Gremillet, and M. E. Dieckmann, “Multidimensional electron beam-plasma instabilities in the relativistic regime,” Phys. Plasmas 17, 120501 (2010).</footnote>

<footnote n="35">³⁵A. Grassi, M. Grech, F. Amiranoff, F. Pegoraro, A. Macchi, and C. Riconda, “Electron Weibel instability in relativistic counter-streaming plasmas with flow-aligned external magnetic fields,” Phys. Rev. E 95, 023203 (2017).</footnote>

<footnote n="36">³⁶B. D. Fried, “Mechanism for instability of transverse plasma waves,” Phys. Fluids 2, 337 (1959).</footnote>

<footnote n="37">³⁷A. Inglebert, T. Ghizzo, A. Réveillé, D. Del Sarto, P. Bertrand, and F. Califano, “A multi-stream Vlasov modeling unifying relativistic Weibel-type instabilities,” Euro. Phys. Lett. 95, 45002 (2011).</footnote>

<footnote n="38">³⁸F. Califano, R. Prandi, F. Pegoraro, and S. V. Bulanov, “Magnetic-field generation and wave-breaking in collisionless plasmas,” Phys. Plasmas 60, 331 (1998).</footnote>

<footnote n="39">³⁹M. Sarrat, D. Del Sarto, and A. Ghizzo, “Fluid description of Weibel-type instabilities via full pressure tensor dynamics,” Europhys. Lett. 115, 45001 (2016).</footnote>

<footnote n="40">⁴⁰M. Sarrat, D. Del Sarto, and A. Ghizzo, “A pressure tensor description for the time-resonant Weibel instability,” J. Plasma Phys. 83, 705830103 (2017).</footnote>

<footnote n="41">⁴¹A. Ghizzo, M. Sarrat, and D. Del Sarto, “Vlasov models for kinetic Weibel-type instabilities,” J. Plasma Phys. 83, 705830101 (2017).</footnote>

<footnote n="42">⁴²P. Bertrand, M. Albrecht-Marc, T. Réveillé, and A. Ghizzo, “Vlasov models for laser-plasma interaction,” Transport Theory Stat. Phys. 34, 103 (2005).</footnote>

<footnote n="43">⁴³M. Sarrat, A. Ghizzo, D. Del Sarto, and L. Sarrat, “Parallel implementation of a relativistic semi-Lagrangian Vlasov-Maxwell solver,” Europhys. J. D 71, 271 (2017).</footnote>

<footnote n="44">⁴⁴G. A. Askar'Yan, S. V. Bulanov, F. Pegoraro, and A. M. Pukhov, “Magnetic interaction of self-focusing channels and fluxes of electromagnetic radiation, their coalescence, the accumulation of energy, and the effect of external magnetic fields,” JETP Lett. 60, 251 (1994); available at http://jetpletters.ru/ps/0/article_20333.shtml.</footnote>

<footnote n="45">⁴⁵M. Sarrat, “Physique des instabilités de type Weibel,” Ph.D. thesis (Université de Lorraine, 2017).</footnote>

<footnote n="46">⁴⁶F. Pegoraro, S. V. Bulanov, F. Califano, and M. Lontano, “Nonlinear development of the Weibel instability and magnetic field generation in collisionless plasmas,” Phys. Scr. T63, 262 (1996).</footnote>

<footnote n="47">⁴⁷B. Basu, “Moment equation description of Weibel instability,” Phys. Plasmas 9, 5131 (2002).</footnote>

<footnote n="48">⁴⁸We recall that the conservation of the distribution function of a species of charge *q* in a magnetic field associated with a vector potential *A*(*x*, *t*) in the absence of collisions is established in a phase-space where the canonical conjugated variables are (*x*, *P*). The passage from the canonical momentum *P*</footnote>

<footnote>= *p* + *qA*(*x*, *t*) to the kinematic momentum *p* that allows the immediate rewriting of Vlasov equation for *f*(*x*, *p*, *t*) as a Vlasov equation for *f*(*x*, *p*, *t*) in the “phase-space” (*x*, *p*) (or for *f*(*x*, *v*, *t*)) in the (*x*, *v*) configuration space for the non-relativistic limit *p* = *mv* is allowed by the equivalence of the infinitesimal volumes *dxdP* = *dxdp*, which grants the validity of Liouville’s theorem when performing the change of variables *P* → *p* (cf., e.g., Ref. 45, §2.2.1.3).</footnote>

<footnote n="49">⁴⁹A. Ghizzo and P. Bertrand, “On the multistream approach of relativistic Weibel instability. I. Linear analysis and specific illustrations,” Phys. Plasmas 54, 082109 (2013).</footnote>

<footnote n="50">⁵⁰F. Califano, F. Pegoraro, and S. V. Bulanov, “Spatial structure and time evolution of the Weibel instability in collisionless inhomogeneous plasmas,” Phys. Rev. E 56, 963 (1997).</footnote>

<footnote n="51">⁵¹A. Ghizzo, “On the multistream approach of relativistic Weibel instability. II. Bernstein-Green-Kruskal-type waves in magnetic trapping,” Phys. Plasmas 20, 082110 (2013).</footnote>

<footnote n="52">⁵²A. Ghizzo, “On the multistream approach of relativistic Weibel instability. III. Comparison with full-kinetic Vlasov simulations,” Phys. Plasmas 20, 082111 (2013).</footnote>

<footnote n="53">⁵³A. Ghizzo, D. Del Sarto, and M. Sarrat, “Low- and high-frequency nature of oblique filamentation modes. I. Linear theory,” Phys. Plasmas 27, 072103 (2020).</footnote>

<footnote n="54">⁵⁴S. Kingsep, K. V. Chukbar, and V. V. Yan'kov, “Electron magnetohydrodynamics,” Phys. Rep. 243, 215-315 (1990).</footnote>

<footnote n="55">⁵⁵D. Del Sarto, A. Tenerani, and F. Pegoraro, “Magneto-elastic’ waves in an anisotropic magnetized plasma,” Plasma Phys. Controlled Fusion 59, 045002 (2017).</footnote>

<footnote n="56">⁵⁶A. Macmahon, “Finite gyro-radius corrections to the hydromagnetic equations for a Vlasov plasma,” Phys. Fluids 8, 1840 (1965).</footnote>

<footnote n="57">⁵⁷A. B. Mikhailovskii and A. I. Smolyakov, “Theory of low-frequency magneto-sonic solitons,” Zh. Eksp. Teor. Fiz. 88, 189 (1985); available at http://www.jetp.ras.ru/cgi-bin/dn/e_061_01_0109.pdf.</footnote>

<footnote n="58">⁵⁸D. Del Sarto, F. Pegoraro, and F. Califano, “Electron parallel compressibility in the nonlinear development of two-dimensional collisionless magnetohydrodynam-ic reconnection,” Mod. Phys. Lett. B 20, 931-961 (2006).</footnote>

<footnote n="59">⁵⁹M. G. Haines, “Magnetic-field generation in laser fusion and hot-electron transport,” Can. J. Phys. 64, 912-919 (1986).</footnote>

<footnote n="60">⁶⁰M. G. Haines, “Saturation mechanisms for the generated magnetic field in non-uniform laser-matter irradiation,” Phys. Rev. Lett. 78, 254 (1997).</footnote>

<footnote n="61">⁶¹L. Biermann, “Über den Ursprung der Magnetfelder auf Sternen und im interstellaren Raum,” Z. Naturforsch. 5a, 65 (1950).</footnote>

<footnote n="62">⁶²K. M. Schoeffler, N. F. Loureiro, R. A. Fonseca, and L. O. Silva, “The generation of magnetic fields by the Biermann battery and the interplay with the Weibel instability,” Phys. Plasmas 23, 056304 (2016).</footnote>

<footnote n="63">⁶³A. Inglebert, A. Ghizzo, T. Réveillé, P. Bertrand, and F. Califano, “Electron temperature anisotropy instabilities represented by superposition of streams,” Plasma Phys. Controlled Fusion 19, 122109 (2012).</footnote>

<footnote n="64">⁶⁴M. Navet and P. Bertrand, “‘Multiple water-bag’ model and Landau damping,” Phys. Lett., A 34, 117-118 (1971).</footnote>

<footnote n="65">⁶⁵P. Bertrand and M. R. Feix, “Non linear electron plasma oscillation: The ‘water bag’ model,” Phys. Plasmas 28, 68 (1968).</footnote>

<footnote n="66">⁶⁶P. H. Yoon and R. C. Davidson, “Exact analytical model of the classical Weibel instability in a relativistic anisotropic plasma,” Phys. Rev. A 35, 2718 (1987).</footnote>

<footnote n="67">⁶⁷P. Bertrand, D. Del Sarto, and A. Ghizzo, The Vlasov Equation 1. History and General Properties (ISTE and Wiley, London, 2019).</footnote>

<footnote n="68">⁶⁸M. Lazar, M. E. Dieckmann, and S. Poedts, “Resonant Weibel instability in counterstreaming plasmas with temperature anisotropies,” J. Plasma Phys. 76, 49 (2010).</footnote>

<footnote n="69">⁶⁹L. Franci, E. Papini, A. Micera, G. Lapenta, P. Hellinger, D. Del Sarto, D. Burgess, and S. Landi, “Anisotropic electron heating in turbulence-driven magnetic reconnection in the near-sun solar wind,” Astrophys. J. 936, 27 (2022).</footnote>

<footnote n="70">⁷⁰Y. Yang, W. H. Matthaeus, S. Roy, V. Roytershteyn, T. N. Parashar, R. Bandyopadhyay, and M. Wan, “Pressure-strain interaction as the energy dissipation estimate in collisionless plasma,” Astrophys. J. 929, 142 (2022).</footnote>

<footnote n="71">⁷¹Y. Yang, F. Pecora, W. H. Matthaeus, S. Roy, M. E. Cuesta, A. Chasapis, T. Parashar, R. Bandyopadhyay, D. J. Gershman, B. L. Giles, and J. L. Buchr,</footnote>
"Quantifying the agyrotropy of proton and electron heating in turbulent plasmas," *Astrophys. J.* **944**, 148 (2023).

<footnote n="72">⁷²A. Bret, M.-C. Firpo, and C. Deutsch, "Characterization of the initial filamentation of a relativistic electron beam passing through a plasma," *Phys. Rev. Lett.* **94**, 115002 (2005).</footnote>

<footnote n="73">⁷³A. Bret, M.-C. Firpo, and C. Deutsch, "Electromagnetic instabilities for relativistic beam-plasma interaction in whole k space: Nonrelativistic beam and plasma temperature effects," *Phys. Rev. E* **72**, 016403 (2005).</footnote>

<footnote n="74">⁷⁴M. Tzoufras, C. Ren, F. Tsung, J. W. Tonge, W. B. Mori, M. Fiore, R. A. Fonseca, and L. O. Silva, "Space-charge effects in the current-filamentation or Weibel instability," *Phys. Rev. Lett.* **96**, 105002 (2006).</footnote>

<footnote n="75">⁷⁵A. Achterberg and J. Wiersma, "The weibel instability in relativistic plasmas—I. Linear theory," *Astron. Astrophys.* **475**, 1–18 (2007).</footnote>

<footnote n="76">⁷⁶A. Bret, L. Gremillet, and J. C. Bellido, "How really transverse is the filamentation instability?," *Phys. Plasmas* **14**, 032103 (2007).</footnote>

<footnote n="77">⁷⁷R. Kumar, D. Eichler, and M. Gedalin, "Electron heating in a relativistic, Weibel-unstable plasma," *Astrophys. J.* **806**, 165 (2015).</footnote>

<footnote n="78">⁷⁸A. Ghizzo and D. Del Sarto, "Momentum transfer driven by fluctuations in relativistic counter-propagating electron beams," *Plasma Phys. Controlled Fusion* **61**, 055007 (2021).</footnote>

<footnote n="79">⁷⁹P. H. Yoon, "Relativistic Weibel instability," *Phys. Plasmas* **14**, 024504 (2007).</footnote>

<footnote n="80">⁸⁰H. H. Kaang, C.-M. Ryu, and P. H. Yoon, "Nonlinear saturation of relativistic Weibel instability driven by thermal anisotropy," *Phys. Plasmas* **16**, 082103 (2009).</footnote>

<footnote n="81">⁸¹M. Lazar, R. Schlickeiser, and S. Poedts, "On the existence of Weibel instability in a magnetized plasma. I. Parallel wave propagation," *Phys. Plasmas* **16**, 012106 (2009).</footnote>

<footnote n="82">⁸²M. Lazar, R. Schlickeiser, R. Wielebinski, and S. Poedts, "Cosmological effects of Weibel-type instabilities," *Astrophys. J.* **693**, 1133 (2009).</footnote>

<footnote n="83">⁸³M. Ghorbanalilu, S. Sadegzadeh, Z. Ghaderi, and A. R. Niknam, "Weibel instability for a streaming electron, counterstreaming e-e, and e-p plasmas with intrinsic temperature anisotropy," *Phys. Plasmas* **21**, 052102 (2014).</footnote>

<footnote n="84">⁸⁴L. Gremillet, D. Bénisti, E. Lefebvre, and A. Bret, "Linear and nonlinear development of oblique beam-plasma instabilities in the relativistic kinetic regime," *Phys. Plasmas* **14**, 040704 (2007).</footnote>

<footnote n="85">⁸⁵A. Bret, L. Gremillet, D. Bénisti, and E. Lefebvre, "Exact relativistic kinetic theory of an electron-beam-plasma system: Hierarchy of the competing modes in the system-parameter space," *Phys. Rev. Lett.* **100**, 205008 (2008).</footnote>

<footnote n="86">⁸⁶V. Skoutnev, A. Hakim, J. Juno, and J. M. TenBarge, "Temperature-dependent saturation of Weibel-type instabilities in counter-streaming plasmas," *Astrophys. J. Lett.* **872**, L28 (2019).</footnote>

<footnote n="87">⁸⁷A. Ghizzo and D. Del Sarto, "Low- and high-frequency nature of oblique filamentation modes. II. Vlasov-Maxwell simulations of collisionless heating process," *Phys. Plasmas* **27**, 072104 (2020).</footnote>

<footnote n="88">⁸⁸A. Ghizzo, D. Del Sarto, and H. Betar, "Collisionless heating driven by Vlasov filamentation in a counterstreaming beams configuration," *Phys. Rev. Lett.* **131**, 035101 (2023).</footnote>

<footnote n="89">⁸⁹B. D. Fried and S. D. Conte, *The Plasma Dispersion Function. The Hilbert Transform of the Gaussian* (Elsevier, 1961).</footnote>