Reversible energy conversion between magnetic and kinetic energies has been recently demonstrated in a system of counterstreaming electron beams [see A. Ghizzo et al., Phys. Rev. Lett. 131, 035101 (2023)]. During the first step of the instability, the growth of a currentdriven filamentation magnetic field is observed when propagative oblique solutions are considered, followed by the reversal of energy transfer from magnetic to kinetic energy in a second step. This highlights a new physical mechanism of the Vlasov equation: the enhancement of filamentation of the distribution function in the presence of the phase synchronization of the Van Kampen eigenmodes. This gives rise to a bifurcation toward selforganization and to a strong plasma heating. This new plasma heating mechanism possibly provides a new perspective on the role played by the filamentation in phase space in the relativistic regime of Weibeltype instabilities.
I. INTRODUCTION
In a configuration characterized by two counterpropagating electron beams, Weibeltype instabilities can couple with electrostatic unstable modes by generating the socalled oblique instabilities (OIs). Recently, in Ref. 1, we have shown via kinetic simulations—performed with both semiLagrangian (SL) and particleincell (PIC) Vlasov–Maxwell codes—how the spatial filamentation induced by these modes is capable of affecting, both linearly and nonlinearly, the phasespace filamentation of the distribution function by synchronizing the phases of the Van Kampen modes. In this work, we further discuss the role played by phasespace filamentation in this process.
Two filamentation mechanisms are involved: on the one hand, the filamentation of the distribution function in the velocity space, an intrinsic property of the Vlasov equation, and on the other hand, the spatial filamentation of currents, related to the current filamentation instability (CFI), a kind of Weibel instability. This leads to largeamplitude fluctuations of the distribution function in the phase space, to which hereafter we refer as to a “reinforced filamentation.” This process may lead to a change in the saturation regime of the OIs, where the energy stored in the magnetic field can be transferred back to plasma particles via a kinetic heating mechanism first discussed in Ref. 1. Such a process differs from the magnetic reconnection, which is forbidden in the $ 2 D 2 V$ geometry considered in that work. A global reorganization takes place in which collisionless wave–particle interactions transfer the energy stored in the beams and in the magnetic field, into the internal energy of the plasma.
The amplification of magnetic fields induced by temperature or momentum anisotropy and the conversion of the magnetic energy into kinetic or internal energy represent two complementary features of the dynamics of a magnetic field in a collisionless plasma. The first process is exemplified by Weibeltype instabilities. Previous linear studies of Weibeltype instabilities^{2–14} have identified different classes of electron beamplasma instabilities: the current filamentation instability^{15–19} (CFI), the Weibel instability (WI), the electrostatic twostream instability (TSI), and finally the oblique instability (OI).^{20,21} The special case of two counterstreaming electrons beams is particularly relevant to astrophysical plasmas, e.g., in the gammaray burst production scenario^{22} or in interpenetrating plasma flows.^{23,24} WI is driven by thermal anisotropy in the plasma provided that the perpendicular temperature $ T \u22a5$ is larger than the parallel temperature ( $ T \u2225$ here), where the symbols $\u22a5$ and $\u2225$ denote the perpendicular and parallel direction with respect the wave vector k. While TSI modes are longitudinal (i.e., with $k$ aligned to the beams), CFI modes are perpendicular (i.e., with $k$ perpendicular to the beams). In these systems, we usually observe the conversion of kinetic energy into magnetic energy.
The reverse process of energy transfer where the magnetic energy is converted into kinetic energy through acceleration or heating of charged particles is usually met in magnetic reconnection.^{25–33} Such a process involves a topology change of field lines, which leads to a new equilibrium configuration of lower magnetic energy. In ideal magnetohydrodynamics (MHD), magnetic field lines consistently move with the plasma due to the frozenin flux principle, ensuring they cannot break or tear apart. However, satellite observations (see, e.g., Refs. 24 and 34) indicate the occurrence of rapid global magnetic reconnection processes, suggesting the presence of anomalous efficient dissipation of magnetic energy. How such physical processes occur and how fast line breaking takes place are today open questions.
These two energy transfer mechanisms are not only complementary but strongly correlated with turbulence and with associated direct cascades, magnetic reconnection, and collisionless plasma heating.^{34–45}
This paper is devoted to extending ideas about the turbulence cascade into the deep kinetic regime. The importance of the energy conversion implied by magnetic reconnection in the turbulence cascade has been already pointed out by several authors,^{46–49} so as the role played by phasemixing in Ref. 50. In this context, it is of interest to develop a better understanding of the basic mechanisms of kinetic heating, which encompasses processes of conversion from magnetic into kinetic energy that are more fundamental than magnetic reconnection and that may become important in kinetic plasma turbulence. This study is based on a twodimensional (2D) configuration consisting of two counterpropagating electrons beams with a linearly polarized electromagnetic field. Such a configuration is unstable to CFI, and it can develop a turbulent cascade toward large wave vectors. In order to take into account this turbulent cascade in a selfconsistent way with the CFI instability, we focus on a particular class of these instabilities: the oblique modes.
Recently, in Refs. 51 and 52, two different classes of OI modes have been investigated by solving the corresponding dispersion relation: the first one, referred here as “the nonpropagative branch,” exhibits the generation and the growth of a quasistatic magnetic field from CFI, in which the kinetic energy is converted into magnetic energy. However, by solving numerically the dispersion relation, a second branch of solutions has been identified in Ref. 51. This is the propagative branch of oblique modes. While the beginning of the dynamics, driven by this propagative OI mode, is very similar to the dynamics of the nonpropagative branch, a transition toward a different dynamics is observed in its nonlinear evolution. When several propagative OI modes are excited, this dynamics leads to the emergence of a directlike cascade to smaller scales. Here, we will focus on the physical mechanism that leads to the reversion of the energy transfer, associated with this directlike cascade. This mechanism, first proposed in Ref. 1, is linked to the emergence of a spatial filamentation of the distribution function. This is a direct result of the turbulent cascade, which produces thinner and thinner filaments in the configuration space. In the strong relativistic regime, such a process induces an enhancement of the filamentation of the distribution function (here referred to as the “reinforced” filamentation). This mechanism can induce deep modifications in the energy transfer and can lead to a strong stochastic heating. The emergence of these effects introduces asymmetries (elements of irreversibility) between two fundamental aspects of the Vlasov model: the filamentation of the distribution function in links with the time reversibility of the Vlasov equation and the phase synchronization mechanism of Van Kampen modes.^{53} The analysis of phase synchronization of oscillators (Kuramoto's approach) can help to identify and characterize these aspects.
The formulation of the Vlasov–Maxwell system as a (nonlinear) Hamiltonian Kuramoto system, involving differential equations, highlights the role of Van Kampen mode synchronization in the heating process. A linear analysis shows that the amount of phasestrophy (a kind of enstrophy extended to phase space) plays the role of an “energy” in the linear regime, and demonstrates a connection between mode synchronization and filamentation of the distribution function. Quantifying the phasestrophy flux density represents however a useful analysis tool even in the nonlinear regime: it provides an estimate of the importance of the synchronization process (cf. also Ref. 1) and displays a secular growth when Van Kampen mode synchronization becomes global.
This paper addresses the statistical properties of energy transfer across scales, recognizing the possible significance of a direct cascade, induced by nonlinear interactions, as well as the role of phase synchronization of Van Kampen modes at the kinetic scales. Phasespace filamentation driven by (global) phase synchronization can trigger an alternative channel of the conversion between kinetic and magnetic energy. With respect to our previous work,^{1} we investigate here in greater detail the features of the phasespace filamentation involved in this complex mechanism.
The paper is organized as follows. In Sec. II, we briefly recall the Vlasov–Maxwell (VM) model and we discuss in Sec. III several of its fundamentals aspects, mainly those related to the concepts of phase synchronization of Van Kampen modes. We also discuss the notion of “the enhanced” filamentation process of the distribution function in the presence of transition toward global phase synchronization. Section IV shows the main differences between freestreaming and reinforced filamentation. Section V presents the numerical particleincell and semiLagrangian (Vlasov) schemes and discusses the results of the study of two classes of oblique Weibeltype modes, propagating and nonpropagating. Section VI is devoted to the discussion of a new collisionless plasma heating scenario that originates from the growth of propagative oblique Weibel modes in the presence of large fluctuations. Finally, conclusions are presented in Sec. VII.
II. THE VLASOV–MAXWELL MODEL AND THE PHASESPACE FILAMENTATION
A. The physical model
B. Filamentation of f in the phase space vs spatial filamentation of currents
One of the fundamental properties of the Vlasov equation (1) concerns the conservation of any functional of f, noted $ G ( f )$, $G$ being an arbitrary function of f. This conservation leads to the socalled Casimir invariants $ C ( f ) = \u222b d 3 x d 3 p G ( f ) = const$. Their conservation follows from the invariance of the Vlasov equation under time reversal and from the fact that $ G ( f )$ commutes with the Hamiltonian. From a mathematical point of view, the asymptotic convergence of f holds only in the “weak” sense (see Refs. 54–56) and the derivatives of the distribution f with respect the velocity variables grow quickly in time. Indeed, linear Landau damping suggests that a “homogenization” phenomenon takes place leading to an entropypreserving kinetic dissipation where kinetic fast oscillations globally compensate each other in an averaging process in phase space. This leads to the phenomenon of the filamentation of f in the velocity space (or equivalently the momentum space), which is one aspect of filamentation, namely, the “kinematic filamentation” of the distribution function in the velocity space, one of the fundamental properties of the Vlasov equation. This kinematic feature can be modified or even amplified by CFI. CFI induces a filamentation of current densities in configuration space, which also leads to the production of thinner and thinner filaments of the current densities and, consequently, of the distribution function itself.
While the generation of thinner filaments goes on, the information transfer from the small to the large wavenumbers follows from the energy transfer between scales associated with this filamentation process. Information is usually conserved in a continuum velocity space (i.e., when the size of the elementary cell tends to zero) and the entropy is exactly conserved. It is such a mathematical property that “guarantees” the time reversibility of the Vlasov equation in the kinetic scales, as met for instance in the echoes problem in Ref. 57.
Furthermore, from a physical point of view, the kinematic filamentation is a real physical phenomenon in plasmas that might be responsible of the collisionless plasma heating as it, e.g., observed in the Earth's magnetosheath in Refs. 34 and 45 or in interpenetrating plasmas in Ref. 24. From Eq. (1), it can be seen that the value of f cannot change along a characteristic. Often the characteristics of the Vlasov equation mix together phasespace regions where the value of f is significantly different. Steep gradients of f are this way generated. In simulations that are performed on a (Eulerian) fixed grid in phase space, the grid almost inevitably becomes too coarse as this fine graining develops, which leads to the growth of the entropy. From a numerical point of view, this “filamentation” problem was well known in Eulerian Vlasov simulation and has been often discussed in Refs. 58–63.
On the another hand, oblique Weibeltype instabilities can induce smallscale, largeamplitude fluctuations in the distribution function, leading to what we refer to as “dynamical filamentation.” This phenomenon may alter the typical filamentation property of the Vlasov equation in velocity space. Consequently, the “kinematic filamentation” is enhanced, and we can correlate it with the variation in the phasestrophy flux within the Kuramoto description.
III. SYNCHRONIZATION VERSUS FILAMENTATION IN THE VLASOV FORMALISM
A. Synchronization aspects in the VM system
In Eq. (11), the first term in the righthand side of the equation corresponds to the intrinsic filamentation process of f in the momentum space (coming from the freestreaming term), which is related to the Van Kampen modes.^{53} This term plays the role of the natural frequency ω_{n} of the nth oscillator in Eq. (A3) of the Hamiltonian formulation of the Kuramoto model (see Appendix A). In a similar way, the parameter Kr, product of the coupling strength K by the order parameter r (cf. the Hamilton–Jacobi equation (A2) in the Appendix A), is replaced by the term $  R k ( p , t )  /  f k ( p , t )  =  U k $. This depends explicitly of the inverse of the amplitude $  f k $, whose evolution satisfies Eq. (10).
.  Kuramoto's model .  Vlasov model . 

Element  oscillator n  Van Kampen eigenmode k 
H_{0}  $ \u2211 n \omega n J n$  $ \u2212 \u222b d 3 p d 3 k  f k  k . p / m \gamma $ 
δH  $ K r \u2211 n J n \u2009 sin ( \Theta n \u2212 \phi n )$  $ \u222b d 3 p d 3 k  U k   f k  \u2009 sin ( \Theta k \u2212 \phi k )$ 
Phase  $ \phi n$  $ \phi k$ 
Action  J_{n}  $ J k =  f k ( p , t ) $ 
Natural frequency  ω_{n}  $ \omega k = \u2212 k . p / ( m \gamma )$ 
Entropy density  $ k B \u2009 ln \u2009 J n$  $ k B \u2009 ln \u2009  W k  = \u2212 k B \u2009 ln \u2009  f k ( p , t ) $ 
Coupling strength  Kr  $  U k  \u2261  R k  /  f k $ 
.  Kuramoto's model .  Vlasov model . 

Element  oscillator n  Van Kampen eigenmode k 
H_{0}  $ \u2211 n \omega n J n$  $ \u2212 \u222b d 3 p d 3 k  f k  k . p / m \gamma $ 
δH  $ K r \u2211 n J n \u2009 sin ( \Theta n \u2212 \phi n )$  $ \u222b d 3 p d 3 k  U k   f k  \u2009 sin ( \Theta k \u2212 \phi k )$ 
Phase  $ \phi n$  $ \phi k$ 
Action  J_{n}  $ J k =  f k ( p , t ) $ 
Natural frequency  ω_{n}  $ \omega k = \u2212 k . p / ( m \gamma )$ 
Entropy density  $ k B \u2009 ln \u2009 J n$  $ k B \u2009 ln \u2009  W k  = \u2212 k B \u2009 ln \u2009  f k ( p , t ) $ 
Coupling strength  Kr  $  U k  \u2261  R k  /  f k $ 
B. The “reinforced” filamentation concept

The magnetic contribution remains dominant over the electrical energy (given the results observed in the numerical simulations).
 The two source terms $ S \u03f5 , c ; k$ and $ S \u03f5 , s ; k$ are negligible, as a first approximation, in particular if the system keeps a certain symmetry according to the component p_{y} of the momentum, due to the presence of the term $ \u2211 j \alpha j ( C j / m \gamma ) F 0 j$, which leads to the electric mean quasineutrality along p_{y}. Note that a strong symmetry breaking in p_{y} leads to the transfer of the momentum in p_{y} and thus to a modification of the conservation of the total quantity $ \u03f5 m , k + \u03f5 \phi , k$, integrated in the configuration space. Equation (17), for a given wave vector k, leads to the conservation of the sum of the magnetic energy spectral density, defined by $ \u03f5 m , k = 1 2 \epsilon 0 c 2  B z , k  2$ and the phasestrophy^{73} flux density $ \u03f5 \phi , k$, satisfying Eq. (18). In the absence of the source term $ S \u03f5 , s ; k$, Eq. (21) reads as$ 1 2 \epsilon 0 c 2  B z , k  2 \u2202 \theta z , k \u2202 t = \u222b d 3 p p x m \gamma  \delta f k  2 2 \u2202 F 0 \u2202 p x ( \u2202 \phi k \u2202 t + k . p m \gamma ).$
If $ \u2202 \theta z , k / \u2202 t \u223c \omega B$, where ω_{B} represents the characteristic frequency of the magnetic mode, i.e., the mode associated with the Weibeltype instability, we recognize two possible cases:

The case of a stationary magnetic mode, with $ \omega B \u223c 0$, induced by a CFItype instability. The synchronization process takes here in a weak coupling regime, without resonance, to a partially synchronized state.

The case of a Weibeltype instability, oblique in nature. Here, the direct cascade associated with oblique modes leads inexorably to a broadened spectrum in wave vectors and to a “reinforced” filamentation regime. The magnetic mode that is induced by the OI corresponds to a propagative mode of nonzero frequency $ \omega B \u2260 0$. A synchronization induced by the resonant term $ \u2202 \phi k / \u2202 t + k . p / m \gamma \u2243 0$ may occur in a neighborhood of some specific mode, which we can call $ k 0$. A transition to a global synchronization becomes possible and the term $ \u2202 \phi k 0 ( p , t ) / \u2202 t + p . k 0 / m \gamma $ tends to zero. Due to the singularity associated with it, this is usually called a secular term since, being the product between a bounded quantity, close to zero but finite and an ever growing one (related to the integrand in the definition (18) of the phasestrophy), it is doomed to diverge. We thus obtain a singularity in the phasestrophy flux $ \u03f5 \phi $. Considering the invariance of $ \u03f5 \phi + \u03f5 m = const$, approximatively valid also in the linear regime (and in the absence of source terms $ S \u03f5 , c ; k$ and $ S \u03f5 , s ; k$, cf. Fig. 6, below), thus implies a decrease of the magnetic energy, which is transformed into kinetic energy.
IV. FREESTREAMING FILAMENTATION VERSUS REINFORCED FILAMENTATION
In principle, the “reinforced” filamentation differs from the usual “freeparticle” filamentation. In the presence of a strong mode coupling, such as the one which can be highlighted by the momentum transfer measured by the phasestrophy introduced in Sec. III, the “reinforced filamentation” is the main mechanism responsible to an entropy production observed in numerical experiments. The entropy increase induced by coarsegraining mechanism, once the size of the filaments reaches the elementary (numerical) grid, is here “enhanced” by the instabilities at play, which may induce a thinning/stretching of the filaments at a rate faster than that, which is given by the “kinematic” filamentation associated with the free streaming of particles. At the same time, selforganization processes in the phasespace mediated by synchronization of oscillatory modes (e.g., in a Kuramotolike scenario) can further affect the filamentation by inducing, e.g., the emergence of filaments at a specific spatial scale. With the notion of “reinforced” or “enhanced” filamentation, we here mean the combination of all these possible effects, in contraposition to the kinematic filamentation given by the free streaming of particles.
In this sense, the variation of entropy related to the reinforced filamentation is also connected to the bifurcation, which involves the phase synchronization. A detailed investigation of the entropy “anomalous” evolution reported in Refs. 1 and 52 goes beyond the scope of this present study and will be reported in a future, dedicated work.
We turn here the attention to the characterization of the features or properties, which give a “reinforced” character to the filamentation of f in the velocity space, especially in connection with a Kuramototype synchronization. In order to analyze in greater detail the coupling between the filamentation of f and the phase synchronization mechanism, and their mutual impact on plasma dynamics as well, we consider the example of a kinematic filamentation. To this end, we consider the free streaming of a Gaussian distribution in the coordinate space, while preserving its Maxwellian distribution. It should be noted that a purely homogeneous initial condition in configuration space is not sufficient to induce a filamentation mechanism of f in the absence of the Lorentz force, so it is necessary here to introduce a spatial inhomogeneity in the distribution. Figure 1 (top frame) shows, at a given instant, the filamentation process of the distribution function in the reduced phase space $ x \u2212 p x$.
The stretching process of the filament leads inexorably to thinner and thinner filaments whose length increases in the phase space. This dynamics is perfectly reproduced by the Vlasov electromagnetic (VLEM) solver,^{78} in which the force term in Vlasov equation was suppressed. Numerical integration shows that the entropy is perfectly conserved in the simulation up to the seventh decimal place. This filament stretching mechanism also corresponds to a phase synchronization of the Van Kampen modes. Figure 1, bottom frame, shows a phase mapping $ \phi ( x , p x , t )$ of the corresponding distribution obtained from a Hilbert transform of the distribution function: we clearly observe locking at a phase $ \phi = \pi / 2$ in regions of phase space where the filament is stretching, while the phase remains close to zero in regions of nearzero density (e.g., in the tail of the distribution, but also between the various filaments in the plasma bulk).
The complex envelop and the fast oscillating phase have been obtained by using a spatial Hilbert transform of the reduced distribution function $ g ( x , p x , t ) = \u222b d p y f ( x , L y / 2 , p x , p y )$ in the $ x \u2212 p x$ phasespace plane, according to the following procedure: the signal g is Fourier transformed in x and the positive k_{x} modes are multiplied by $ + i$, while the negative k_{x} modes are multiplied by –i. The resulting function is Fourier transformed back to obtain g_{Hilbert}, which results in a Hilbert transform. Then, the complex amplitude of the $ g = a ( x , p x , t ) e i \phi ( x , p x , t ) + c . c .$ distribution is given by $ a = g 2 + g Hilbert 2$, while the phase $\phi $ obeys $ tan \u2009 \phi = g Hilbert / g$.
The phaselocking displayed in Fig. 1 appears to be an intrinsic property of the freestreaming system, which completely follows from its kinematics evolution (i.e., it is mediated by the $ v \u2202 f / \u2202 x$ term). In the presence of external forces, which induce modifications of the phasespace distribution function that go beyond the action of the $ v \u2202 f / \u2202 x$ freestreaming term, the phaselocking can give rise to more complex selforganization processes, which may possibly lead to the emergence of localized filaments with some initial characteristic thickness. This will be discussed later in Secs. V and VI, when nonlinear simulations of OI modes will be discussed.
We conclude this section by noting that the behavior of f in this kinematic filamentation process is remarkably similar to that found for the baker transformation, socalled because of its similarity to the kneading and folding operations involved in the preparation of pastry dough. The baker transformation involves a compression in the p_{x} direction followed by a stretching in x. It must be pointed out that the baker transformation has an entropy that is constant for all time and which is equal to the entropy of the initial distribution (a property that results of the Kolmogorov automorphism property^{79} of the baker transformation).
More recently, the possibility of the emergence of a bifurcation in Vlasov systems has been proposed in Ref. 80, when it is driven by an inhomogeneity in the density, where the basic conditions of weak resonance are possible. As we have seen, the spatial inhomogeneity is required also by the kinematic filamentation due to the free streaming of particles. The generalization of this process, in the presence of a nonzero Lorentz force, leads to a profound change in the nature of the filamentation of the distribution function in phase space, involving a reorganization of Van Kampen modes into clusters of synchronized modes.
V. NUMERICAL EXPERIMENTS ON THE NONLINEAR DYNAMICS OF A SINGLE OI MODE
A. Filamentation aspects in PIC and semiLagrangian numerical schemes
In this section, we first recall some general features of the integration scheme of Lagrangian (i.e., PIC) and semiLagrangian solvers, especially in relation to the notion of phasespace filamentation. We will then discuss some numerical results obtained with both types of numerical schemes, in which the filamentation process is highlighted. The simulation cases we consider here are initialized with the same configuration unstable to OI modes, which we will later (Sec. VI) show to lead to a “new” kind of kinetic heating, occurring in a cascadelike turbulent process, when a broad band of unstable modes is excited. However, in order to show and discuss some fundamental features of the nonlinear OI dynamics, which act as “ingredients” of the more complex kinetic heating in a turbulencelike scenario, in the simulation cases considered in this section we will perturb this configuration with a single wavenumber corresponding to the most unstable OI mode compatible with the size of the simulation box.
Filamentation is a real phenomenon in solutions of the Vlasov equation; it is not a spurious side effect of the numerical schemes. The source of filamentation lies in the treatment of the freestreaming term $ \u2202 f / \u2202 t + v \u2202 f / \u2202 x = 0$, whose exact solution in the Fourier space reads as $ f ( k , v , t ) = f ( k , v , 0 ) e ikvt$ and which clearly show that $ f ( k , v , t )$ can oscillate with respect to v at the frequency kt. A first solution for the correct treatment of filamentation of the Vlasov equation was finally given by Knorr and Cheng in Ref. 81, which first introduced the concept of time splitting. This numerical scheme is now known as the semiLagrangian (SL) scheme (see Refs. 78, 82–84), which is based, for each time step, on the two following sequences:

first, transport either f itself (by using the semiLagrangian backward characteristics)

then compute the value of f at the origin of the characteristics by reconstructing the distribution function from the different mesh point values using a cubic spline interpolation technique.
Particle codes involve a Lagrangian formulation:^{85–87} the particleincell (PIC) technique can be regarded as a discretization of the phase space in terms of a superposition of moving elements, each representing a “cloud” of physical particles, usually referred as superparticles. In PIC codes, superparticle trajectories are computed from the electromagnetic fields, prescribed on a fixed grid whose typical size is of the order of the Debye length. At the end of the time step, the charge of each superparticle is redistributed among the neighboring grid points, allowing one to solve the Maxwell equations. Thus, while the SL scheme invariably implies a loss of information, the particle method involves a smoothing of information, which, on the other hand, is in principle numerically conserved during the evolution of plasma. Numerical smoothing techniques, for instance those consisting in redistributing the charge of superparticle among the neighboring points in space, efficiently decrease the individual effects introduced by the grid (which are unwanted since associated with a numerical noise), but they do not change the scaling of these effects. In order to suppress them, we must smooth the electromagnetic field components not only over the interparticle distance but also at the scale of the Debye length. As a drawback, however, we thus begin to modify the possible real collective behavior of the plasma.
Freestreaming filamentation is well resolved in PIC codes and the Lagrangian scheme is not subject to noise due to the filamentation process. This makes the PIC model a very welladapted scheme to study the mechanism of filamentation amplification, but a price must be payed in terms of the corresponding amplified thermal noise.
In the following, both SL (which hereafter we will sometimes refer to as “Vlasov”) and PIC codes are used next to simulate the process of energy transfer. The spatial meshes of both numerical schemes are chosen to be small enough to fully describe both the wave numbers corresponding to the chosen oblique modes and the excitation of largek oblique modes associated with possible direct cascades. The number of grid points of the spatial configuration is chosen in order to obtain a very accurate description of spatial filamentation, if present in the system, i.e., with a sampling of order of $ N x N y \u223c ( L / \Delta x ) d x \u223c 256 2$, L being the typical system length along a given direction and d_{x} the spatial dimension of the problem (here d_{x} = 2).
Both SL Vlasov and PIC codes are initialized with the same sampling in the configuration space. To assume an equivalent description in the momentum space with the two different numerical methods, the ratio between the momentum space sampling used in the (SL) VLEM code^{78} N_{Vlas} to the number of superparticles by cell N_{part} (used in the SMILEI code) is chosen as $ N Vlas / N part = g PIC N p 2 \u223c 1$. This corresponds to a graininess parameter of the order of $ g PIC \u223c 10 \u2212 4$ for the PIC solver SMILEI.^{88} Although both numerical approaches solve the Vlasov equation, the Eulerian or semiLagrangian description corresponds to a statistical approach where the graininess parameter g initially tends toward zero. In the particlesincell approach, the reintroduction of finitesize superparticles leads to a PIC specific graininess parameter g_{PIC} in principle small but finite. In the simulations we have carried out, we have $ N p x N p y = N p 2 = 128 2$ in the momentum space in the SL scheme, which is equivalent to use 128^{2} particles by cell in the SMILEI code (corresponding to a value of the inverse of the graininess parameter of $ g PIC \u2243 6.1 \xd7 10 \u2212 5$). The numerical parameters used in the PIC/Vlasov comparison are detailed in Appendix C. Space variables being normalized in $ d e = c / \omega p$ units, times to the inverse of the plasma frequency $ \omega p \u2212 1$, and energy densities to $ n 0 m c 2$.
B. Vlasov–Maxwell modeling of oblique Weibeltype instability: OI in a counterstreaming configuration
We now discuss some numerical results, obtained with the two complementary numerical approaches, a Lagrangian approach, based on the use of (PIC) SMILEI code and finally the VLEM solver. In both codes, we used the same spatial grid with a mesh size $ \Delta x , \u2009 \Delta y$ of the order of onehundredth of the Debye length and an equivalent number of particles by cell (in the PIC code) or of grid point in the pspace (for the SL VLEM code) as discussed in Sec. III A.
Used code .  Vlem (SL) .  Smilei (PIC) . 

Box length in x  $ L x = 2.51327 d e$  $ L x = 2.51327 d e$ 
Box length in y  $ L y = \pi d e$  $ L y = \pi d e$ 
Dominant wave vector  $ ( k x d e = 5 , k y d e = 2 )$  $ ( k x d e = 5 , k y d e = 2 )$ 
Corresponding mode numbers  $ ( n x = 2 , n y = 1 )$  $ ( n x = 2 , n y = 1 )$ 
Momentum sampling/particle  $ N p x N p y = 128 2$  $ N part by cell = 128 2$ 
Graininess parameter  g_{PIC} = 0  $ g PIC = 6.1 \xd7 10 \u2212 5$ 
Used code .  Vlem (SL) .  Smilei (PIC) . 

Box length in x  $ L x = 2.51327 d e$  $ L x = 2.51327 d e$ 
Box length in y  $ L y = \pi d e$  $ L y = \pi d e$ 
Dominant wave vector  $ ( k x d e = 5 , k y d e = 2 )$  $ ( k x d e = 5 , k y d e = 2 )$ 
Corresponding mode numbers  $ ( n x = 2 , n y = 1 )$  $ ( n x = 2 , n y = 1 )$ 
Momentum sampling/particle  $ N p x N p y = 128 2$  $ N part by cell = 128 2$ 
Graininess parameter  g_{PIC} = 0  $ g PIC = 6.1 \xd7 10 \u2212 5$ 
1. Temporal evolution of the kinetic and electromagnetic energy densities
Guided by the linear analysis of the dispersion relation of OI (see Refs. 9,10, and 51) and summarized in Table II, numerical simulations have been carried out in the weakly relativistic regime to investigate the nonlinear dynamics of this class of solutions up to its nonlinear saturation regime. SL Vlasov–Maxwell simulations, shown in Fig. 2, highlight the different behavior in the temporal dynamics of the propagative OI class with respect to the nonpropagative OI modes. A comparison between the SL VLEM code and the PIC SMILEI code is shown in Fig. 3, which focuses on the temporal evolution of the magnetic and electric energy density components, obtained with the SL VLEM code (left column) and with the SMILEI code (right column). The magnetic energy density shown in Fig. 2 is reproduced to facilitate a direct comparison.
The energy evolution shown in Fig. 3 allows us to identify a fundamental feature of the linear evolution of the OI mode, in which the generation of the magnetic field B_{z}, growing at the expenses of the momentum anisotropy, is coupled to the increase in the electrostatic field, which is typical of the twostream instability (TSI). In particular, the kinetic energy density $ \u03f5 K = \u222b d 3 x V \u222b d 3 p m c 2 ( \gamma \u2212 1 ) f$ (top frame in Fig. 2) is transferred to the magnetic energy density $ \u03f5 m = \u222b d 3 x V 1 2 \epsilon 0 c 2 B z 2$, shown in the bottom frame, and to the total electric energy density defined by $ \u03f5 e = \u222b d 3 x V 1 2 \epsilon 0 ( E x 2 + E y 2 )$ (shown in the middle frame). In the evolution of the nonpropagative OI, (i.e., with $ R e ( \omega ) = 0$), not shown here (see Ref. 1 for more details), the electrostatic field contribution remains negligible. In the propagative OI mode shown in Fig. 2, a first peak of the electrostatic energy density is observed at time $ t \u223c 40 \omega p \u2212 1$ with a maximum amplitude of the order of $ \u03f5 m / n 0 m c 2 \u2243 0.012$, followed by a decrease in the magnetic energy (bottom frame). A second peak of the electron energy density, associated with a reversal of the magnetic energy transfer, appears for $ t \omega p \u2265 100$. As we will see later, in a dedicated section, this process is linked to phasespace filamentation.
A second simulation was run with the same parameters (see Table II), using the SMILEI code. The contribution of the electric energy densities, i.e., the $ \u03f5 e , x$ and $ \u03f5 e , y$ contributions of the total electric energy density ϵ_{e} along the corresponding x and y directions, is shown in Fig. 3. Comparison with Fig. 2 shows that the PIC and Vlasov solvers exhibit an identical behavior of the energy transfer and the dynamics of magnetic energy remains identical in both simulations. Also, the maximum amplitude of the y component of the electric energy obtained from the SMILEI code, shown on bottom right frame in Fig. 3, is very close to that obtained by the VLEM solver (shown on bottomleft frame).
It can be noticed that, taking into account the level of noise specific to the Lagrangian scheme, the level of fluctuations remains important for the E_{y} component of the electric field and this happens despite the use of a large number of particles per cell. However, in spite of the amplified thermal noise, linked to the grain effect reintroduced in the PIC scheme, the growth of E_{y} remains perfectly visible.
2. Contourplots of the magnetic field
During the nonlinear stage of the OI evolution, a transition is observed toward a selforganized state, in which the magnetic field B_{z} displays the characteristic spatial modulation observed in the evolution of the CFI instability. This transition can be identified in Fig. 4 at the beginning of the saturation. Figure 4 displays the dynamics of the (normalized) zcomponent of the magnetic field $ e B z / m \omega p$ in the configuration space, for the OI case (obtained from the VLEM code). Results of the SMILEI code, shown in Fig. 5 for the same plasma state (propagative OI), show an even more complex dynamics. The selforganization begins with a transition toward the CFIlike magnetic configuration, where the plasma adjusts itself in wave vector so to give a metastable state dominated by magnetic trapping (the electric contribution being negligible at that time).
In the nonlinear saturation regime shown in the PIC simulations, however, one can observe an oscillatory process in the dynamics of the magnetic field where the field oscillates, depending on the xspace direction, between two possible asymptotic states. In the CFIlike magnetic configuration, the magnetic field presents a homogeneous structure along the y direction (this type of structure is similar to the one observed at the bottom right, at time $ t \omega p = 40$, in Fig. 4). It appears quite early in the time evolution (until after the first magnetic energy peak in Figs. 2 and 3). The subsequent slight decrease and oscillation of the magnetic field can be interpreted as linked to the existence of a bifurcation between the two oblique modes CFI and OI (note that the OI mode presents an inhomogeneity according to the y direction). This phenomenon has been observed in simulations performed with both PIC and Vlasov codes. In the right column of Fig. 5, we show this oscillation process at two different times chosen on long scales in the saturation regime (simulation performed with SMILEI). The amplitude of the oscillations is however smaller in the VLEM model. This seems to indicate that the Lagrangian integration scheme amplifies the phenomenon because of its larger numerical noise.
3. Time evolution of the phasestrophy flux density
As mentioned in relation to Eq. (18), the phasestrophy $ S \phi = \u222b d 3 x V \u222b d 3 p ( \delta f ) 2 2 \u2202 F 0 \u2202 p x$, introduced in Ref. 73 in the context of the gyrokinetic modeling for studying the turbulence and the transport in tokamaks, is a kind of L^{2} norm representing the phasespace density autocorrelation function linked to the growth of filaments and therefore to the “kinematic” filamentation.
Equations (17) and (21) give an insight on how the energy transfer takes place in the presence of a phase synchronization process by showing a link with the phasestrophy flux density $ \u03f5 \phi , k$ and the source terms $ S \u03f5 , c ; k$ and $ S \u03f5 , s ; k$. When the source terms are both zero, the resulting conditions (19), i.e., $ \u03f5 m + \u03f5 \phi = const$, highlight an energy transfer occurring in the linear stage of the OI, which takes place in combination with the energy conservation $ \u03f5 m + \u03f5 K = const$ (assuming that the electric energy density ϵ_{e} is negligible). This is also verified in the nonlinear regime. The variation of $ \u03f5 \phi $, during the growth of the OI, can be interpreted in terms of momentum transfer,^{74} a process already observed in tokamak turbulence,^{73} of which we take $ \u03f5 \phi , k$ as a measure for the “enhanced” filamentation mechanism associated with the mode k. In Fig. 6, we have shown in the same figure the time evolution of the (total) phasestrophy flux density $ \u03f5 \phi $, of the magnetic energy density ϵ_{m,} and of their sum $ \u03f5 \phi + \u03f5 m$, in the linear phase of the OI, i.e., for a time $ t \omega p \u2264 35$. We observe that the quantity $ \u03f5 \phi + \u03f5 m$ is very well conserved in the linear regime of the OI, at least as long as the level of synchronization remains low. This is consistent with the results already discussed in Refs. 1 and 74. As synchronization takes place, at the end of the linear stage, the term $ \u222b d 3 p p x m \gamma  \delta f k  2 2 \u2202 F 0 \u2202 p x ( \u2202 \phi k \u2202 t + k . p m \gamma )$ in the righthand member of Eq. (22), although finite, leads to a growth in the flux density of phasestrophy. In particular, an abrupt change of several order of magnitudes is measured when the quantity $ \u2202 \phi k \u2202 t + k . p m \gamma $ tends toward zero. As already discussed in Ref. 1, this occurs when a global synchronization takes place. Although this happens in the advanced nonlinear stage of the instability (cf. Fig. 7), in which Eq. (19) is formally not valid anymore, the growth of the enstrophy evaluated according to Eq. (18) is still a reliable proxy for the momentum exchanged during the synchronization process. Indeed, Eq. (18), is defined regardless of the ordering of $  \delta f k  2 / ( \u2202 F 0 / \u2202 p x )$. One may also look at Fig. 2 of Ref. 1, where a striking difference in the phasestrophy evolution is visible, between the simulation case considered here (i.e., Fig. 7, here), and another simulation case in which no synchronization takes place (i.e., in which only nonpropagative oblique modes are excited).
Such conservation is thus verified only in the absence of synchronization [when the phase synchronization takes place, the source term $ S \u03f5 , c ; k$, which contains some terms like $ cos \u2009 ( \theta y , k \u2212 \phi k )$ in Eq. (B10), becomes nonnull]. In that case, a transfer of energy is possible leading to an increase in the phasestrophy or the phasestrophy flux density (in absolute value). The time evolution of the phasestrophy flux density $  \u03f5 \phi $ is shown in Fig. 6. The phasestrophy flow is not a Casimir invariant in the strict sense; however, it evolves little in the absence of a synchronization mechanism. Its divergence is instead associated with the “bifurcation” process that leads to a global synchronization in the Kuramoto model.
In Fig. 7, we thus observe two major events. The first one, at $ t \omega p \u223c 50$, is related to the transition to the CFIlike state, with $ \u03f5 K \u2192 \u03f5 m$. A second event takes place at a later time, for $ t \omega p \u223c 150 \u2212 250$, and is linked to the reversal of the energy transfer, i.e., $ \u03f5 m \u2192 \u03f5 K + \u03f5 e , y$: in the presence of an OI instability, the phases of the Van Kampen modes are locked (since the phase $ \Theta k \u2212 \phi k$ tend to zero, and the source terms $ S \u03f5 , c ; k$ and $ S \u03f5 , s ; k$ are modified), which leads to a growth of the quantity $  \u03f5 \phi , k $.
The first growth in the phasestrophy density flux around time $ t \omega p \u223c 40 \u2212 100$ results directly from the OI and from the associated conversion of the kinetic energy of the particle beams into magnetic energy, the electrical energy being negligible in this case. The second growth phase, occurring at time $ \omega p \u223c 150 \u2212 250$, results from the reversal in the energy transfer process, where the magnetic energy is converted into kinetic energy due to the enhanced filamentation. This process can be observed directly from Eq. (22). During the phase synchronization, the term on the righthand side of Eq. (22), i.e., $ ( \u2202 \phi k \u2202 t + k . p m \gamma )$ in the integrand, tends toward zero, due to the transition associated with the global phase synchronization mechanism mentioned above in relation to Fig. 7. In particular, given that the lefthand term $ \u03f5 m \u2202 \theta z , k \u2202 t$ of Eq. (22) remains finite, this transition is accompanied by a divergence of the term $ p x m \gamma  \delta f k  2 2 \u2202 F 0 \u2202 p x$ and therefore of the flux of the phasestrophy density. At a microscopic level, the reversibility of the energy transfer is compatible with the time reversibility of the Vlasov equation, which is responsible, e.g., both of the reversibility of the kinematic filamentation and of the phenomenon of the plasma echoes.
C. Role of phasespace filamentation
In order to provide a more detailed view of the physical filamentation mechanism involved, we have represented in Fig. 8 a plot of the distribution function in the $ x \u2212 p x$ phase space, evaluated at $ y = L y / 2$ and averaged along p_{y} (the corresponding simulation is carried out with the VLEM solver).
The focus here is on the phasespace dynamics of the distribution function f at the beginning of the OI. The filamentation process takes place since the beginning of the first peak in the phasestrophy growth, which occurs around the times $ t \omega p = 45 \u2212 50$. At these times, the formation of two filaments in the vicinity of the X points is observed. These structures, which appear in the form of a very narrow filament, stretching progressively in the phase space, are characteristic of the process of the “reinforced” filamentation of f. The adjective “reinforced” is used in contraposition to the “kinematic” filamentation due to the free convection of the distribution function, so to stress the fact that here the instabilities at play enhance the rate at which filamentation occurs in the phase space. These filamentary structure can also split in two, as observed, at time $ t \omega p = 50$, in the vicinity of point X (top right), and can lead to the appearance of two or more “tilting stripes,” in the phase synchronization processes. Note that these tilting stripes are related to the enhanced filamentation process and differ from the analogous structures (Sec. IV), which are usually observed during the (linear) plasma dynamics of both the Vlasov and Kuramoto models (cf. Sec. IV).
These filamented structures have been already observed in the continuum limit of the Kuramoto's model in Ref. 89, in the case of an infinite population of oscillators. Such pairwise filamentary structures (often characterized by a thin critical layer in phase space) are a consequence of the phase synchronization. This synchronization pattern (which is also shown in Fig. 8 on top right, close to the X points) emerges also during the bifurcation toward a globally synchronized state in the case of the discrete Kuramoto model from a continuum spectrum.
Figure 9 shows the formation of a thin filament obtained in a simulation performed with the SMILEI code: the emergence of this thin filamented structure is visible from time $ t \omega p \u223c 30 \u2212 40$, and it follows the “stretching” process of the vortex structure, located in the central region.
This stretching of the filament leads to thinner filaments whose length increases while their thickness decreases (following the mass conservation imposed by the Liouville theorem) while tending to a zero thickness (which is comparable to a Dirac type distribution). This dynamics is perfectly reproduced by the SMILEI code and is shown in the left frames of Fig. 9. In the top right frame of the same figure, at time $ t \omega p = 40$, the filament located in the vicinity of an “Xtype” point (i.e., a bifurcation point in the phase space) disappears at $ t \omega p = 44$. The “Otype” point where there is a phasespace vortex (e.g., at time $ t \omega p = 42$) is also stretched and deformed in time, but it literally splits into two parts, two curved filaments which mutually move away from each other and thus form a kind of boundary layer. These filaments keep on extending and getting thinner over time. The Xtype points appear to be unstable with respect to the phasespace dynamics and quickly lead to the emergence of further filaments, which are perfectly visible in the bottom right of Fig. 9, at $ t \omega p = 52$.
VI. COLLISIONLESS PLASMA HEATING IN A TURBULENCELIKE CASCADE INITIALIZED WITH A COUNTERSTREAMING BEAM CONFIGURATION
We have shown that the global energy exchange between kinetic and magnetic energy, derived from the Vlasov equation, is mediated by the coupling mechanism between the phase synchronization of Van Kampen modes (whose dynamics is described by basic oscillators in the Fourier spectrum modeled in a meanfield Kuramoto framework) and a “reinforced” filamentation process (which differs from the usual kinematic filamentation process) of the Vlasov equation. As previously mentioned in Sec. III, the $ \u2212 k . p / m \gamma $ free advection term in the time derivative of the phase in Eq. (11) can trigger an alternative “channel” of the conversion between kinetic and electromagnetic energy, leading, in the regime of “reinforced filamentation” discussed above, to the reversal of the energy transfer. The role of phasespace filamentation, and in particular, of a spatially sheared flow, in the anisotropic conversion between order kinetic energy of the plasma to its internal (thermal) energy had been already put in evidence in Refs. 38, 90, and 91.
This process provides a new perspective on the collisionless dissipation mechanisms and on the collisionless plasma cascade. In this section, we discuss a last numerical experiment, showing that the collisionless plasma heating induced by the OI nonlinear dynamics is an efficient mechanism of energy transfer in a turbulentlike scenario in which several modes are excited. We will also show that the nonlinear dynamics associated with this process may lead to a global selforganization of the system, due to a bifurcationtype mechanism, which can be modeled in the meanfield (Hamiltoniantype) Kuramoto approach.
We will show indeed that a global reorganization takes place in which collisionless wave–particle interactions transfer the “free” ordered kinetic energy stored in the highenergy beams into disordered kinetic energy, by increasing the internal energy of the plasma and thus heating it. The first step of the process consists in the generation and growth of an electromagnetic field excited by the Weibeltype instability (i.e., OI modes) at the expenses of the initial ordered kinetic energy of the beams. This magnetic energy acts as a temporary “reservoir” of potential energy, which is later reconverted into disordered kinetic energy, i.e., internal energy of the plasma. In the turbulencelike scenario, this conversion from magnetic to internal energy occurs if a bifurcation in the Vlasov–Maxwell system takes place when an initial perturbation encompassing a broad spectrum of wave numbers is applied.
A. Numerical experiments on the nonlinear dynamics of a broad spectrum of OI modes
.  Case OI with broad spectrum . 

Used code  Vlem (SL) 
Box length in x  $ L x = 2 \pi d e$ 
Box length in y  $ L x = 2 \pi d e$ 
Wave numbers excited  $ n x , n y \u2208 [ 1 , 15 ]$ 
Momentum sampling/ particle  $ N p x N p y = 128 2$ 
Magnetic perturbation amplitude  $ e B 0 / m \omega p = 5 \xd7 10 \u2212 5$ 
Corresponding figures  Figs. 10–12 
.  Case OI with broad spectrum . 

Used code  Vlem (SL) 
Box length in x  $ L x = 2 \pi d e$ 
Box length in y  $ L x = 2 \pi d e$ 
Wave numbers excited  $ n x , n y \u2208 [ 1 , 15 ]$ 
Momentum sampling/ particle  $ N p x N p y = 128 2$ 
Magnetic perturbation amplitude  $ e B 0 / m \omega p = 5 \xd7 10 \u2212 5$ 
Corresponding figures  Figs. 10–12 
Figure 10 shows the temporal evolution of the different energy densities: kinetic ϵ_{K} (top), electric ϵ_{e} (middle), and finally magnetic ϵ_{m} (bottom). The physical parameters are identical to those used in the previous simulations corresponding to the excitation of a single OI mode, discussed in Sec. V. Differently from those cases, here a large spectrum is initially excited.
The use of an extended spectrum in wave vectors allows to excite a wide range of modes. These include several modes in both the nonpropagative and propagative branches of OI modes, but also the standard TSI and CFI instabilities. This initialization, which is likely more realistic than the cases considered for theoretical simplicity in Sec. V, also allows for a direct cascade process. In Ref. 52, we have indeed shown that the propagative branch of the OI modes is characterized by a saturation mechanism (different from the magnetic trapping typical of the nonpropagative OI branch occurring at $ k d e \u2273 1$), which is dominated by electrostatic trapping and by the nonlinear excitation of increasingly larger wavenumbers, up to values $ k d e \u226b 1$.
The time evolution of the energy densities displayed in Fig. 10 shows the emergence of four different phases of this dynamics:

step (I): a first stage is characterized by an exact conservation of the energy contributions, corresponding to the initial equilibrium in the time interval $ [ 0 , 17 \omega p \u2212 1 ]$,

step (II): a second stage is characterized by the growth of the electromagnetic energy in both the electric (middle frame) and magnetic components (lower frame). This is associated with a corresponding decrease of the kinetic energy (upper panel). This second phase of the instability gives rise to an energy transfer of kind $ \u03f5 K \u2192 \u03f5 m + \u03f5 e$, which saturates at time $ t \omega p \u223c 35$ in Fig. 10.

step (III): the reversal of the energy transfer (i.e., $ \u03f5 m \u2192 \u03f5 K ; \u03f5 e \u226a \u03f5 K , \u03f5 m$) takes place in the third phase of the plasma dynamics, in the time interval $ [ 40 , 100 ] \omega p \u2212 1$ and is characterized by a slow decrease of both the electric ϵ_{e} and magnetic ϵ_{m} energy densities.

step (IV): the last phase corresponds to the saturation regime where the kinetic energy has reached its maximum level in amplitude, close to its initial level and where the magnetic energy exhibits a lowfrequency oscillation at a low level of amplitude.
Two main stages can be recognized in the four steps identified above, in relation to the crucial role played by the OI modes: a first, main stage (constituted by steps (I) and (II)) is characterized by the growth of the most unstable OI modes. This is associated with an energy transfer from the kinetic ϵ_{K} to electromagnetic $ \u03f5 e + \u03f5 m$ energy density, till saturation, observed at time $ t \omega p \u223c 40$. A second stage (constituted by steps III and IV) initiates with the saturation of the most unstable OI modes. This gives rise, due to an “augmented” filamentation process, to the reversal of the energy transfer from magnetic to kinetic, leading to a strong plasma heating process. Mass and energy are well conserved during the entire simulation, their relative variations being of the order of $ 4.00 \xd7 10 \u2212 4$ and $ 1.20 \xd7 10 \u2212 4$, respectively. In this case in which a large spectrum of oblique modes is excited, the amount of kinetic energy transferred to magnetic energy remains quite small.
Figure 11 illustrates the dynamics of the normalized magnetic field component, $ e B z / m \omega p$ at different times (left column,) together with the plot of the electron density, normalized to n_{0}, at the same instants (right column). In the strong nonlinear regime, where the phase synchronization is expected to be quasiglobal (i.e., leading to a resonant filamentation process, forced by the synchronization of Van Kampen modes), the excitation of a wide spectrum leads also to the excitation (and thus “linear” growth) of the nonpropagative OI mode (with wavenumber n_{x} = 2 and n_{y} = 1). This is shown in Fig. 11, at time $ t \omega p = 34$, i.e., at a time close to the beginning of the saturation phase (bottomleft frame). Even if, during the whole process the electric field component remains weak (cf. the relative amount of electrostatic energy in the center frame of Fig. 10, its impact on saturation of the OI is not negligible. In the representation of the electron density in Fig. 11, we observe the emergence of thin filaments associated with a transition from a dominance of lowk OI modes to a dominance of largek OI modes, all corresponding to the propagative branch of OI, as it was already observed in Ref. 52. Such a transition is linked to the amplification of fluctuations of f due to a coupling with a process of filamentation that takes place in both the velocity and the coordinate spaces.
In Fig. 12, we show the contours of the distribution function in the pspace. The distribution function is concentrated in space, i.e., located at $ x = L x / 2$ and $ y = L y / 2$. The plasma heating scenario, in which the initial ordered energy associated with the beams is finally converted into thermal (random) energy of the electron, is visible in the broadening of the initially biMaxwellian distribution function over time. The top left panel represents the electron beams behavior, at time $ t \omega p = 16$, and exhibits a symmetry breaking in electron beams, showing the appearance of an inhomogeneous structure. Closer inspection of different patterns in $ x \u2212 p x$ or $ y \u2212 p y$ subphase spaces (not shown here) reveals, as expected, that the electron dynamics exhibits finer substructures consisting of thin filaments at the electrons scales down to the (electron) inertial lengths.
VII. CONCLUSION
The phasespace filamentation of the distribution function and the synchronization of Van Kampen modes are two intrinsic features of the Vlasov–Maxwell system, which we have shown to play a crucial role in the conversion of magnetic energy into internal energy of a collisionless plasma. In this work, we have used two complementary numerical schemes, namely, the PIC and semiLagrangian models, to show with numerical experiments how these features intervene in the nonlinear dynamics of oblique Weibeltype instabilities (Secs. V and VI). In this regard, we have also pointed out some specific features of the concept of “reinforced” or “enhanced” (phasespace) filamentation, which we have introduced in contraposition to the notion of “kinematic” filamentation of a freestreaming bunch of particles (Sec. IV). We have also shown (Sec. III and Appendix B) how the combined effect of this “reinforced filamentation” with synchronization aspects of Van Kampen modes described in a Kuramototype representation of the Fourier modes of the Vlasov plasma can be diagnosed with a detailed examination of the time history of an integral quantity, the “phasestrophy,” which measures the momentum exchanged between beams of cold electrons. This quantity also contributes to a quantity which, during the linear stage of the instability, is a constant of motion, the sum of the phasestrophy and of the magnetic energy density $ \u03f5 \phi + \u03f5 m = const$.
We have applied these tools to better characterize (Sec. VI) a new kinetic heating mechanism, first identified in Ref. 1, which can be relevant to the energy transfer in a turbulentlike energy cascade. This process takes place during the nonlinear dynamics of oblique Weibeltype modes, thanks to the role played by their propagative branch, which was identified in Refs. 51 and 52. This branch displays a saturation mechanism which—differently from other Weibeltype modes—is dominated by electrostatic trapping rather than by magnetic trapping, and which naturally leads via nonlinear coupling to the excitation of higher wave number OI modes. This process resembles a classical turbulence cascade. A remarkable feature of this kinetic heating process is the reversal of the energy transfer from the kinetic to the magnetic energy components that is typical of Weibeltype modes. Interestingly, the nontrivial conversion from magnetic to kinetic energy is sometimes also observed in nonideal plasmas, as for instance in magnetic reconnection. However, in the geometry here considered, magnetic reconnection cannot take place [the suppression of the electron dynamics along the zcomponent forbids the generation of magnetic fields in the (x, y) plane]. Thus, the kinetic heating related to the reinforced filamentation studied here, intrinsic to the Vlasov equation, is more general than and remains independent of the magnetic reconnection.
However, if one looks at the role of phasemixing induced by phasespace filamentation, which several previous works have suggested to be at the basis of kinetic magnetic reconnection itself (see Refs. 92–96), this kinetic heating process seems to be more fundamental than magnetic reconnection itself. In particular, two interesting elements deserve to be investigated in future works: if and when magnetic reconnection in a kinetic regime can originate from this kinetic heating process; and which elements of the kinetic heating mechanism, here identified to occur in the nonlinear evolution of a counterstreaming beam configuration, may be generally relevant to magnetic reconnection and to kinetic turbulence at the electron scales.
ACKNOWLEDGMENTS
The authors are indebted to the IDRIS computational center, Orsay, France, for computer time allocation on their computers and access to the cluster EXPLOR of the University of Lorraine (project “Mechanisms of energy conversion in collisionless plasmas” 2019M4XXX0978). This work was granted to the HPC resources (Grants A0130507290 and A0150507290) made by GENCI (Grand Equipement National de Calcul Intensif). Partial fundings are from the French Federation for Magnetic Fusion Studies (FRFCM) through the FRFCM AAP2023 Evolution of current sheets in lowcollision plasmas. During the completion of this work, one of the authors (H. Betar) has received financial support from the AIM4EP project (ANR21CE300018), funded by the French National Research Agency (ANR). 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
A. Ghizzo: Software (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal). D. Del Sarto: Conceptualization (equal); Methodology (equal); Validation (equal); Writing – review & editing (equal). H. Betar: Methodology (equal); Resources (equal); Software (equal); Visualization (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request. Different codes may be available.