# Collisionless heating in Vlasov plasma and turbulence-driven filamentation aspects

A. Ghizzo, D. Del Sarto, H. Betar

► To cite this version:

A. Ghizzo, D. Del Sarto, H. Betar. Collisionless heating in Vlasov plasma and turbulence-driven filamentation aspects. Physics of Plasmas, 2024, 31 (7), pp.072109. 10.1063/5.0205253. hal-04655073

HAL Id: hal-04655073

https://hal.science/hal-04655073v1

Submitted on 21 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
# Collisionless heating in Vlasov plasma and turbulence-driven filamentation aspects

A. Ghizzo

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

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

D. Del Sarto

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

H. Betar

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

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

Check for updates

Phys. Plasmas 31, 072109 (2024)

https://doi.org/10.1063/5.0205253

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

View
Online

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

Export
Citation

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

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

<figure><img src="image_11.png" /></figure>
# Collisionless heating in Vlasov plasma and turbulence-driven filamentation aspects

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

Submitted: 25 February 2024 · Accepted: 30 June 2024 ·

Published Online: 18 July 2024

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

View Online

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

Export Citation

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

CrossMark

A. Ghizzo,<fnref n="1" />,ª ID D. Del Sarto,<fnref n="1" /> ID and H. Betar<fnref n="2" /> ID

## AFFILIATIONS

## ABSTRACT

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 current-driven 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 self-organization 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 Weibel-type instabilities.

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

## I. INTRODUCTION

In a configuration characterized by two counter-propagating electron beams, Weibel-type instabilities can couple with electrostatic unstable modes by generating the so-called oblique instabilities (OIs). Recently, in Ref. 1, we have shown via kinetic simulations—performed with both semi-Lagrangian (SL) and particle-in-cell (PIC) Vlasov-Maxwell codes—how the spatial filamentation induced by these modes is capable of affecting, both linearly and nonlinearly, the phase-space filamentation of the distribution function by synchronizing the phases of the Van Kampen modes. In this work, we further discuss the role played by phase-space 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 large-amplitude 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 2D2V 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 Weibel-type instabilities. Previous linear studies of Weibel-type instabilities$^{2–14}$ have identified different classes of electron beam-plasma instabilities: the current filamentation instability$^{15–19}$ (CFI), the Weibel instability (WI), the electrostatic two-stream 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 gamma-ray 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_⊥$ is larger than the parallel temperature ($T_‖$ here), where the symbols $⊥$ and $‖$ 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

<footnote n="1">¹IJL-UMR 7168, Campus ARTEM, 2 allee Andre Guinier, BP 50840, 54011 Nancy Cedex, France</footnote>

<footnote n="2">²Aix-Marseille Universite, CNRS, Centrale Marseille, M2P2, UMR 7340 Marseille, France and Renaissance Fusion, 38600 Fontaine, France</footnote>

<footnote>ª) Author to whom correspondence should be addressed: alain.ghizzo@univ-lorraine.fr</footnote>
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.²⁵⁻³³ 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 frozen-in 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.³⁴⁻⁴⁵

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,⁴⁶⁻⁴⁹ so as the role played by phase-mixing 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 two-dimensional (2D) configuration consisting of two counter-propagating 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 self-consistent 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 non-propagative branch,” exhibits the generation and the growth of a quasi-static 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 non-propagative 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 direct-like cascade to smaller scales. Here, we will focus on the physical mechanism that leads to the reversion of the energy transfer, associated with this direct-like 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.³³ 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 entrophy 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. Phase-space 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,¹ we investigate here in greater detail the features of the phase-space 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 free-streaming and reinforced filamentation. Section V presents the numerical particle-in-cell and semi-Lagrangian (Vlasov) schemes and discusses the results of the study of two classes of oblique Weibel-type modes, propagating and non-propagating. 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 PHASE-SPACE FILAMENTATION

### A. The physical model

Assuming a linearly polarized electromagnetic wave in the form of two electric and magnetic field components $E = (E_x, E_y, 0)$ and $B = (0, 0, B_z)$, our model uses a two-dimensional (2D) configuration where the “simplified” electron Vlasov equation for a phase-space distribution function $f = f(x, y, p_x, p_y, t)$ reads

$$ \frac{\partial f}{\partial t} + \frac{p_x}{m\gamma} \frac{\partial f}{\partial x} + \frac{p_y}{m\gamma} \frac{\partial f}{\partial y} + e \left( E_x + \frac{p_y B_z}{m\gamma} \right) \frac{\partial f}{\partial p_x} + e \left( E_y - \frac{p_x B_z}{m\gamma} \right) \frac{\partial f}{\partial p_y} = 0. \quad (1) $$

The corresponding electromagnetic field ($E, B$) obeys the Maxwell equations
$$ \frac{\partial E_x}{\partial t} = c^2 \frac{\partial B_z}{\partial y} - \frac{J_x}{\epsilon_0}, \qquad (2) $$

$$ \frac{\partial E_y}{\partial t} = -c^2 \frac{\partial B_z}{\partial x} - \frac{J_y}{\epsilon_0}, \qquad (3) $$

$$ \frac{\partial B_z}{\partial t} = \frac{\partial E_x}{\partial y} - \frac{\partial E_y}{\partial x}, \qquad (4) $$

together with the Poisson law ($divE = \rho/\epsilon_0$), where the electron current density $J = (J_x, J_y, 0)$ and the charge density $\rho$ are defined by

$$ J = e \int \frac{p}{m\gamma} fd^3p \quad \text{and} \quad \rho = e \int fd^3p - en_0. \qquad (5) $$

Here, $n_0$, $e<0$, and $m$ are the fixed background ion density, the elementary charge, and rest mass of electron, while $\gamma = \sqrt{1+p^2/m^2c^2}$ represents the Lorentz factor. $J$ and $\rho$ verify together the continuity equation

$$ \frac{\partial \rho}{\partial t} + \text{div}J = 0 \qquad (6) $$

that expresses the local charge conservation.

## 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 $\mathcal{G}(f)$, $\mathcal{G}$ being an arbitrary function of $f$. This conservation leads to the so-called Casimir invariants $\mathcal{C}(f) = \int d^3x d^3p \mathcal{G}(f) = \text{const}$. Their conservation follows from the invariance of the Vlasov equation under time reversal and from the fact that $\mathcal{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 entropy-preserving 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 phase-space 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 Weibel-type instabilities can induce small-scale, large-amplitude 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 phaseostrophy flux within the Kuramoto description.

# III. SYNCHRONIZATION VERSUS FILAMENTATION IN THE VLASOV FORMALISM

## A. Synchronization aspects in the VM system

Expressing the Vlasov equation (1) in the following form:

$$ \frac{\partial f}{\partial t} + \frac{p}{m\gamma} \cdot \frac{\partial f}{\partial \mathbf{x}} = R(\mathbf{x}, \mathbf{p}, t), \qquad (7) $$

where

$$ R(\mathbf{x}, \mathbf{p}, t) = -e \left( E + \frac{\mathbf{p} \times \mathbf{B}}{m\gamma} \right) \cdot \frac{\partial f}{\partial \mathbf{p}}. \qquad (8) $$

The Fourier transform of Eq. (7) reads

$$ \frac{\partial f_k}{\partial t} = - \frac{i \mathbf{k} \cdot \mathbf{p}}{m\gamma} f_k + R_k(\mathbf{p}, t). \qquad (9) $$

Thus, the Fourier component $f_k(\mathbf{p}, t)$ of the distribution function $f$ can be described, from an oscillator point of view, by its real amplitude $|f_k(\mathbf{p}, t)|$ and its phase $\varphi_k(\mathbf{p}, t)$. By separating real and imaginary parts, Eq. (9) leads to the set of differential equations:

$$ \frac{\partial |f_k|}{\partial t} = |R_k| \cos (\Theta_k - \varphi_k), \qquad (10) $$

$$ \frac{\partial \varphi_k}{\partial t} = - \frac{\mathbf{k} \cdot \mathbf{p}}{m\gamma} + \frac{|R_k|}{|f_k|} \sin (\Theta_k - \varphi_k), \qquad (11) $$

where $|R_k|$ and $\Theta_k = \Theta_k(\mathbf{p}, t)$ define the amplitude and the corresponding phase of the quantity $R_k = |R_k| e^{i\Theta_k(\mathbf{p}, t)}$. Equation (11) describes the Kuramoto dynamics (see Refs. 64–67), in which the angles are coupled through the sine function. The synchronization model was used in Refs. 68 and 69 to study the Landau damping in the Vlasov-Poisson model. Extension to the gyrokinetic modeling was proposed in Refs. 70 and 71 for the study of the ion-temperature-gradient instability in tokamaks.

In Eq. (11), the first term in the right-hand side of the equation corresponds to the intrinsic filamentation process of $f$ in the momentum space (coming from the free-streaming term), which is related to the Van Kampen modes.53 This term plays the role of the natural frequency $\omega_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(\mathbf{p}, t)|/|f_k(\mathbf{p}, t)| = |U_k|$. This depends explicitly of the inverse of the amplitude $|f_k|$, whose evolution satisfies Eq. (10).

The set of Eqs. (10) and (11) for the Vlasov–Kuramoto approach is equivalent to the set of equations in the Hamilton–Jacobi formalism given in (A2) and (A3). Table I summarizes the correspondence between the two models. By rewriting the Hamiltonian in the form

$$ \mathcal{H} = - \int d^3 p d^3 k \frac{\mathbf{k} \cdot \mathbf{p}}{m^\gamma} J_k + \int d^3 p d^3 k |U_k| J_k \sin (\Theta_k - \varphi_k), $$

and by introducing an action $J_k = |f_k|$, the set of Eqs. (10) and (11) reads

$$ \dot{J}_k = - \frac{\partial \mathcal{H}}{\partial \varphi_k} = |U_k| J_k \cos (\Theta_k - \varphi_k), \quad (12) $$

$$ \dot{\varphi}_k = \frac{\partial \mathcal{H}}{\partial J_k} = - \frac{\mathbf{k} \cdot \mathbf{p}}{m^\gamma} + |U_k| \sin (\Theta_k - \varphi_k). \quad (13) $$

Using the Maupertuis principle, we clearly see that it is the quantity $|f_k(\mathbf{p}, t)|$ that plays the role of the action $J_n$ in Eq. (A2) and that it is possible to define, for each Van Kampen mode of wave vector $\mathbf{k}$, an entropy density $S_k$ written in terms of the logarithm of the action density as

$$ S_k = -k_B \ln |f_k(\mathbf{p}, t)| = k_B \ln W_k. \quad (14) $$

Here, $W_k = 1/|f_k(\mathbf{p}, t)|$ denotes a “number” of configurations or micro-states characterized by the wavevector $\mathbf{k}$ in the $[\mathbf{k}, \mathbf{k} + d\mathbf{k}]$ interval. Thus, in Eq. (14), the entropy density $S_k$ counts the number of effective Fourier “modes” in an interval $[\mathbf{k}, \mathbf{k} + d^3\mathbf{k}]$, exactly as the thermodynamical (Boltzmann’s) entropy $S_B = k_B \ln W$ counts the number of microscopic configurations $W$ corresponding to a macroscopic state. This provides an interpretation of synchronization in the strong coupling regime in terms of entropy: the “network” of Van Kampen oscillators tends to synchronize by “maximizing” their entropy production, a feature already observed in the Hamiltonian framework of the Kuramoto model, as previously mentioned in Ref. 72.

## B. The “reinforced” filamentation concept

We can consider without loss of generality a simplified system with wavevector $\mathbf{k} = (k_x, k_y, 0)$. In order to study the influence of CFI on the transition toward spontaneous synchronization, it is illuminating to perform a linear analysis of the Vlasov–Maxwell set of Eqs. (1)–(5): we expand the distribution functions $f$, Fourier transformed with respect to $\mathbf{x}$, through a linearization of $f$ around the equilibrium $F_0(\mathbf{p}) = n_0 \sum_j \alpha_j F_{0j}$ as

$$ f = n_0 \sum_j \alpha_j F_{0j} + \delta f_k(\mathbf{p}, t) e^{i(k_x x + k_y y)}. \quad (15) $$

$F_{0j}$ describes here the initial counterstreaming electron beam distribution. We have considered the case of a weakly relativistic distribution function of two-shifted Maxwellians of velocity $C_2 \simeq 0.90mc = -C_1$ and of temperature $T_1 = T_2 = 6$ keV, the distribution being isotropic. The use of oblique propagative modes allows us to study the filamentation and synchronization processes from an innovative point of view, because the dynamics of these modes occurs at the interface between the concept of instability (CFI in a simplified configuration and the cascade process) and the concept of fully developed turbulence, which takes place during the wave–particle interaction processes.

The Maxwellian distribution of each stream depends on the drift momentum $C_j$ according to

$$ F_{0j} = A \exp\left(-\frac{p_x^2}{2m^2 c^2 \beta_x^2}\right) \exp\left(-\frac{(p_y + C_j)^2}{2m^2 c^2 \beta_y^2}\right), \quad (16) $$

where we have assumed the current neutrality condition $n_0 \sum_{j=1,2} \alpha_j F_{0j} (C_j/m^\gamma) = 0$. Moreover, we only consider a two-counterstreaming case without perpendicular/longitudinal temperature anisotropy, i.e., we assume here the condition $\beta_x = \beta_y = \beta = v_{th}/c$; $v_{th}$ being the thermal velocity. With the definitions of $\delta f_k = |\delta f_k| e^{i\varphi_k}$ for the perturbed distribution function, $E_{x,k} = |E_{x,k}| e^{i\theta_{x,k}}$, $E_{y,k} = |E_{y,k}| e^{i\theta_{y,k}}$ for the electric field components, and $B_{z,k} = |B_{z,k}| e^{i\theta_{z,k}}$ for the magnetic field component, a little algebra leads to

$$ \frac{\partial}{\partial t} (\epsilon_{em,k} + \epsilon_{\phi,k}) = S_{c,c,k}. \quad (17) $$

<table><caption>**TABLE I.** Correspondence table between (the Hamiltonian) Kuramoto’s model and the Vlasov model. We have introduced the notation $R(x, p, t) = -e(E + p \times B/(m^\gamma)) \cdot \nabla_p f$, where $f = f(x, p, t)$ is the electron distribution function and $f_k(\mathbf{p}, t) = |f_k| e^{i\varphi_k(p, t)}$ is the Fourier transform of $f$. A corresponding Hamiltonian can be built in the form $H = H_0 + \delta H$, where $H_0 = - \int d^3 p d^3 k |f_k| k \cdot \mathbf{p} / m^\gamma$; $\delta H = \int d^3 p d^3 k |U_k| J_k \sin(\Theta_k - \varphi_k)$, where $|U_k| = |R_k|/|f_k|$ is independent of the quantity $|f_k|$.</caption><thead><tr><th></th><th>Kuramoto's model</th><th>Vlasov model</th></tr></thead><tbody><tr><td>Element</td><td>oscillator n</td><td>Van Kampen eigenmode k</td></tr><tr><td>H<sub>0</sub></td><td>&sum;<sub>n</sub> &omega;<sub>n</sub>J<sub>n</sub></td><td>- &int; d<sup>3</sup>p d<sup>3</sup>k |f<sub>k</sub>| k &sdot; p / m<sup>&gamma;</sup></td></tr><tr><td>&delta;H</td><td>Kr &sum;<sub>n</sub> J<sub>n</sub> sin(&Theta;<sub>n</sub> - &varphi;<sub>n</sub>)</td><td>&int; d<sup>3</sup>p d<sup>3</sup>k |U<sub>k</sub>| |f<sub>k</sub>| sin(&Theta;<sub>k</sub> - &varphi;<sub>k</sub>)</td></tr><tr><td>Phase</td><td>&varphi;<sub>n</sub></td><td>&varphi;<sub>k</sub></td></tr><tr><td>Action</td><td>J<sub>n</sub></td><td>J<sub>k</sub> = |f<sub>k</sub>(p, t)|</td></tr><tr><td>Natural frequency</td><td>&omega;<sub>n</sub></td><td>&omega;<sub>k</sub> = -k &sdot; p / (m<sup>&gamma;</sup>)</td></tr><tr><td>Entropy density</td><td>k<sub>B</sub> ln J<sub>n</sub></td><td>k<sub>B</sub> ln |W<sub>k</sub>| = -k<sub>B</sub> ln |f<sub>k</sub>(p, t)|</td></tr><tr><td>Coupling strength</td><td>Kr</td><td>|U<sub>k</sub>| = |R<sub>k</sub>| / |f<sub>k</sub>|</td></tr></tbody></table>
Here, $\epsilon_{em,k}$ is the electromagnetic energy spectral density defined by $\epsilon_{em,k} = \frac{1}{2}\epsilon_0(|E_{x,k}|^2 + |E_{y,k}|^2 + c^2|B_{z,k}|^2)$ and $\epsilon_{\phi,k}$ is a generalization of the phasestrophy flux density, for the mode $k$, defined by

$$ \epsilon_{\phi,k} = - \int d^3 p \frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2\partial p_x}. \qquad (18) $$

This quantity has been introduced by analogy with the “phasestrophy” integral, abbreviation of “phase-space entrophy,” namely, the scale-invariant self-similarity of $\delta f^2$, introduced in Ref. 73 for tokamak turbulence, which provides measure of a two-point phase space density correlation. In Ref. 74, it was shown how a slightly different form of Eq. (18), in which the derivative of $F_0$ was evaluated on the beam initial velocity [cf. Eq. (20) therein], provides a useful tool to quantify the momentum exchanged between electromagnetically interacting electron beams. In Ref. 1, Eq. (18) was used for the same purpose. Details of the derivation of Eq. (18) are given in Appendix B. From (B23), by neglecting the electric field contribution and in the absence of source term $S_{c,ck}$, Eq. (17) reads as after integrating over the $k$ space:

$$ \epsilon_m + \epsilon_\phi = \text{const.} \qquad (19) $$

It must be noted that Eq. (19) formally holds only in the linear limit, due to the $|\delta f| \ll F_0$ constraint, using which Eq. (17) is obtained. Nevertheless, it has been shown in Ref. 1 to provide a reliable measure of the momentum exchanged between electron beams also during the nonlinear dynamics (the same for the analogous integral named $\delta V_y$ of Ref. 74), since Eq. (18) is well defined regardless of the amplitude of $|\delta f|^2$. Equation (19) is the generalization to the electromagnetic case of the expression derived by Kruskal and Oberman in Ref. 75 (see also Refs. 76 and 77) for the “energy” conservation in the electrostatic case. In an equivalent way, a connection can be established between the different phases $\theta_{x,k}$, $\theta_{y,k}$, and $\theta_{z,k}$, of the electromagnetic field components and the phase $\phi_k$ introduced in the dynamics of the distribution function, which evolve following Eq. (B25) as

$$ \frac{1}{2} \epsilon_0 \left( |E_{x,k}|^2 \frac{\partial \theta_{x,k}}{\partial t} + |E_{y,k}|^2 \frac{\partial \theta_{y,k}}{\partial t} - c^2 |B_{z,k}|^2 \frac{\partial \theta_{z,k}}{\partial t} \right) + \int d^3 p \frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2\partial F_0} \left( \frac{\partial \phi_k}{\partial t} + \frac{\mathbf{k} \cdot \mathbf{p}}{m\gamma} \right) = S_{c,ck}. \qquad (20) $$

By neglecting the electric contribution of the electric field, one obtains from Eq. (20) the condition:

$$ -\frac{1}{2} \epsilon_0 c^2 |B_{z,k}|^2 \frac{\partial \theta_{z,k}}{\partial t} + \int d^3 p \frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2\partial F_0} \left( \frac{\partial \phi_k}{\partial t} + \frac{\mathbf{k} \cdot \mathbf{p}}{m\gamma} \right) = S_{c,ck}. \qquad (21) $$

The source terms $S_{c,ck}$ and $S_{c,sk}$ used in Eqs. (17) and (21) are given in the Appendix B by Eqs. (B22) and (B26). The system of Eqs. (17) and (21) allows to highlight the physical mechanisms leading to the reversal of the energy transfer observed in the system during the numerical simulations. We make the following assumptions:

* The magnetic contribution remains dominant over the electrical energy (given the results observed in the numerical simulations).

* The two source terms $S_{c,ck}$ and $S_{c,sk}$ 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 $\sum_j \alpha_j (C_j/m\gamma) F_{0j}$, which leads to the electric mean quasi-neutrality 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 $\epsilon_{m,k} + \epsilon_{\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 $\epsilon_{m,k} = \frac{1}{2} \epsilon_0 c^2 |B_{z,k}|^2$ and the phasestrophy$^{73}$ flux density $\epsilon_{\phi,k}$, satisfying Eq. (18). In the absence of the source term $S_{c,sk}$, Eq. (21) reads as

$$ \frac{1}{2} \epsilon_0 c^2 |B_{z,k}|^2 \frac{\partial \theta_{z,k}}{\partial t} = \int d^3 p \frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2\partial p_x} \left( \frac{\partial \phi_k}{\partial t} + \frac{\mathbf{k} \cdot \mathbf{p}}{m\gamma} \right). \qquad (22) $$

If $\partial \theta_{z,k}/\partial t \sim \omega_B$, where $\omega_B$ represents the characteristic frequency of the magnetic mode, i.e., the mode associated with the Weibel-type instability, we recognize two possible cases:

* The case of a stationary magnetic mode, with $\omega_B \sim 0$, induced by a CFI-type instability. The synchronization process takes here in a weak coupling regime, without resonance, to a partially synchronized state.

* The case of a Weibel-type 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 non-zero frequency $\omega_B \neq 0$. A synchronization induced by the resonant term $\partial \phi_k / \partial t + \mathbf{k} \cdot \mathbf{p} / m\gamma \approx 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 $\partial \phi_{k_0} (\mathbf{p}, t) / \partial t + \mathbf{p} \cdot \mathbf{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 $\epsilon_\phi$. Considering the invariance of $\epsilon_\phi + \epsilon_m = \text{const}$, approximatively valid also in the linear regime (and in the absence of source terms $S_{c,ck}$ and $S_{c,sk}$, cf. Fig. 6, below), thus implies a decrease of the magnetic energy, which is transformed into kinetic energy.

## IV. FREE-STREAMING FILAMENTATION VERSUS REINFORCED FILAMENTATION

In principle, the “reinforced” filamentation differs from the usual “free-particle” 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 coarse-graining 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, self-organization processes in the phase-space mediated
by synchronization of oscillatory modes (e.g., in a Kuramoto-like 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 Kuramoto-type 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 - 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,⁷⁸ 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 $\varphi(x, p_x, t)$ of the corresponding distribution obtained from a Hilbert transform of the distribution function: we clearly observe locking at a phase $\varphi = \pi/2$ in regions of phase space where the filament is stretching, while the phase remains close to zero in regions of near-zero 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) = \int dp_y f(x, L_y/2, p_x, p_y)$ in the $x - p_x$ phase-space 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\varphi(x, p_x, t)} + c.c.$ distribution is given by $a = \sqrt{g^2 + g_{Hilbert}^2}$, while the phase $\varphi$ obeys $\tan \varphi = g_{Hilbert}/g$.

The phase-locking displayed in Fig. 1 appears to be an intrinsic property of the free-streaming system, which completely follows from its kinematics evolution (i.e., it is mediated by the $v\partial f/\partial x$ term). In the presence of external forces, which induce modifications of the phase-space distribution function that go beyond the action of the $v\partial f/\partial x$ free-streaming term, the phase-locking can give rise to more complex self-organization processes, which may possibly lead to the

<figure><img src="image_15.png" /><figcaption>FIG. 1. Formation of phase-space filaments as an initial Gaussian distribution evolves according the free-streaming motion, i.e., the 2D advection equation $\partial f/\partial t + v \cdot \nabla f = 0$ (where $v = p/m_y$ is the electron velocity) in the projected phase space $(x\omega_p/c, p_x/mc)$ (on top). The corresponding map of the phase of the distribution function $\tilde{f}(x, p_x, t) = \int dp_y f(x, y = L_y/2, p_x, p_y)$ is plotted on bottom frame at the same time $t\omega_p = 32$. Note that the phase $\varphi(x, p_x, t)$ is locked at $\pi/2$ inside the filaments, even at a very weak level of the distribution function. The simulation has been performed using the VLEM solver.</figcaption></figure>

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, so-called 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⁷⁹ 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 reso-
nance 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
non-zero Lorentz force, leads to a profound change in the nature of
the filamentation of the distribution function in phase space, involving
a re-organization 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 semi-Lagrangian
numerical schemes

In this section, we first recall some general features of the integra-
tion scheme of Lagrangian (i.e., PIC) and semi-Lagrangian solvers,
especially in relation to the notion of phase-space 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 cascade-like
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 turbulence-like 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 free-streaming term
∂f /∂t + v∂f /∂x = 0, whose exact solution in the Fourier space reads
as f (k, v, t) = f (k, v, 0)eikvt 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 con-
cept of time splitting. This numerical scheme is now known as the
semi-Lagrangian (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 semi-Lagrangian 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:⁸⁵⁻⁸⁷ the particle-
in-cell (PIC) technique can be regarded as a discretization of the phase
space in terms of a superposition of moving elements, each represent-
ing a “cloud” of physical particles, usually referred as superparticles. In
PIC codes, superparticle trajectories are computed from the electro-
magnetic 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 col-
lective behavior of the plasma.

Free-streaming 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 well-adapted 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
large-k 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<sub>x</sub>N<sub>y</sub> ~ (L/Δx)<sup>d<sub>x</sub></sup> ~ 256<sup>2</sup>,
L being the typical system length along a given direction and d<sub>x</sub> the
spatial dimension of the problem (here d<sub>x</sub> = 2).

Both SL Vlasov and PIC codes are initialized with the same sam-
pling 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<sup>78</sup> N<sub>Vlas</sub> to the number of superparticles by cell N<sub>part</sub> (used
in the SMILEI code) is chosen as N<sub>Vlas</sub>/N<sub>part</sub> = g<sub>PIC</sub>N<sub>p</sub><sup>2</sup> ~ 1. This cor-
responds to a graininess parameter of the order of g<sub>PIC</sub> ~ 10<sup>-4</sup> for the
PIC solver SMILEI.<sup>88</sup> Although both numerical approaches solve the
Vlasov equation, the Eulerian or semi-Lagrangian description corre-
sponds to a statistical approach where the graininess parameter g ini-
tially tends toward zero. In the particles-in-cell approach, the
reintroduction of finite-size superparticles leads to a PIC specific grain-
iness parameter g<sub>PIC</sub> in principle small but finite. In the simulations we
have carried out, we have N<sub>p,x</sub>N<sub>p,y</sub> = N<sub>p</sub><sup>2</sup> = 128<sup>2</sup> in the momentum
space in the SL scheme, which is equivalent to use 128<sup>2</sup> particles by cell
in the SMILEI code (corresponding to a value of the inverse of the
graininess parameter of g<sub>PIC</sub> ≃ 6.1 × 10<sup>-5</sup>). The numerical parame-
ters used in the PIC/Vlasov comparison are detailed in Appendix C.
Space variables being normalized in d<sub>s</sub> = c/ω<sub>p</sub> units, times to the
inverse of the plasma frequency ω<sub>p</sub><sup>-1</sup>, and energy densities to n<sub>0</sub>mc<sup>2</sup>.

B. Vlasov-Maxwell modeling of oblique Weibel-type 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 Δx, Δy of the order of one-hundredth of the Debye length and an equivalent number of particles by cell (in the PIC code) or of grid point in the **p**-space (for the SL VLEM code) as discussed in Sec. III A.

The two symmetrical electron beams have normalized velocities
$\beta_1 = v_{y1}/c = -0.67$ and $\beta_2 = v_{y2}/c = 0.67$, counter-propagating in
the y direction. Charge and current neutrality are initially ensured by imposing $n_1 + n_2 = n_0$ (where $n_1 = n_2$) and $n_1\beta_1 + n_2\beta_2 = 0$, for all simulations corresponding to the non-propagative OI and to the propagative branch of OI. The initial distribution function is composed of two drifted Maxwellian distributions $f(\mathbf{x}, \mathbf{p}, t=0) = \sum_j n_j F_{0j}(\mathbf{p})$, where

$$F_{0j} = \frac{1}{2\pi\beta_{th}^2} \exp\left(-\frac{p_x^2 + (p_y + C_j)^2}{2m^2c^2\beta_{th}^2}\right), \quad (23)$$

where $\beta_{th} = \sqrt{k_B T_{eq}/mc^2}$ is the thermal velocity normalized to the light velocity c. For each propagating electron beam, we choose the same temperature $T_{eq} = 6keV$. To start the instability, a perturbation on a single mode ($n_x, n_y$) is introduced in the z-component of the magnetic field in the form of $\delta B_z = B_0 \sin(n_x \Delta k_x x + n_y \Delta k_y y)$. This allows us to excite the most unstable mode of wave vector components ($k_x, k_y, 0$) for each class of solution of the OI. The corresponding values are given in Table II. In both simulations, we use a time step of $\Delta t\omega_p = 0.005$.

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

<table><caption>**TABLE II.** Physical parameters used in numerical simulations performed for the study of the propagative branch, in the case of excitation of a single unstable wave number.</caption><thead><tr><th>Used code</th><th>Vlem (SL)</th><th>Smilei (PIC)</th></tr></thead><tbody><tr><td>Box length in x</td><td>L<sub>x</sub> = 2.51327d<sub>e</sub></td><td>L<sub>x</sub> = 2.51327d<sub>e</sub></td></tr><tr><td>Box length in y</td><td>L<sub>y</sub> = πd<sub>e</sub></td><td>L<sub>y</sub> = πd<sub>e</sub></td></tr><tr><td>Dominant wave vector</td><td>(k<sub>x</sub>d<sub>e</sub> = 5, k<sub>y</sub>d<sub>e</sub> = 2)</td><td>(k<sub>x</sub>d<sub>e</sub> = 5, k<sub>y</sub>d<sub>e</sub> = 2)</td></tr><tr><td>Corresponding mode numbers</td><td>(n<sub>x</sub> = 2, n<sub>y</sub> = 1)</td><td>(n<sub>x</sub> = 2, n<sub>y</sub> = 1)</td></tr><tr><td>Momentum sampling/ particle</td><td>N<sub>p_x</sub>N<sub>p_y</sub> = 128<sup>2</sup></td><td>N<sub>part by cell</sub> = 128<sup>2</sup></td></tr><tr><td>Graininess parameter</td><td>g<sub>PIC</sub> = 0</td><td>g<sub>PIC</sub> = 6.1 &times; 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>**FIG. 2.** Evolution of the kinetic $\epsilon_K$ (top frame), electric $\epsilon_e$ (middle frame), and magnetic $\epsilon_m$ (bottom frame) energy densities of the simulation case of Table II (evolution of a single propagative OI mode). The simulation has been performed using the (SL) VLEM solver.</figcaption></figure>

field, which is typical of the two-stream instability (TSI). In particular, the kinetic energy density $\epsilon_K = \int \frac{d^3x}{V} \int d^3p m c^2 (\gamma - 1) f$ (top frame in Fig. 2) is transferred to the magnetic energy density $\epsilon_m = \int \frac{d^3x}{V} \frac{1}{2} \epsilon_0 c^2 B_z^2$, shown in the bottom frame, and to the total electric energy density defined by $\epsilon_e = \int \frac{d^3x}{V} \frac{1}{2} \epsilon_0 (E_x^2 + E_y^2)$ (shown in the middle frame). In the evolution of the non-propagative OI, (i.e., with Re($\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 \sim 40\omega_p^{-1}$ with a maximum amplitude of the order of $\epsilon_m/n_0 m c^2 \simeq 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
<figure><img src="image_17.png" /><figcaption>**FIG. 3.** Evolution of the magnetic energy density $\epsilon_m$ (upper frames) of the corresponding x-component of the electric energy density $\epsilon_{ex}$ (middle frames) and of the y component of the electrostatic energy (lower frames). The left column corresponds to the simulation performed using the (SL) VLEM code, while the right column shows the results obtained from (PIC) SMILEI code. Note that the level of noise is bigger along the y component in the PIC numerical scheme in comparison to the semi-Lagrangian VLEM approach (middle frames).</figcaption></figure>

$t\omega_p \ge 100$. As we will see later, in a dedicated section, this process is linked to phase-space 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 $\epsilon_{ex}$ and $\epsilon_{ey}$ contributions of the total electric energy density $\epsilon_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 bottom-left 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. Contour-plots of the magnetic field

During the nonlinear stage of the OI evolution, a transition is observed toward a self-organized 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) z-component of the magnetic field $eB_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 self-organization begins with a transition toward the CFI-like 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 x-space direction, between two possible asymptotic states. In the CFI-like 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

<figure><img src="image_18.png" /><figcaption>FIG. 4. Series of contour-plots of the magnetic field component $eB_z/m\omega_p$ at different times chosen in the beginning of the saturation regime of OI (propagative branch). The corresponding time history of the magnetic energy density is shown in Fig. 2 (bottom frame). The simulation has been carried out using the SL VLEM code. The same simulation case studied using the PIC code is shown in Fig. 5. A toward a CFI-like magnetic pattern is clearly visible in the series of plots, which show that a homogeneous state along the y direction at time $t\omega_p = 40$ (bottom right frame), typical of the CFI, is attained. The simulation has been performed using the (SL) VLEM solver.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 5. Contour-plot of the magnetic field component $eB_z/m\omega_p$, using the same physical parameters of Fig. 4. Results have been obtained from a simulation performed with the PIC code SMILEI. The beginning of the simulation displays results very similar to those obtained with the semi-Lagrangian code (cf. Fig. 4). Here, the emphasis is on the long time dynamics of the magnetic field: we observe an oscillation between two different states characteristic of the OI and CFI magnetic field evolution, respectively. This oscillating behavior was also observed in the semi-Lagrangian simulation but with a smaller frequency.</figcaption></figure>

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 phaseatrophy flux density

As mentioned in relation to Eq. (18), the phaseatrophy $S_\phi = \int \frac{d^3x}{V} \int d^3p \frac{(\delta f)^2}{2\epsilon_{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 phase-space density auto-correlation 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 phaseatrophy flux density $\epsilon_{\phi,k}$ and the source terms $S_{c,c,k}$ and $S_{c,s,k}$. When the source terms are both zero, the resulting conditions (19), i.e., $\epsilon_m + \epsilon_\phi = \text{const}$, highlight an energy transfer occurring in the linear stage of the OI, which takes place in combination with the energy conservation $\epsilon_m + \epsilon_K = \text{const}$ (assuming that the electric energy density $\epsilon_e$ is negligible). This is also verified in the nonlinear regime. The variation of $\epsilon_\phi$, during the growth of the OI, can be interpreted in terms of momentum transfer,⁷⁴ a process already observed in tokamak turbulence,⁷³ of which we take $\epsilon_{\phi,k}$ as a measure for the “enhanced” filamentation mechanism associated with the
<figure><img src="image_20.png" /><figcaption>FIG. 6. Plot of the phasestrophy flux density $\epsilon_{\phi}$, of the magnetic energy density $\epsilon_m$, and of their sum $\epsilon_{\phi} + \epsilon_m$, vs time, in the linear regime of the OI instability. The simulation has been performed using the (SL) VLEM solver.</figcaption></figure>

mode *k*. In Fig. 6, we have shown in the same figure the time evolution of the (total) phasestrophy flux density $\epsilon_{\phi}$, of the magnetic energy density $\epsilon_m$, and of their sum $\epsilon_{\phi} + \epsilon_m$, in the linear phase of the OI, i.e., for a time $t\omega_p \le 35$. We observe that the quantity $\epsilon_{\phi} + \epsilon_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 $\int d^3\mathbf{p} \frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2\frac{\partial f_k}{\partial p_x}} \left( \frac{\partial \phi_k}{\partial t} + \frac{k\mathbf{p}}{m\gamma} \right)$ in the right-hand 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 $\frac{\partial \phi_k}{\partial t} + \frac{k\mathbf{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 entrophy 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

<figure><img src="image_21.png" /><figcaption>FIG. 7. Time history of the phasestrophy flux density $\epsilon_{\phi}$. This must be compared with the previous results in Fig. 1. Results were obtained from a simulation performed by the VLEM solver. Reproduced with permission from Ghizzo et al., Phys. Rev. Lett. 131, 035101 (2023). Copyright 2024 American Physical Society.</figcaption></figure>

the ordering of $|\delta f_k|^2/(\partial F_0/\partial 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 non-propagative 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_{e,c,k}$, which contains some terms like $\cos(\theta_{y,k} - \varphi_k)$ in Eq. (B10), becomes non-null]. 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 $|\epsilon_{\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 \sim 50$, is related to the transition to the CFI-like state, with $\epsilon_K \to \epsilon_m$. A second event takes place at a later time, for $t\omega_p \sim 150-250$, and is linked to the reversal of the energy transfer, i.e., $\epsilon_m \to \epsilon_K + \epsilon_{y\gamma}$: in the presence of an OI instability, the phases of the Van Kampen modes are locked (since the phase $\Theta_k - \varphi_k$ tend to zero, and the source terms $S_{e,c,k}$ and $S_{e,s,k}$ are modified), which leads to a growth of the quantity $|\epsilon_{\phi,k}|$.

The first growth in the phasestrophy density flux around time $t\omega_p \sim 40-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 \sim 150-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 right-hand side of Eq. (22), i.e., $(\frac{\partial \phi_k}{\partial t} + \frac{k\mathbf{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 left-hand term $\epsilon_m \frac{\partial \theta_{z,k}}{\partial t}$ of Eq. (22) remains finite, this transition is accompanied by a divergence of the term $\frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2\frac{\partial f_k}{\partial 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 phase-space 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-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 phase-space 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-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,
<figure><img src="image_22.png" /><figcaption>FIG. 8. Series of plots of the distribution function in the $x - p_x$ phase space in the case evolution of a single propagative OI mode (see Table II). The simulation has been carried out by the (SL) VLEM code.</figcaption></figure>

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 \sim 30-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

<figure><img src="image_23.png" /><figcaption>FIG. 9. Contour-plots of the distribution function in the $x - p_x$ phase space for the nonlinear evolution of a propagative OI mode (cf. Table II), obtained with the PIC SMILEI solver. These results highlight the (nonlinear) process of “enhanced filamentation” forced by phase synchronization of Van Kampen modes. A progressively thinner filament is produced during the stretching mechanism until it disappears at time $t\omega_p = 44$, on the top right frame. A critical layer structure is formed at time $t\omega_p = 46$ (middle right frame), which is characteristic of a phase synchronization process. Finally, at the X-point, on the bottom right frame, we observe the resurgence of a strong filamentation process. We have added the symbols of the X- and O-point in the bottom-left frame for the sake of clarity.</figcaption></figure>

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 “X-type” point (i.e., a bifurcation point in the phase space) disappears at $t\omega_p = 44$. The “O-type” point where there is a phase-space 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 X-type points appear to be unstable with respect to the phase-space 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 TURBULENCE-LIKE 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 mean-field 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 $-\mathbf{k}\cdot\mathbf{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 phase-space 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 turbulent-like scenario in which several modes are excited. We will also show that the nonlinear dynamics associated with this process may lead to a global self-organization of the system, due to a bifurcation-type mechanism, which can be modeled in the mean-field (Hamiltonian-type) Kuramoto approach.

We will show indeed that a global re-organization takes place in which collisionless wave-particle interactions transfer the “free” ordered kinetic energy stored in the high-energy 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 Weibel-type 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 re-converted into disordered kinetic energy, i.e., internal energy of the plasma. In the turbulence-like 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

A last simulation (we refer here as OI with broad spectrum) is now carried out, by using the VLEM code, using an initial magnetic field perturbation in the form

$$ \delta B_z = \sum_{n_x, n_y=1}^{15} B_0 \sin(n_x \Delta k_x x + n_y \Delta k_y y + \varphi_{n_x, n_y}), \quad (24) $$

where $\varphi_{n_x, n_y}$ is a random phase, and $\Delta k_x$ and $\Delta k_y$ are the fundamental wavevector components. Physical parameters, used in this simulation, are summarized in Table III.

<table><caption>**TABLE III.** Numerical and physical parameters used in the semi-Lagrangian VLEM code for a simulation case in which the magnetic perturbation is introduced to excite a broad spectrum in wavenumbers, by thus allowing for the coupling between stationary OI and propagative OI modes.</caption><thead><tr><th></th><th>Case OI with broad spectrum</th></tr></thead><tbody><tr><td>Used code</td><td>Vlem (SL)</td></tr><tr><td>Box length in x</td><td>L<sub>x</sub> = 2&pi;d<sub>e</sub></td></tr><tr><td>Box length in y</td><td>L<sub>y</sub> = 2&pi;d<sub>e</sub></td></tr><tr><td>Wave numbers excited</td><td>n<sub>x</sub>, n<sub>y</sub> &isin; [1, 15]</td></tr><tr><td>Momentum sampling/ particle</td><td>N<sub>p_x</sub>N<sub>p_y</sub> = 128<sup>2</sup></td></tr><tr><td>Magnetic perturbation amplitude</td><td>eB<sub>0</sub>/m&omega;<sub>p</sub> = 5 &times; 10<sup>-5</sup></td></tr><tr><td>Corresponding figures</td><td>Figs. 10–12</td></tr></tbody></table>

Figure 10 shows the temporal evolution of the different energy densities: kinetic $\epsilon_K$ (top), electric $\epsilon_e$ (middle), and finally magnetic $\epsilon_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 non-propagative 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 non-propagative OI branch occurring at $kd_e \ge 1$), which is dominated by electrostatic trapping and by the nonlinear excitation of increasingly larger wavenumbers, up to values $kd_e \gg 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^{-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 $\epsilon_K \to \epsilon_m + \epsilon_e$, which saturates at time $t\omega_p \sim 35$ in Fig. 10.

* step (III): the reversal of the energy transfer (i.e., $\epsilon_m \to \epsilon_K; \epsilon_e \ll \epsilon_K, \epsilon_m$) takes place in the third phase of the plasma dynamics, in the time interval $[40, 100]\omega_p^{-1}$ and is characterized by a slow decrease of both the electric $\epsilon_e$ and magnetic $\epsilon_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 low-frequency 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
<figure><img src="image_24.png" /><figcaption>FIG. 10. Illustration of the time evolution of counterstreaming beam configuration in which a wide spectrum of unstable modes is excited on the $B_z$ component (cf. Table III): the reversal of the energy transfer ($\epsilon_m + \epsilon_e \rightarrow \epsilon_K$) takes place just after the saturation regime of the OI. Almost all of the magnetic energy is eventually transferred back into kinetic energy in the form of stochastic heating. A residual part of the magnetic energy (bottom figure) but also of the electrostatic energy (middle figure) remains over a long timescale: the electromagnetic field is not completely attenuated. A large amount of the ordered kinetic energy initially stored in the two electron beams is transformed into heat during a stochastic process. The simulation has been performed with the VLEM code.</figcaption></figure>

growth of the most unstable OI modes. This is associated with an energy transfer from the kinetic $\epsilon_K$ to electromagnetic $\epsilon_e + \epsilon_m$ energy density, till saturation, observed at time $t\omega_p \sim 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 \times 10^{-4}$ and $1.20 \times 10^{-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, $eB_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 quasi-global (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 non-propagative 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 (bottom-left 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 low-k OI modes to a dominance of large-k 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 $p$-space. The distribution function is concentrated in space, i.e.,

<figure><img src="image_25.png" /><figcaption>FIG. 11. Illustration of the plasma dynamics for the simulation case of Table II. In the left column, plots of the magnetic field component $eB_z/m\omega_p$ at two different times showing the amplification and the growth of the magnetic field in the form of tilting stripes. These are characteristic features of the growth of an oblique instability. On the right column, the corresponding plot of the electron density in the configuration space is shown at the same times, exhibiting the emergence of small-scale (spatial) filamentation structure driven by the “reinforced” filamentation mechanism. The results were obtained from a simulation performed by the VLEM solver.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>FIG. 12. Series of representation of the distribution in the *p*-space at different times in simulation case of Table III, performed by using the VLEM code.</figcaption></figure>

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 bi-Maxwellian 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 - p_x$ or $y - p_y$ sub-phase 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 phase-space 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 semi-Lagrangian models, to show with numerical experiments how these features intervene in the nonlinear dynamics of oblique Weibel-type instabilities (Secs. V and VI). In this regard, we have also pointed out some specific features of the concept of “reinforced” or “enhanced” (phase-space) filamentation, which we have introduced in contraposition to the notion of

“kinematic” filamentation of a free-streaming 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 Kuramoto-type 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 $\epsilon_{\phi} + \epsilon_{m} = \text{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 turbulent-like energy cascade. This process takes place during the nonlinear dynamics of oblique Weibel-type 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 Weibel-type 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 $\text{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 Weibel-type modes. Interestingly, the non-trivial conversion from magnetic to kinetic energy is sometimes also observed in non-ideal 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 $z$-component 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 phase-mixing induced by phase-space 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 (FR-FCM) through the FR-FCM AAP-2023 Evolution of current sheets in low-collision plasmas. During the completion of this work, one of the authors (H. Betar) has received financial support from the AIM4EP
project (ANR-21-CE30-0018), 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.

## APPENDIX A: THE INFINITE-N KURAMOTO MODEL IN A HAMILTONIAN APPROACH

In providing a Hamiltonian description of a conservative system of N coupled oscillators, we may identify a new set of action-angle variables, say $J_n$, $\varphi_n$ for each oscillator $n$ with $1 \le n \le N$. Following Ref. 72, we can introduce a Hamiltonian in the form

$$ \mathcal{H} = \mathcal{H}_0 + \delta\mathcal{H} = \sum_n J_n \omega_n + Kr \sum_n J_n \sin(\bar{\varphi} - \varphi_n), \quad (A1) $$

where the first term in the r.h.s. describes the uncoupled oscillators having the natural frequencies $\omega_n$ while the second takes into account the coupling via the coupling parameter $K$. The mean-field (Hamiltonian) Kuramoto model reads

$$ \dot{J}_n = - \frac{\partial \mathcal{H}}{\partial \varphi_n} = Kr J_n \cos(\bar{\varphi} - \varphi_n), \quad (A2) $$

$$ \dot{\varphi}_n = \Omega_n = \frac{\partial \mathcal{H}}{\partial J_n} = \omega_n + Kr \sin(\bar{\varphi} - \varphi_n). \quad (A3) $$

Equation (A3) constitutes the Kuramoto equation for the phase of the n-th oscillator, while (A2) describes the dynamics of action $J_n$. Following Refs. 72 and 97 we can introduce and entropy $S_n$ for each oscillator $n$, in the form

$$ S_n = k_B \ln J_n, \quad (A4) $$

the total entropy of the system S being defined by the sum of all the $S_n$ terms. It can be shown (Ref. 72) that this model yields a positive entropy production, compatible with the violation of Casimirs induced by the coarse-graining of sub-grid phase-space filaments. In Refs. 1 and 52, a numerical decrease in entropy has been numerically measured, in relation to the violation of Casimirs induced by phase-space filamentation, which is apparently compatible with a self-organization process related to a Kuramoto-type global synchronization, discussed in Secs. III and VI. These numerical effects will be discussed in a forthcoming, dedicated work.

## APPENDIX B: DERIVATION OF THE INVARIANT QUANTITIES USED TO CHARACTERIZE THE "ENHANCED" FILAMENTATION

We here provide details about the derivation of invariant-type quantities in the linearization of Vlasov-Maxwell equations in the case of a counter-propagating beams system, which are introduced in Sec. III. We expand the distribution function in Fourier transform in $x$, through a linearization process around the initial equilibrium $F_0(\mathbf{p}) = n_0 \sum_{j=1,2} \alpha_j F_{0j}(\mathbf{p})$ in the form

$$ f = n_0 \sum_{j=1}^{2} \alpha_j F_{0j} + \sum_{k} \delta f_k(\mathbf{p}, t) \exp i(k_x x + k_y y), \quad (B1) $$

and the electromagnetic field component

$$ E = \sum_k E_k(t) \exp i(k_x x + k_y y) $$

and

$$ B_z = \sum_k B_{z,k}(t) \exp i(k_x x + k_y y). \quad (B2) $$

Linearizing the Vlasov equation (1) leads to

$$ \begin{aligned} \frac{\partial \delta f_k}{\partial t} = & -i \frac{\mathbf{k} \cdot \mathbf{p}}{m_\gamma} \delta f_k + e \left( \sum_{j \in \{x, y\}} \frac{p_j E_{j,k}}{m^2 c^2 \beta_j^2} \right) F_0(\mathbf{p}) \\ & + F_{y,k} \frac{n_0 \sum_{j=1,2} \alpha_j C_j F_{0j}}{m^2 c^2 \beta_y^2} + en_0 \left( \sum_{j=1,2} \frac{p_x p_y}{m^3 \gamma c^2} B_{z,k} \right) \\ & \times \sum_{j=1,2} \alpha_j F_{0j}(\mathbf{p}) \left( \frac{1}{\beta_x^2} - \frac{1}{\beta_y^2} \right), \end{aligned} \quad (B3) $$

where the Fourier component of the Lorentz force component $F_y$ is defined by $F_{y,k} = e(E_{y,k} - p_x B_{z,k}/m\gamma)$. Assuming that the current neutrality condition $n_0 \sum_j \alpha_j (C_j/m\gamma) F_{0j} = 0$ is initially satisfied, the Fourier components of linearly polarized electromagnetic field, i.e., the components $E_{x,k}$, $E_{y,k}$ and $B_{z,k}$, obey to

$$ \frac{\partial E_{x,k}}{\partial t} = ik_y c^2 B_{z,k} - \frac{e}{\varepsilon_0} \int d^3 p \frac{p_x}{m\gamma} \delta f_k(\mathbf{p}, t), \quad (B4) $$

$$ \frac{\partial E_{y,k}}{\partial t} = -ik_x c^2 B_{z,k} - \frac{e}{\varepsilon_0} \int d^3 p \frac{p_y}{m\gamma} \delta f_k(\mathbf{p}, t), \quad (B5) $$

$$ \frac{\partial B_{z,k}}{\partial t} = ik_y E_{x,k} - ik_x E_{y,k}. \quad (B6) $$

The last term in the right-hand side of (B3) disappears when $\beta_x = \beta_y = \beta_{th}$: this condition would correspond to suppress the Weibel instability, which is driven by a temperature anisotropy (we recall that, for counterstreaming beams along the direction $y$, the Weibel instability is excited when $\beta_x > \beta_y$). Thus, without loss of generality, we can consider only a situation of two counterstreaming electron beams having the equal (normalized) thermal velocities $\beta_{th}$ (i.e., with $\beta_{th} = v_{th}/c$ and $v_{th} = \sqrt{k_B T_e / m}$), supposing for them Maxwellian velocity distributions along $y$. The perturbed distribution function, $\delta f_k(\mathbf{p}, t)$, in the Fourier space reads
$$ \frac{\partial \delta f_k}{\partial t} = -i \frac{\mathbf{k} \cdot \mathbf{p}}{m_y} \delta f_k + e \sum_{j \in \{x,y\}} p_j E_{j,k} \sum_j \frac{\alpha_j F_{0j}}{m^2 c^2 \beta_{th}^2} + F_{y,k} \frac{n_0 \sum_j \alpha_j C_j F_{0j}}{m^2 c^2 \beta_{th}^2}. \quad (B7) $$

We now introduce the exponential notation $\delta f_k = |\delta f_k(\mathbf{p}, t)| \times \exp i\phi_k(\mathbf{p}, t)$, where $\phi_k$ is the phase and $|\delta f_k|$ is the modulus of the complex quantity $\delta f_k$, and a similar notation for the different electric field components $E_{j,k} = |E_{j,k}| \exp i\theta_{j,k}(t)$ for $j \in \{x, y\}$ and $B_{z,k} = |B_{z,k}| \exp i\theta_{z,k}(t)$ for the $B_z$ component of the magnetic field. Separating the real and imaginary parts in Eq. (B7) allows us to separate it in the two sets of equations:

$$ \frac{\partial |\delta f_k|}{\partial t} = \sum_j \frac{e \alpha_j F_{0j}}{m^2 c^2 \beta_{th}^2} \sum_{j \in \{x,y\}} p_j |E_{j,k}| \cos(\theta_{j,k} - \phi_k) + |F_{y,k}| \frac{n_0 \sum_j \alpha_j C_j F_{0j}}{m^2 c^2 \beta_{th}^2}, \quad (B8) $$

$$ \frac{\partial \varphi_k}{\partial t} |\delta f_k| = - \frac{\mathbf{k} \cdot \mathbf{p}}{m_y} |\delta f_k| + \sum_j \frac{e \alpha_j F_{0j}}{m^2 c^2 \beta_{th}^2} \sum_{j \in \{x,y\}} p_j |E_{j,k}| \sin(\theta_{j,k} - \varphi_k) + |F_{y,k}| \frac{n_0 \sum_j \alpha_j C_j F_{0j}}{m^2 c^2 \beta_{th}^2}. \quad (B9) $$

The real quantities $|F_{y,k}|_c$ and $|F_{y,k}|_s$ represent the “Lorentz force” contributions of “cosine” and “sine” functions:

$$ |F_{y,k}|_c = e \left( |E_{y,k}| \cos(\theta_{y,k} - \varphi_k) - \frac{p_x |B_{z,k}|}{m_y} \cos(\theta_{z,k} - \varphi_k) \right), \quad (B10) $$

$$ |F_{y,k}|_s = e \left( |E_{y,k}| \sin(\theta_{y,k} - \varphi_k) - \frac{p_x |B_{z,k}|}{m_y} \sin(\theta_{z,k} - \varphi_k) \right). \quad (B11) $$

To complete the analysis of the linearized Vlasov-Maxwell system, we must also determine the time evolution of the amplitudes and phases of the first-order perturbations of the electromagnetic field components. By thus separating the real and imaginary contributions of field components in Eqs. (B4)-(B6), we can write the field amplitudes as

$$ \frac{\partial |E_{x,k}|}{\partial t} = -k_y c^2 |B_{z,k}| \sin(\theta_{z,k} - \theta_{x,k}) - \frac{e}{\epsilon_0} \int d^3 p \frac{p_x}{m_y} |\delta f_k| \cos(\varphi_k - \theta_{x,k}), \quad (B12) $$

$$ \frac{\partial |E_{y,k}|}{\partial t} = k_x c^2 |B_{z,k}| \sin(\theta_{z,k} - \theta_{y,k}) - \frac{e}{\epsilon_0} \int d^3 p \frac{p_y}{m_y} |\delta f_k| \cos(\varphi_k - \theta_{y,k}), \quad (B13) $$

$$ \frac{\partial |B_{z,k}|}{\partial t} = -k_y |E_{x,k}| \sin(\theta_{x,k} - \theta_{z,k}) + k_x |E_{y,k}| \sin(\theta_{y,k} - \theta_{z,k}), \quad (B14) $$

and the corresponding phases as

$$ |E_{x,k}| \frac{\partial \theta_{x,k}}{\partial t} = k_y c^2 |B_{z,k}| \cos(\theta_{z,k} - \theta_{x,k}) - \frac{e}{\epsilon_0} \int d^3 p \frac{p_x}{m_y} |\delta f_k| \sin(\varphi_k - \theta_{x,k}), \quad (B15) $$

$$ |E_{y,k}| \frac{\partial \theta_{y,k}}{\partial t} = -k_x c^2 |B_{z,k}| \cos(\theta_{z,k} - \theta_{y,k}) - \frac{e}{\epsilon_0} \int d^3 p \frac{p_y}{m_y} |\delta f_k| \sin(\varphi_k - \theta_{y,k}), \quad (B16) $$

$$ |B_{z,k}| \frac{\partial \theta_{z,k}}{\partial t} = k_y |E_{x,k}| \cos(\theta_{x,k} - \theta_{z,k}) - k_x |E_{y,k}| \cos(\theta_{y,k} - \theta_{z,k}). \quad (B17) $$

Using the set of Eqs. (B12)-(B14), it is possible to determine the quantity $|E_{x,k}|\partial_t|E_{x,k}| + |E_{y,k}|\partial_t|E_{y,k}| + c^2|B_{z,k}|\partial_t|B_{z,k}|$, which can write as follows:

$$ \begin{aligned} |E_{x,k}| \frac{\partial |E_{x,k}|}{\partial t} + |E_{y,k}| \frac{\partial |E_{y,k}|}{\partial t} + c^2 |B_{z,k}| \frac{\partial |B_{z,k}|}{\partial t} \\ = - \frac{e}{\epsilon_0} \int d^3 p \left[ p_x |E_{x,k}| \cos(\varphi_k - \theta_{x,k}) + p_y |E_{y,k}| \cos(\varphi_k - \theta_{y,k}) \right] \frac{|\delta f_k|}{m_y}, \end{aligned} \quad (B18) $$

while the quantity into bracket is defined from (B8) to be

$$ \begin{aligned} & \sum_{j \in \{x,y\}} e p_j |E_{j,k}| \cos(\theta_{j,k} - \varphi_k) \\ = & \frac{m^2 c^2 \beta_{th}^2}{n_0 \sum_j \alpha_j F_{0j}} \left[ \frac{\partial |\delta f_k|}{\partial t} - \frac{n_0 \sum_j \alpha_j C_j F_{0j}}{m^2 c^2 \beta_{th}^2} |F_{y,k}|_c \right]. \end{aligned} \quad (B19) $$

This leads us to the final expression

$$ \begin{aligned} & \frac{\partial}{\partial t} \left( \frac{1}{2} \epsilon_0 |E_{x,k}|^2 + \frac{1}{2} \epsilon_0 |E_{y,k}|^2 + \frac{1}{2} \epsilon_0 c^2 |B_{z,k}|^2 \right) \\ = & - \int \frac{d^3 p}{m_y} \frac{m^2 c^2 \beta_{th}^2}{n_0 \sum_j \alpha_j F_{0j}} \left[ |\delta f_k| \frac{\partial |\delta f_k|}{\partial t} - |\delta f_k| |F_{y,k}|_c \frac{n_0 \sum_j \alpha_j C_j F_{0j}}{m^2 c^2 \beta_{th}^2} \right]. \end{aligned} \quad (B20) $$

This in turn leads us to the energy-phaseatrophy conservation law

$$ \frac{\partial}{\partial t} (\epsilon_{e,k} + \epsilon_{m,k} + \epsilon_{\phi,k}) = S_{e,c;k}, \quad (B21) $$

where the source term $S_{e,c;k}$ is defined by

$$ S_{e,c;k} = \int d^3 p \frac{|\delta f_k|}{n_0 \sum_j \alpha_j F_{0j}} \sum_j \frac{C_j}{m_y} F_{0j} |F_{y,k}|_c. \quad (B22) $$

In the absence of source term $S_{e,c;k} = 0$, and by integrating over the $k$ space, Eq. (B21) leads to the invariant quantity

$$ \epsilon_e + \epsilon_m + \epsilon_\phi = \text{const.} \quad (B23) $$

To complete the analysis, we can also calculate the quantity $|E_{x,k}|^2 \partial_t \theta_{x,k} + |E_{y,k}|^2 \partial_t \theta_{y,k} - c^2 |B_{z,k}|^2 \partial_t \theta_{z,k}$ by multiplying the set of Eqs. (B15)-(B17) defining the temporal dynamics of phases $\theta_{t,k}$, for $j \in \{x, y, z\}$, by the corresponding field component amplitudes. Combining the obtained equations with (B9), we obtain
$$ 
\begin{aligned}
& \frac{1}{2} \varepsilon_0 \left( |E_{x,k}|^2 \frac{\partial \theta_{x,k}}{\partial t} + |E_{y,k}|^2 \frac{\partial \theta_{y,k}}{\partial t} - c^2 |B_{z,k}|^2 \frac{\partial \theta_{z,k}}{\partial t} \right) \\
& = - \int \frac{d^3 p}{m\gamma} \left[ \frac{|\delta f_k|^2}{2} \times \left( \frac{\partial \varphi_k}{\partial t} + \frac{\mathbf{k} \cdot \mathbf{p}}{m\gamma} \right) \right] \frac{m^2 c^2 \beta_{th}^2}{n_0 \sum_j \alpha_j F_{0j}} \\
& \quad - \int d^3 p |\delta f_k|^2 |F_{y,k}|_s \frac{n_0 \sum_j \alpha_j C_j F_{0j}}{m\gamma \sum_j \alpha_j F_{0j}}.
\end{aligned}
\quad (B24) $$

Finally, we obtain

$$ 
\begin{aligned}
& \frac{1}{2} \varepsilon_0 \left( |E_{x,k}|^2 \frac{\partial \theta_{x,k}}{\partial t} + |E_{y,k}|^2 \frac{\partial \theta_{y,k}}{\partial t} - c^2 |B_{z,k}|^2 \frac{\partial \theta_{z,k}}{\partial t} \right) \\
& + \int d^3 p \frac{p_x}{m\gamma} \frac{|\delta f_k|^2}{2} \frac{\partial \varphi_k}{\partial p_x} \left( \frac{\partial \varphi_k}{\partial t} + \frac{\mathbf{k} \cdot \mathbf{p}}{m\gamma} \right) = S_{\varepsilon,s;k},
\end{aligned}
\quad (B25) $$

where the source term $S_{\varepsilon,s;k}$ is defined as

$$ S_{\varepsilon,s;k} = - \int d^3 p \frac{|\delta f_k|}{n_0 \sum_j \alpha_j F_{0j}} \sum_j \frac{C_j}{m\gamma} F_{0j} |F_{y,k}|_s. \quad (B26) $$

Equations (B23) and (B25) constitute the set of invariant quantities that are used in Sec. III to characterize the “reinforced” filamentation mechanism driven by the global synchronization of Van Kampen modes.

## APPENDIX C: PIC VERSUS VLASOV CODES

The numerical semi-Lagrangian scheme, used for the Vlasov code, solves the characteristic equations in the backward directions so as the Lagrangian scheme of the PIC code solves it for the single particle. Thus, it is clear that “pushing one particle” in the PIC approach or interpolating the distribution at a given (particle) location from the foot of the characteristics need the same numerical cost in Ref. 98, it was estimated it to be close to the ratio of two quantities $N_{Vlas}/N_{part}$, where $N_{Vlas}$ is the total number of mesh points in the phase space in the SL case, and $N_{part}$ is the total number of “superparticles.” We have

$$ N_{Vlas} = \left(\frac{L}{\Delta x}\right)^{d_x} N_p^{d_p}; \quad N_{part} = g_{PIC}^{-1} \left(\frac{L}{\Delta x}\right)^{d_x}. \quad (C1) $$

Note that this $g_{PIC} \equiv (L/\Delta x)^{d_x}/N_{part}$ provides an “ab initio” finite discrepancy between the Vlasov equation integrated by PIC schemes and its continuum, mathematical limit. Its effect on the filamentation dynamics and on synchronization is thus a priori not easily quantifiable, if not through comparison with Eulerian schemes (which are however affected by the usual numerical limitations imposed by the grid sampling).

Let us then suppose (e.g., based on our “experience” with the VLEM code) that for an accurate enough description of the filamentation in the one-dimensional problem of interest here, we need $\Delta x \sim \lambda_{De}/100$ and a momentum sampling in momentum space, at least of $N_p \sim 100$. This constraint is severe and imposes to choose a large number of particles in the PIC code, i.e., $N_{part} \sim N_{Vlas}$ to so to guarantee an equivalent numerical accuracy (at least initially) in both models (see Sec. V B). We see that $N_{Vlas}/N_{part}$ can be written

<table><caption>TABLE IV. Computational cost: PIC code vs Vlasov code: ratio $N_{Vlas}/N_{part}$ for a sampling of $N_p = 100$ points in the momentum space.</caption><tr><td>g<sub>PIC</sub></td><td>10<sup>-2</sup></td><td>10<sup>-4</sup></td><td>10<sup>-6</sup></td></tr><tr><td>N<sub>Vlas</sub>/N<sub>part</sub> = g<sub>PIC</sub>N<sub>p</sub><sup>2</sup></td><td>100</td><td>1</td><td>10<sup>-2</sup></td></tr></table>

as a function of the quantity $g_{PIC} = 1/(n_0 \lambda_{De}^{d_x})$, which represents the graininess parameter (related to the particle discreteness of the PIC approach (for $d_x = 2$ and $d_p = 2$), $N_{part}$ can now be put in the form $N_{part} = g_{PIC}^{-1}(L/\Delta x)^{d_x}$, leading to

$$ \frac{N_{Vlas}}{N_{part}} = g_{PIC} N_p^{d_p} = g_{PIC} N_p^2. \quad (C2) $$

In Table IV, we give the expected numerical cost between the two different approaches for a physical problem, characterized by a typical graininess parameter $g_{PIC}$.

## REFERENCES

<footnote n="1">¹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="2">²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="3">³A. Bret, I. Gremillet, D. Benisti, 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="4">⁴A. Bret, L. Gremillet, and D. Benisti, “Exact relativistic kinetic theory of the full unstable spectrum of an electron-beam-plasma system with Maxwell–Juttner distribution functions,” Phys. Rev E **81**, 036402 (2010).</footnote>

<footnote n="5">⁵A. Bret, “Weibel, two-stream, filamentation, oblique, bell, Buneman…Which one grows faster,” Astro. Phys. J. **699**, 990 (2009).</footnote>

<footnote n="6">⁶A. Bret and C. Deutsch, “Beam plasma electromagnetic instabilities in a smooth density gradient: Application to the fast ignition scenario,” Phys. Plasmas **12**, 102702 (2005).</footnote>

<footnote n="7">⁷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="8">⁸L. O. Silva, R. A. Fonseca, J. W. Tonge, W. B. Mori, and J. M. Dawson, “On the role of the purely transverse Weibel instability in fast ignitor scenarios,” Phys. Plasmas **9**, 2458 (2002).</footnote>

<footnote n="9">⁹A. Achterberg and J. Wiersma, “The Weibel instability in relativistic plasmas. I. Linear theory,” Astron. Astrophys. **475**, 1 (2007).</footnote>

<footnote n="10">¹⁰P. H. Yoon, “Electromagnetic Weibel instability in a fully relativistic bi-Maxwellian plasma,” Phys. Fluids B **1**, 1336 (1989).</footnote>

<footnote n="11">¹¹P. H. Yoon, “Relativistic Weibel instability,” Phys. Plasmas **14**, 024504 (2007).</footnote>

<footnote n="12">¹²U. Schaefer-Rolffs, I. Lerche, and R. Schlickeiser, “The relativistic kinetic Weibel instability: General arguments and specific illustrations,” Phys. Plasmas **13**, 012107 (2006).</footnote>

<footnote n="13">¹³R. Schlickeiser, “Covariant kinetic dispersion theory of linear waves in anisotropic plasmas. I. General dispersion relations, bi-Maxwellian distribution and nonrelativistic limits,” Phys. Plasmas **11**, 5532 (2004).</footnote>

<footnote n="14">¹⁴U. Schaefer-Rolffs and R. Schlickeiser, “Covariant kinetic dispersion theory of linear waves in anisotropic plasmas. II. Comparison of covariant and noncovariant growth rates of the nonrelativistic Weibel instability,” Phys. Plasmas **12**, 022104 (2005).</footnote>

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

<footnote n="16">¹⁶R. Lee and M. Lampe, “Electromagnetic instabilities, filamentation and focusing of relativistic electron beams,” Phys. Rev. Lett. **31**, 1390 (1973).</footnote>

<footnote n="17">¹⁷F. Califano, N. Attico, F. Pegoraro, G. Bertin, and S. V. Bulanov, “Fast formation of magnetic island in a plasma in the presence of counterstreaming electrons,” Phys. Rev. Lett. **86**, 5293 (2001).</footnote>
<footnote n="18">¹⁸M. Tzoufras, C. Ren, F. S. Tsung, J. W. Tonge, W. M. Pori, 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="19">¹⁹J. R. Peterson, S. Glenzer, and F. Fiua, "Magnetic field amplification by a non-linear electron streaming instability," Phys. Rev. Lett. 126, 215101 (2021).</footnote>

<footnote n="20">²⁰A. Bret, M. E. Dieckmann, and C. Deutsch, "Oblique electromagnetic instabilities for a hot relativistic beam interacting with a hot and magnetized plasma," Phys. Plasmas 13, 082109 (2006).</footnote>

<footnote n="21">²¹L. Gremillet, D. Benisti, E. Lefebvre, and A. Bret, "Linear and nonlinear development of oblique beam-plasma instabilities in the relativistic regime," Phys. Plasmas 14, 040704 (2007).</footnote>

<footnote n="22">²²M. V. Medvedev and A. Loeb, "Generation of magnetic fields in relativistic shock of gamma-ray burst sources," Astrophys. J. 526, 697 (1999).</footnote>

<footnote n="23">²³C. M. Huntington, F. Fiua, J. S. Ross, A. B. Zylstra, R. P. Drake, D. H. Froula, G. Gregori, N. L. Kugland, C. C. Kuranz, M. C. Levy, C. K. Li, J. Meinecke, T. Morita, R. Petrasso, C. Plechaty, B. A. Renington, D. D. Ryutoc, Y. Sakawa, A. Splitkovsky, H. Takabe, and H. S. Park, "Observation of magnetic field generation via the Weibel instability in interpenetrating plasma flows," Nat. Phys. 11, 173 (2015).</footnote>

<footnote n="24">²⁴G. F. Swadling, C. Brulsema, F. Fiua, D. P. Higginson, C. M. Huntington, H. S. Park, B. B. Pollock, W. Rozmus, H. G. Rinderknecht, J. Katz, A. Birkel, and J. S. Ross, "Measurement of kinetic-scale current filamentation dynamics and associated magnetic fields in interpenetrating plasmas," Phys. Rev. Lett. 124, 215001 (2020).</footnote>

<footnote n="25">²⁵M. Yamada, "Progress in understanding magnetic reconnection in laboratory and space astrophysical plasmas," Phys. Plasmas 14, 058102 (2007).</footnote>

<footnote n="26">²⁶M. Yamada, R. Kulsrud, and H. Ji, "Magnetic reconnection," Rev. Mod. Phys. 82, 603 (2010).</footnote>

<footnote n="27">²⁷H. J. Cai, D. Q. Ding, and L. C. Lee, "Momentum transport near a magnetic X line in collisionless reconnection," J. Geophys. Res.: Space Phys. 99(A1), 35, https://doi.org/10.1029/93JA02519 (1994).</footnote>

<footnote n="28">²⁸M. Hoshino, T. Mukai, T. Terasawa, and I. Shinohara, "Suprathermal electron acceleration in magnetic reconnection," J. Geophys. Res. Space Phys. 106, 25979, https://doi.org/10.1029/2001JA0900052 (2001).</footnote>

<footnote n="29">²⁹X. R. Fu, Q. M. Lu, and S. Wang, "The process of electron acceleration during collisionless magnetic reconnection," Phys. Plasmas 13, 012309 (2006).</footnote>

<footnote n="30">³⁰J. T. Dahlin, J. F. Drake, and M. Swisdak, "The mechanisms of electron heating and acceleration during magnetic reconnection," Phys. Plasmas 21, 092304 (2014).</footnote>

<footnote n="31">³¹S. R. Totorica, T. Abel, and F. Fiua, "Nonthermal electron energization from magnetic reconnection in laser-driven plasmas," Phys. Rev. Lett. 116, 095003 (2016).</footnote>

<footnote n="32">³²F. Pucci, S. Usami, H. Ji, X. Guo, R. Horiuchi, S. Okamura, W. Fox, J. Jara-Almonte, M. Yamada, and J. Yoo, "Energy transfer and electron energization in collisionless magnetic reconnection for different guide-field intensities," Phys. Plasmas 25, 122111 (2018).</footnote>

<footnote n="33">³³A. Chien, L. Gao, S. Zhang, H. Ji, E. G. Blackman, W. Daughton, A. Stanier, A. Le, F. Guo, R. Follett et al., "Non-thermal electron acceleration from magnetically driven reconnection in a laboratory plasma," Nat. Phys. 19, 254 (2023).</footnote>

<footnote n="34">³⁴R. Dandyopadhyay, W. H. Matthaeus, T. N. Parashar, Y. Yang, A. Chasapis, B. L. Giles, D. J. Gershman, C. J. Pollock, C. T. Russel, R. J. Strangeway, R. B. Torbert, T. E. Moore, and J. L. Burch, "Statistics of kinetic dissipation in the Earth's magnetosheath: MMS observations," Phys. Rev. Lett. 124, 255101 (2020).</footnote>

<footnote n="35">³⁵J. P. Eastwood, M. V. Goldman, T. D. Phan, J. E. Stawarz, P. A. Cassak, J. F. Drake, D. Newman, B. Lavraud, M. A. Shay, R. E. Regun, J. L. Burch, D. J. Gershman, B. L. Giles, P. A. Lindqvist, R. B. Torbert, R. J. Strangeway, and C. T. Russel, "Energy flux densities near the electron dissipation region in asymmetric magnetopause reconnection," Phys. Rev. Lett. 125, 265102 (2020).</footnote>

<footnote n="36">³⁶W. H. Matthaeus and S. L. Lamkin, "Turbulent magnetic reconnection," Phys. Fluids 29, 2513 (1986).</footnote>

<footnote n="37">³⁷H. Karimabadi, W. Daughton, and K. B. Quest, "Role of electron temperature anisotropy in the onset of magnetic reconnection," Geophys. Res. Lett. 31, L18801, https://doi.org/10.1029/2004GL020791 (2004).</footnote>

<footnote n="38">³⁸D. Del Sarto, F. Pegoraro, and F. Califano, "Pressure anisotropy and small spatial scales induced by velocity shear," Phys. Rev. E 93, 053203 (2016).</footnote>

<footnote n="39">³⁹D. Del Sarto and F. Pegoraro, "Shear-induced pressure anisotropization and correlation with fluid vorticity in a low collisionality plasma," Mon. Not. R. Astron. Soc. 475, 181 (2018).</footnote>

<footnote n="40">⁴⁰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).</footnote>

<footnote n="41">⁴¹F. Pucci, W. H. Matthaeus, A. Chasapis, S. Servidio, L. Sorriso-Valvo, V. Olshevsky, D. L. Newman, M. V. Goldman, and G. Lapenta, "Generation of turbulence in colliding reconnection jets," Astrophys. J. 867, 10 (2018).</footnote>

<footnote n="42">⁴²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).</footnote>

<footnote n="43">⁴³W. H. Matthaeus, "Turbulence in space plasmas: Who needs it?," Phys. Plasmas 28, 032306 (2021).</footnote>

<footnote n="44">⁴⁴R. Bandyopadhyay, Y. Yang, W. H. Matthaeus, T. N. Parashar, V. Roytershteyn, A. Chasapis, D. J. Gershman, B. L. Giles, and J. L. Burch, "Collisional-like dissipation in collisionless plasmas," Phys. Plasmas 30, 080702 (2023).</footnote>

<footnote n="45">⁴⁵Y. V. Khotyaintsev, D. B. Graham, K. Steinvall, L. Alm, A. Vaivads, A. Johlander, C. Norgren, W. Li, A. Divin, H. S. Fu, K. J. Hwang, J. L. Burch, N. Ahmadi, O. Le Contel, D. J. Gershman, C. T. Russel, and R. B. Torbert, "Electron heating by Debye-scale turbulence in guide-field reconnection," Phys. Rev. Lett. 124, 045101 (2020).</footnote>

<footnote n="46">⁴⁶L. Franci, S. S. Cerri, F. Califano, S. Landi, E. Papini, A. Verdini, L. Matteini, F. Jenko, and P. Hellinger, "Magnetic reconnection as a driver for a sub-ion-scale cascade in plasma turbulence," Astrophys. J. Lett. 850, L16 (2017).</footnote>

<footnote n="47">⁴⁷S. S. Cerri and F. Califano, "Reconnection and small-scale fields in 2D-3V hybrid-kinetic driven turbulence simulations," New J. Phys. 19, 025007 (2017).</footnote>

<footnote n="48">⁴⁸N. F. Louriero and S. Boldyrev, "Role of magnetic reconnection in magnetohydrodynamic turbulence," Phys. Rev. Lett. 118, 245101 (2017).</footnote>

<footnote n="49">⁴⁹A. Tenerani and M. Velli, "Spectral signatures of recursive magnetic field reconnection," Mon. Not. R. Astron. Soc. 491, 4267 (2020).</footnote>

<footnote n="50">⁵⁰S. S. Cerri, M. Kunz, and F. Califano, "Dual phase-space cascades in 3D hybrid-Vlasov-Maxwell turbulence," Astrophys. J. Lett. 856, L13 (2018).</footnote>

<footnote n="51">⁵¹A. Ghizzo and D. D. Sarto, "Low- and high-frequency nature of oblique filamentation modes. I. Linear theory," Phys. Plasmas 27, 072103 (2020).</footnote>

<footnote n="52">⁵²A. Ghizzo and D. D. 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="53">⁵³N. G. Van Kampen, "On the theory of stationary waves in plasmas," Physica 21, 949 (1955).</footnote>

<footnote n="54">⁵⁴R. J. DiPerna and P. L. Lions, "Global weak solutions of Vlasov-Maxwell systems," Commun. Pure Appl. Math. 42, 729 (1989).</footnote>

<footnote n="55">⁵⁵G. Rein, "Global weak solutions to the relativistic Vlasov-Maxwell system revisited," Commun. Math. Sci. 2, 145 (2004).</footnote>

<footnote n="56">⁵⁶C. Villani, "Particle systems and nonlinear Landau damping," Phys. Plasmas 21, 030901 (2014).</footnote>

<footnote n="57">⁵⁷J. H. Malmberg, C. B. Wharton, R. W. Gould, and T. M. O'Neil, "Observation of plasma wave echoes," Phys. Fluids 11, 1147 (1968).</footnote>

<footnote n="58">⁵⁸F. C. Grant and M. R. Feix, "Transition between Landau and Van Kampen treatments of the Vlasov equation," Phys. Fluids 10, 696 (1967).</footnote>

<footnote n="59">⁵⁹F. C. Grant and M. R. Feix, "Fourier-Hermite solutions of the Vlasov equations in the linearized limit," Phys. Fluids 10, 1356 (1967).</footnote>

<footnote n="60">⁶⁰G. Knorr, "Plasma simulation with few particles," J. Comput. Phys. 13, 165 (1973).</footnote>

<footnote n="61">⁶¹P. Bertrand, D. Del Sarto, and A. Ghizzo, The Vlasov Equation I: History and General Properties (ISTE Ltd, London and Wiley and Sons, 2019).</footnote>

<footnote n="62">⁶²A. J. Klimas, "A method for overcoming the velocity space filamentation problem in collisionless plasma model solution," J. Comput. Phys. 68, 202 (1987).</footnote>

<footnote n="63">⁶³A. F. Vinas and A. J. Klimas, "Flux-balance Vlasov simulation with filamentation filtration," J. Comput. Phys. 375, 983 (2018).</footnote>

<footnote n="64">⁶⁴Y. Kuramoto, Chemical Oscillations, Waves and Turbulence (Springer-Verlag, Berlin, 1984).</footnote>

<footnote n="65">⁶⁵J. A. Acebron, L. L. Bonilla, C. J. Perez Vicente, F. Ritort, and R. Spligler, "The Kuramoto model: A simple paradigm for synchronization phenomena," Rev. Mod. Phys. 77, 137 (2005).</footnote>

<footnote n="66">⁶⁶S. H. Strogatz, "From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators," Phys. D 143, 1 (2000).</footnote>
<footnote n="67">⁶⁷S. H. Strogatz, R. E. Mirollo, and P. C. Matthews, “Coupled nonlinear oscillators below the synchronization threshold: Relaxation by generalized Landau damping,” *Phys. Rev. Lett.* **68**, 2730 (1992).</footnote>

<footnote n="68">⁶⁸Z. B. Guo, Y. H. Wang, and S. K. Xu, “Phase condensation and evaporation: Another look at wave-particle interactions and Landau damping,” *Phys. Rev. E* **101**, 030201 (2010).</footnote>

<footnote n="69">⁶⁹S. Xu, Z. B. Guo, and O. D. Gurcan, “Transition from linear Landau damping to nonlinear Bernstein–Greene–Kruskal modes via phase synchronization,” *Phys. Rev. E* **103**, 023208 (2021).</footnote>

<footnote n="70">⁷⁰A. Ghizzo and D. D. Sarto, “The model of particles modes I. A paradigm for phase synchronization in tokamak turbulence,” *Phys. Plasmas* **29**, 042506 (2022).</footnote>

<footnote n="71">⁷¹A. Ghizzo and D. D. Sarto, “The model of particles modes II. Transition to a fishbone-like state triggered by global synchronization and energetic particles,” *Phys. Plasmas* **29**, 042507 (2022).</footnote>

<footnote n="72">⁷²V. Garcia-Morales, J. Pellicer, and J. A. Manzanares, “Thermodynamics based on the principle of least abbreviated action: Entropy production in a network of coupled oscillators,” *Ann. Phys.* **323**, 1844 (2008).</footnote>

<footnote n="73">⁷³P. H. Diamond, S. I. Itoh, and K. Itoh, *Modern Plasma Physics, Volume 1, Physical Kinetics of Turbulent Plasmas* (University Press, Cambridge, 2010), p. 309.</footnote>

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

<footnote n="75">⁷⁵M. D. Kruskal and R. C. Oberman, “On the stability of plasma in static equilibrium,” *Phys. Fluids* **1**, 275 (1958).</footnote>

<footnote n="76">⁷⁶P. J. Morrison, “The energy of perturbations for Vlasov plasmas,” *Phys. Plasmas* **1**, 1447 (1994).</footnote>

<footnote n="77">⁷⁷P. J. Morrison and B. A. Shadwick, “Canonization and diagonalization of an infinite dimensional noncanonical Hamiltonian system: Linear Vlasov theory,” *Acta Phys. Pol. A: Gen. Phys.* **85**, 759 (1994).</footnote>

<footnote n="78">⁷⁸M. Sarrat, A. Ghizzo, D. Del Sarto, and L. Serrat, “Parallel implementation of a relativistic semi-Lagrangian Vlasov–Maxwell solver,” *Eur. Phys. J. D* **71**, 271 (2017).</footnote>

<footnote n="79">⁷⁹M. C. Mackey, “The dynamic origin of increasing entropy,” *Rev. Mod. Phys.* **61**, 981 (1989).</footnote>

<footnote n="80">⁸⁰J. Barre, D. Metivier, and Y. Y. Yamaguchi, “Towards a classification of bifurcations in Vlasov equations,” *Phys. Rev. E* **102**, 052208 (2020).</footnote>

<footnote n="81">⁸¹C. Z. Cheng and G. Knorr, “The integration of the Vlasov equation in configuration space,” *J. Comput. Phys.* **22**, 330 (1976).</footnote>

<footnote n="82">⁸²M. L. Begue, A. Ghizzo, and P. Bertrand, “Two-dimensional Vlasov simulation of Raman scattering and plasma beatwave acceleration on parallel computers,” *J. Comput. Phys.* **151**, 458 (1999).</footnote>

<footnote n="83">⁸³A. Ghizzo, F. Huot, and P. Bertrand, “A non-periodic 2D semi-Lagrangian Vlasov code for laser-plasma interaction on parallel computer,” *J. Comput. Phys.* **186**, 47 (2003).</footnote>

<footnote n="84">⁸⁴E. Sonnendrucker, J. Roche, P. Bertrand, and A. Ghizzo, “The semi-Lagrangian method for the numerical resolution of the Vlasov equations,” *J. Comput. Phys.* **149**, 201 (1998).</footnote>

<footnote n="85">⁸⁵J. M. Dawson, “Computer modeling of plasma: Past, present and future,” *Phys. Plasmas* **2**, 2189 (1994).</footnote>

<footnote n="86">⁸⁶J. M. Dawson, “Particle simulation plasmas,” *Rev. Mod. Phys.* **55**, 403 (1983).</footnote>

<footnote n="87">⁸⁷C. K. Birdsall and A. B. Langdon, *Plasma Physics via Computer Simulation* (Ins. of Phys. Publishing, Bristol, Philadelphia, 1991).</footnote>

<footnote n="88">⁸⁸J. Derouillat, A. Beck, F. Perez, T. Vinci, M. Chiaramello, A. Grassi, M. Fle, G. Bouchard, I. Plotnikov, N. Aunai, J. Dargent, C. Riconda, and M. Grech, “SMILEI: A collaborative, open-source, multi-purpose particle-in-cell code for plasma simulation,” *Comput. Phys. Commun.* **222**, 351 (2018).</footnote>

<footnote n="89">⁸⁹N. J. Balmforth and R. Sassi, “A shocking display of synchrony,” *Phys. D* **143**, 21 (2000).</footnote>

<footnote n="90">⁹⁰P. A. Cassak and M. H. Barbhuiya, “Pressure-strain interaction. I. On compression, deformation, and implications for Pi-D,” *Phys. Plasmas* **29**, 122306 (2022).</footnote>

<footnote n="91">⁹¹P. A. Cassak, M. Barbhuiya, H. H. Liang, and M. R. Argall, “Quantifying energy conversion in higher-order phase space density moments in plasmas,” *Phys. Rev. Lett.* **130**, 085201 (2023).</footnote>

<footnote n="92">⁹²R. Pellat, “About ‘reconnection’ in a collisionless plasma,” *Space Sci. Rev.* **23**, 359 (1979).</footnote>

<footnote n="93">⁹³L. M. Zeleny and A. L. Taktakishvili, “Spontaneous magnetic reconnection mechanisms in plasma,” *Astrophys. Space Sci.* **134**, 185 (1987).</footnote>

<footnote n="94">⁹⁴H. J. Ziegler and H. Wiechen, “Mixing and relaxation in ideal incompressible fluids,” *Phys. Scr.* **T74**, 50 (1998).</footnote>

<footnote n="95">⁹⁵H. Wiechen and H. J. Ziegler, “Magnetic reconnection: On new aspects of the microscopic cause of localized dissipation,” *Phys. Scr.* **T74**, 14 (1998).</footnote>

<footnote n="96">⁹⁶D. Grasso, F. Califano, F. Pegoraro, and F. Porcelli, “Phase-mixing and island saturation in Hamiltonian reconnection,” *Phys. Rev. Lett.* **86**, 5051 (2001).</footnote>

<footnote n="97">⁹⁷A. C. Kalloniatis, “Entropy and stability of phase synchronisation of oscillators on networks,” *Ann. Phys.* **348**, 127 (2014).</footnote>

<footnote n="98">⁹⁸M. R. Feix and P. Bertrand, “A universal model: The Vlasov equation,” *Trans. Theory, Stat. Phys.* **34**, 1 (2005).</footnote>