2.1. Basic equations

The electron distribution function $f=f(x,y,\boldsymbol{p},t)$ obeys the relativistic Vlasov equation

(2.1) \begin{equation} \frac {\partial f}{\partial t}+\frac {\boldsymbol{p}}{m\gamma }\boldsymbol{\cdot}\frac {\partial f}{\partial \boldsymbol{x}}+e\left (\boldsymbol{E}+\frac {\boldsymbol{p\times B}}{m\gamma }\right )\boldsymbol{\cdot}\frac {\partial f}{\partial \boldsymbol{p}}=0, \end{equation}

where the Lorentz factor is given by

(2.2) \begin{equation} \gamma =\sqrt {1+\frac {\boldsymbol{p}^{2}}{m^{2}c^{2}}}, \end{equation}

and $\boldsymbol{p}=m\gamma \boldsymbol{v}$ is the electron kinematic momentum. We choose the Coulomb gauge ( $\mathrm{div}\boldsymbol{A}=0$ , where $\boldsymbol{A}$ is the vector potential) and ions are considered just as a fixed neutralising background. The Vlasov equation (2.1) is self-consistently coupled with the Maxwell equations

(2.3) \begin{equation} \frac {\partial \boldsymbol{E}}{\partial t}-c^{2}\mathrm{rot}\boldsymbol{B}+\frac {\boldsymbol{J}}{\varepsilon _{0}}=0\quad \mathrm{and}\quad \frac {\partial \boldsymbol{B}}{\partial t}+\mathrm{rot}\boldsymbol{E}=0, \end{equation}

together with the Poisson law $\mathrm{div}\boldsymbol{E}=\rho /\varepsilon _{0}$ and the condition $\mathrm{div}\boldsymbol{B}=0$ , where the electron current density $\boldsymbol{J}$ and the charge density $\rho$ are defined by

(2.4) \begin{equation} \boldsymbol{J}=e\int \frac {\boldsymbol{p}}{m\gamma }fd^{3}p\quad \mathrm{and}\quad \rho =e\int fd^{3}p-en_{0}. \end{equation}

Here $n_{0}$ , $e\lt 0$ and $m$ are the fixed background ion density, the elementary charge and electron rest mass, respectively.

2.2. The full kinetic VLEM approach

We first recall the global numerical scheme of the VLEM code implemented on a parallel computer in 1-D or 2-D Cartesian geometry and in 2-D or 3-D momentum space.

From a numerical point of view, the numerical integration of the Vlasov equation (2.1) is based on the use of a semi-Lagrangian technique, detailed in Sonnendrucker et al. (Reference Sonnendrucker, Roche, Bertrand and Ghizzo1999), Begue, Ghizzo & Bertrand (Reference Begue, Ghizzo and Bertrand1999), Ghizzo, Huot & Bertrand (Reference Ghizzo, Huot and Bertrand2003) and Sarrat et al. (Reference Sarrat, Ghizzo, Del Sarto and Serrat2017), allowing the integration of (2.1) along its characteristics. The semi-Lagrangian Vlasov solver implies two main steps.

1. (i) The use of a time-splitting scheme, which allows us to separate the numerical integration in configuration space (i.e. $\boldsymbol{x}$ ) and momentum space ( $\boldsymbol{p}$ ) by solving the following sequence of equations:

(2.5) \begin{equation} \mathrm{step\,(a){:}\;}\frac {\partial f}{\partial t}+\frac {\boldsymbol{p}}{m\gamma }.\frac {\partial f}{\partial \boldsymbol{x}}=0\;\mathrm{for}\;t\in {{\left [t_{n},t_{n+1/2}^{}\right ]}}, \end{equation}

(2.6) \begin{equation} \mathrm{step}\,(b){:}\;\frac {\partial f}{\partial t}+e\left (\boldsymbol{E}+\frac {\boldsymbol{p\times B}}{m\gamma }\right )\boldsymbol{\cdot}\frac {\partial f}{\partial \boldsymbol{p}}=0\;\mathrm{for}\;{{t\in \left [t_{n}^{},t_{n+1}^{}\right ]}}.\\[12pt] \end{equation}

After, step (a) is solved again for a time $t$ in the interval $[t_{n+1/2}^{},t_{n+1}]$ .

1. (ii) For each advection equation (2.5) or (2.6), the integration is realised along the backward characteristics by researching the feet of the characteristics. The Vlasov equation (2.1) is thus integrated in the original phase space by treating the convective term (particle free motion) and the acceleration term in $\boldsymbol{p}$ separately.

It must be mentioned that the time-splitting scheme, first introduced by Cheng & Knorr (Reference Cheng and Knorr1976), presents shortcomings that prevent its direct application for three-dimensional advection in the momentum space $\boldsymbol{p}$ when strong relativistic effects are taken into account (see Huot et al. Reference Huot, Ghizzo, Bertrand, Sonnendrucker and Coulaud2003 for more details).

The VLEM solver uses a domain decomposition only in the $(x,y)$ configuration space plane with the local 1-D cubic spline method introduced by Crouseilles, Latu & Sonnendrucker (Reference Crouseilles, Latu and Sonnendrucker2007) with adapted boundary conditions which make possible a reconstruction of $f$ on the global space domain even on the cell boundaries. For the treatment in the momentum space, VLEM uses a global 2-D or three-dimensional cubic B-spline advection in $\boldsymbol{p}$ , depending on the 1-D (VLEM1D2V) or 2-D (VLEM2D3V) version in the configuration space. Thus, we only decompose the spatial domain into patches (or cells), each patch being devoted to one (message passing interface, MPI). Since the free-advection term (i.e. $\boldsymbol{p}/(m\gamma )$ ) does not depend explicitly on the $x$ or $y$ variables, a splitting technique is also possible in $x$ and $y$ and only the 1-D local spline interpolation technique is necessary in that case. One cell computes its own local cubic spline coefficients by solving linear systems, and Hermite boundary conditions are imposed on all the cells to reconstruct a smooth numerical solution. The domain decomposition in space is also used to solve the Maxwell equations in a self- consistent way with the Vlasov equation. Thus, the full kinetic and relativistic VLEM solver uses a fully parallelised hybrid MPI/OpenMP approach.

3.1. The multi-stream model in one-dimensional geometry

3.1.1. Basic equations

Let us consider a plasma model consisting of a collection of electrons embedded in a 1-D box with a fixed homogeneous neutralising ion background. Although the multi-stream model with one spatial dimension has already been presented by Bertrand et al. (Reference Bertrand, Albrecht-Marc, Reveille and Ghizzo2005) and Ghizzo & Bertrand (Reference Ghizzo and Bertrand2013), it is worth recalling its main features, particularly regarding its extension to higher dimensions.

For this purpose, let us consider the Hamiltonian of a single particle in an electromagnetic field $(\boldsymbol{E},\boldsymbol{B})$ in a 1-D configuration space, specifically in the $x$ direction. Here the 1-D geometry refers only to space, whereas the velocity, momentum and vector potential are three-dimensional. The Hamiltonian takes the form

(3.1) \begin{equation} H=mc^{2}\left [\sqrt {1+\frac {\left (\boldsymbol{P}_{c}-e\boldsymbol{A}\right )^{2}}{m^{2}c^{2}}}-1\right ]+e\phi \left (x,t\right)\!, \end{equation}

where $\phi$ is the electric potential. In (3.1), $\boldsymbol{P}_{c}=\boldsymbol{p}+e\boldsymbol{A}$ is the canonical momentum which depends on the particle momentum $\boldsymbol{p}=m\gamma \boldsymbol{v}$ .

In the Coulomb gauge $\mathrm{div}\boldsymbol{A}=0$ , we can express the vector potential as $\boldsymbol{A}=\boldsymbol{A}_{\perp }(x,t)=(0,A_{y},A_{z})$ . The perpendicular component of the electron’s canonical momentum is here conserved, $\boldsymbol{P}_{c\perp }=\boldsymbol{p}_{\perp }+e\boldsymbol{A}_{\perp }=\text{const}$ . We can express the Hamilton equations in the transverse directions as follows:

(3.2) \begin{equation} \frac {{\rm d}\boldsymbol{P}_{c\perp }}{{\rm d}t}=-\boldsymbol{\nabla }_{\perp }H=0, \end{equation}

while in the longitudinal direction, the equation takes the form

(3.3) \begin{equation} \frac {{\rm d}P_{cx}}{{\rm d}t}=-\frac {\partial H}{\partial x}. \end{equation}

From now on we use $\boldsymbol{C\!}_{j}$ to denote ${\boldsymbol{P}}_{c\perp }$ of a group of particles having the same value of the conserved canonical momentum component:

(3.4) \begin{equation} \boldsymbol{P}_{c\perp }=\boldsymbol{p}_{\perp }+e\boldsymbol{A}_{\perp }={\rm const.}=\boldsymbol{C\!}_{j}. \end{equation}

This conservation law (3.4) allows us to sample any distribution function in the momentum space associated with the translation invariance with a superposition of Dirac $\delta$ distributions, representing different ‘streams’ of particles, each ‘stream’ fulfilling the Vlasov equation (2.1). Thus, the full one-particle distribution function $f(x,\boldsymbol{p},t)$ can be written as a sum of Dirac $\delta$ distributions:

(3.5) \begin{equation} f\!\left (x,p_{x},\boldsymbol{p}_{\perp },t\right )=\sum _{j=-N_{s}}^{N_{s}}\alpha _{j}F_{j}\left (x,p_{x},t\right )\delta \left [\boldsymbol{p}_{\perp }-\left (\boldsymbol{C\!}_{j}-e\boldsymbol{A}_{\perp }\right )\right]\!. \end{equation}

As we will see later, in the generalisation of the model to two spatial dimensions, the parameter $\alpha _{j}$ acts as a normalisation factor for the distribution function, which depends on the nature of the equilibrium distribution under consideration. This allows us to assert that the $j{\rm th}$ distribution $F_{j}$ is normalised to one, i.e. $\int {\rm d}x{\rm d}p_{x}F_{j}(x,p_{x},t)=1$ . By introducing the change of variable $f_{j}(x,p_{x},t)=\alpha _{j}F_{j}(x,p_{x},t)$ , the distribution function $f_{j}$ is now normalised by the coefficient $\alpha _{j}$ and thus satisfies the normalisation condition $\int {\rm d}x{\rm d}p_{x}f_{j}(x,p_{x},t)=\alpha _{j}$ . The global normalisation of the distribution function $f$ imposes the condition $\sum _{j=1,N_\textit{stream}}\alpha _{j}=1$ among the different coefficients $\alpha _{j}$ .

We can now consider a plasma where the electrons are divided in $N_\textit{stream}=2N_{s}+1$ ‘bunches’ or ‘streams’ of particles. Each stream $j$ (for $j = - N_s, - N_s +1, \ldots, 0, 1, \ldots, N_s$ ) has the same initial perpendicular canonical momentum $\boldsymbol{C\!}_{j}$ and its dynamics is described by a reduced Hamiltonian given by

(3.6) \begin{equation} H_{j}\left (x,p_{x},t\right )=mc^{2}\left [\gamma _{j}\left (x,p_{x},t\right )-1\right ]+e\phi \left (x,t\right)\!, \end{equation}

with the corresponding associated Lorentz factor

(3.7) \begin{equation} \gamma _{j}(x,p_{x},t)=\sqrt {1+\frac {p_{x}^{2}}{m^{2}c^{2}}+\frac {(\boldsymbol{C\!}_{j}-e\boldsymbol{A}_{\perp }(x,t))^{2}}{m^{2}c^{2}}}. \end{equation}

Thus the multi-stream model can be reformulated as a summation of $2N_{s}+1$ Vlasov equations describing the dynamics of each cold electron stream. We obtain for each stream $j$

(3.8) \begin{equation} \frac {{\rm d}f_{j}}{{\rm d}t}=\frac {\partial f_{j}}{\partial t}+\frac {\partial H_{j}}{\partial p_{x}}\frac {\partial f_{j}}{\partial x}-\frac {\partial H_{j}}{\partial x}\frac {\partial f_{j}}{\partial p_{x}}=0, \end{equation}

whose derivation using (3.6) and (3.4) leads to the ‘reduced’ Vlasov equations given by

(3.9) \begin{equation} \frac {\partial f_{j}}{\partial t}+\frac {p_{x}}{m\gamma _{j}}\frac {\partial f_{j}}{\partial x}+\left [eE_{x}-\frac {1}{2m\gamma _{j}}\frac {\partial }{\partial x}((\boldsymbol{C\!}_{j}-e\boldsymbol{A}_{\perp }(x,t))^{2})\right ]\frac {\partial f_{j}}{\partial p_{x}}=0. \end{equation}

The set of these $2N_{s}+1$ Vlasov equations is then self-consistently coupled with the Maxwell equations. Let us discuss how this task can be accomplished.

In a 1-D spatial configuration, it is possible to separate the electric field into two parts, i.e. $\boldsymbol{E}=E_{x}\boldsymbol{e}_{x}+\boldsymbol{E}_{\perp }$ , where the longitudinal contribution $E_{x}=-{(\partial \phi}/{\partial x)}$ is a pure electrostatic field which obeys the Maxwell–Gauss equation and the transverse contribution is given by $\boldsymbol{E}_{\perp }=-{(\partial \boldsymbol{A}_{\perp }}/{\partial t)}$ . In the absence of any external magnetic field, $\boldsymbol{B}$ is purely perpendicular and is now given by $\boldsymbol{B}_{\perp }=\mathrm{rot}\boldsymbol{A}_{\perp }$ . The Maxwell equations $\mathrm{rot}\boldsymbol{E}=-{(\partial \boldsymbol{B}}/{\partial t)}$ and $\mathrm{div}\boldsymbol{B}=0$ are satisfied in a natural way, since we have $\mathrm{rot}(-{(\partial \boldsymbol{A}_{\perp }}/{\partial t)})=-({\partial}/{\partial t})(\mathrm{rot}\boldsymbol{A}_{\perp })$ and $\mathrm{div}(\mathrm{rot}\boldsymbol{A}_{\perp })=0$ (using the relation $\boldsymbol{B}=\mathrm{rot}\boldsymbol{A}$ ).

This allows us to simplify the procedure with which the Vlasov equations (3.9) are coupled with the Maxwell equations (2.3). Equations (3.9) are solved in a self-consistent way with the Maxwell–Gauss equation given by

(3.10) \begin{equation} \frac {\partial E_{x}}{\partial x}=\frac {e}{\varepsilon _{0}}\left [\sum _{j=-N_{s}}^{+N_{s}}n_{j}\left (x,t\right )-n_{0}\right]\!, \end{equation}

where $n_{j}$ is the density of the $j{\rm th}$ electron ‘stream’ determined by the relation

(3.11) \begin{equation} n_{j}\left (x,t\right )=\int _{-\infty }^{+\infty }{\rm d}p_{x}f_{j}\left (x,p_{x},t\right)\!. \end{equation}

The Maxwell–Ampère equation $\mathrm{rot}\boldsymbol{B}_{\perp }=\mu _{0}(\boldsymbol{J}_{\perp }+\varepsilon _{0}{(\partial \boldsymbol{E}_{\perp }}/{\partial t)})$ is written only for the perpendicular direction and can be transformed into a wave equation for the vector potential $\boldsymbol{A}_{\perp }$ in the form

(3.12) \begin{equation} \frac {\partial ^{2}\boldsymbol{A}_{\perp }}{\partial t^{2}}-c^{2}\frac {\partial ^{2}\boldsymbol{A}_{\perp }}{\partial x^{2}}=\frac {1}{\varepsilon _{0}}\sum _{j=-N_{s}}^{+N_{s}}\boldsymbol{J}_{\perp j}\left (x,t\right)\!, \end{equation}

where the source term is given by

(3.13) \begin{equation} \boldsymbol{J}_{\perp j}=\frac {e}{m}\left (\boldsymbol{C\!}_{j}-e\boldsymbol{A}_{\perp }\right )\int _{-\infty }^{+\infty }{\rm d}p_{x}\frac {f_{j}}{\gamma _{j}} \end{equation}

Finally, the reduced Vlasov equations given by (3.9) are self-consistently coupled to the Maxwell–Gauss equation (3.10) and the wave equation (3.12). For each population $j$ , the source terms are defined by (3.11) and (3.13).

In this case, the longitudinal component of the Lorentz force is expressed as $F_{Lx}=- ({1}/{2m\gamma _{j}})({\partial }/{\partial x})((\boldsymbol{C\!}_{j}-e\boldsymbol{A}_{\perp }(x,t))^{2})$ , which can be written as

(3.14) \begin{equation} F_{Lx}=\frac {e}{m\gamma _{j}}\left [\left (\boldsymbol{C\!}_{j}-e\boldsymbol{A}_{\perp }\right ).\frac {\partial \boldsymbol{A}_{\perp }}{\partial x}\right ]=e\left (\frac {\boldsymbol{p}_{\perp }}{m\gamma _{j}}\times \boldsymbol{B}_{\perp }\right )_{x}. \end{equation}

This expression involves only the perpendicular momentum of the electron.

From a numerical standpoint, it is feasible to solve the Maxwell equations presented in (2.3) self-consistently alongside the Vlasov equations. This can be achieved using the standard Yee (Reference Yee1966) technique combined with a leapfrog scheme for solving the Maxwell equations. This approach is facilitated by the potential decoupling between the electrostatic and electromagnetic components of the electric field. The transverse vector potential $\boldsymbol{A}_{\perp }$ can be readily computed from the relation

(3.15) \begin{equation} \frac {\partial \boldsymbol{A}_{\perp }}{\partial t}=-\boldsymbol{E}_{\perp }, \end{equation}

which in turn allows for the determination of the Lorentz factor $\gamma _{j}$ via (3.7).

The multi-stream model presents a compelling framework that enables us to succinctly express the invariance of the generalised canonical momentum in the perpendicular direction. Due to the fundamentally nonlinear nature of the model, it holds greater general significance than its applications to the study of Weibel-type instabilities. For example, it can provide a novel approach to determine stationary Vlasov–Maxwell solutions in the relativistic regime. While the overall formalism was introduced by Ghizzo & Bertrand (Reference Ghizzo and Bertrand2013) with a specific focus on the linear regime of Weibel-type instability, a comprehensive analysis of the saturation mechanisms associated with Weibel-type instabilities is detailed in Ghizzo (Reference Ghizzo2013a ). Further features associated to the non-propagative pure Weibel case are discussed in Del Sarto, Ghizzo & Sarrat (Reference Del Sarto, Ghizzo and Sarrat2024).

3.1.2. The case of a single stream

We now focus on the case of a single stream. Because this model, restricted to a 1-D (i.e. 2-D phase-space) system and with a single ‘stream’, is the ‘building block’ of a 2-D multi-stream approach, it is interesting to analyse it and its numerical integration scheme in more detail, so to underscore the advantages and constraints of the method. The Hamiltonian of a single ‘bunch’ of particles with equal energy and canonical momentum takes the form $H=mc^{2}\gamma (x,p_{x},t)+e\phi (x,t)$ and the corresponding reduced Vlasov equation reads

(3.16) \begin{equation} \frac {\partial f}{\partial t}+\frac {p_{x}}{m\gamma }\frac {\partial f}{\partial x}+\left [eE_{x}-\frac {mc^{2}}{2\gamma }\frac {\partial }{\partial x}\left (a^{2}\right )\right ]\frac {\partial f}{\partial p_{x}}=0, \end{equation}

where the Lorentz factor is $\gamma (x,p_{x},t)=\sqrt {1+(p_{x}^{2}/m^{2}c^{2})+a^{2}(x,t)}$ and $a=(e|\boldsymbol{A}_{\perp }|/mc)$ is the normalised amplitude of the perpendicular potential vector $\boldsymbol{A}_{\perp }$ of an electromagnetic wave. It appears that the direct integration using a 2-D full advection of the reduced relativistic Vlasov equation, given by (3.16), is necessary to obtain accurate solutions. It can be easily shown that (3.16) can take a conservative form:

(3.17) \begin{equation} \frac {\partial f}{\partial t}+\frac {\partial }{\partial x}\left (\frac {p_{x}}{m\gamma }f\right )+\frac {\partial }{\partial p_{x}}\left [\left (eE_{x}-\frac {mc^{2}}{2\gamma }\frac {\partial }{\partial x}\left (a^{2}\right )\right )f\right ]=0, \end{equation}

since we have $(\partial/\partial x)(p_{x}/m\gamma (x,p_{x},t))+(\partial/\partial p_{x})(eE_{x}- (mc^{2}/2\gamma (x,p_{x},t))$ $(\partial/\partial x)(a^{2}))=0$ .

Thus, the modification of the Lorentz factor, induced by the Hamiltonian reduction technique, imposes a spatial dependence on the Lorentz factor. This prevents the use of a time-splitting scheme as expressed in (3.17) and as is generally used in ‘Eulerian’ schemes of the Vlasov equation integration. From a numerical point of view this can be penalising as far as parallelisation is concerned. As we shall see in the more general case of a 2-D approach, we can overcome this issue by making a judicious use of the time-splitting technique, provided that each elementary ‘brick’ is kept in the form of a 2-D advection.

Equation (3.17) is completed by the Gauss equation:

(3.18) \begin{equation} \frac {\partial E_{x}}{\partial x}=\frac {e}{\epsilon _{0}}\left [n(x,t)-n_{0}\right]\!, \end{equation}

where

(3.19) \begin{equation} n(x,t)=\int _{-\infty }^{+\infty }f(x,p_{x}){\rm d}p_{x}. \end{equation}

In a 1-D system, the Maxwell equations (2.3) take a transverse form and a leapfrog scheme can be used to solve the magnetic $\boldsymbol{B}$ and electric $\boldsymbol{E}$ components of field in an alternate way. In the Yee (Reference Yee1966) method, the spatial derivatives of field components are evaluated using the simple two-point centred difference scheme. In principle, the derivatives in time are also discretised with second-order centred difference approximations, but the updates of $\boldsymbol{B}$ and $\boldsymbol{E}$ are staggered in time by one quarter time step. Thus we have for the perpendicular electric field:

(3.20) \begin{equation} \boldsymbol{E}_{\perp i}^{n+1/4}=\boldsymbol{E}_{\perp i}^{n-1/4}+\frac {c^{2}\Delta t}{2}\left (\mathrm{rot}\boldsymbol{B}_{\perp }^{n}\right )_{i}-\frac {\Delta t\boldsymbol{J}_{\perp i}^{n}}{2\varepsilon _{0}}, \end{equation}

where the source term is given by

(3.21) \begin{equation} \mathbf{J}_{\bot }(x,t)=-\frac {e^{2}\mathbf{A}_{\bot }}{m_{e}}\int _{-\infty }^{+\infty }\frac {f}{\gamma }{\rm d}p_{x}. \end{equation}

The discretised version of the magnetic field components reads as

(3.22) \begin{equation} \boldsymbol{B}_{\perp i+1/2}^{n+1/2}=\boldsymbol{B}_{\perp i+1/2}^{n}-\frac {\Delta t}{2}\left (\mathrm{rot}\boldsymbol{E}_{\perp }^{n+1/4}\right )_{i+1/2}, \end{equation}

to calculate the vector potential component $\boldsymbol{A}_{\perp }$ leading to

(3.23) \begin{equation} \boldsymbol{A}_{\perp i}^{n+1/2}=\boldsymbol{A}_{\perp i}^{n}-\frac {\Delta t}{2}\boldsymbol{E}_{\perp i}^{n+1/4}. \end{equation}

Numerical details for integrating the ‘rot’ operator are summarised in § 3.2.4.

3.1.3. The semi-Lagrangian technique

While the Eulerian approach consists of discretising the phase space on a fixed grid and applying either finite differences or finite volume or spectral method schemes, the semi-Lagrangian approach uses also the characteristics of the Vlasov equation, i.e. the concept of particle trajectory. This method of integration on an Eulerian grid is referred to here as the ‘Vlasov’ method by comparison with the Lagrangian approach (PIC method) although both methods solve the Vlasov equation.

Let us first recall the principles of the semi-Lagrangian method. By introducing the notation $\boldsymbol{\tilde {x}}=(x,p_{x})$ and $\boldsymbol{U}=((p_{x}/m\gamma),eE_{x}-(mc^{2}/2\gamma)\,(\partial/\partial x)\,(a^{2}))$ for the phase-space coordinates and advection field, the characteristics of (3.16) are solutions of the dynamical system

(3.24) \begin{equation} \frac {{\rm d}\boldsymbol{\tilde {X}}}{{\rm d}t}=\boldsymbol{U}(\boldsymbol{\tilde {X}}(t),t), \end{equation}

where we denote by $\boldsymbol{\tilde {X}}(t;\boldsymbol{\tilde {x}},t_{s})$ the solution at time $t$ whose value is $\tilde {\boldsymbol{x}}$ at time $t_{s}$ . We can now suppress the tilde notation in (3.24) to simplify the presentation. Taking $\boldsymbol{X}(t)$ a solution of (3.24), we have

(3.25) \begin{equation} \frac {{\rm d}}{{\rm d}t}\left (f\!\left (\boldsymbol{X}\!\left (t\right)\!,t\right )\right )=\frac {\partial f}{\partial t}+\frac {{\rm d}\boldsymbol{X}}{{\rm d}t}.\frac {\partial f}{\partial \boldsymbol{x}}=\frac {\partial f}{\partial t}+\boldsymbol{U}\left (\boldsymbol{X}\!\left (t\right )\!, t\right ).\frac {\partial f}{\partial \boldsymbol{x}}=0, \end{equation}

which indicates that $f$ is constant along the characteristics; hence,

(3.26) \begin{equation} f\!\left (\boldsymbol{X}\!\left (t;\boldsymbol{x},t_{s}\right)\!, t\right )=f\!\left (\boldsymbol{X}\!\left (t_{s};\boldsymbol{x},t_{s}\right)\!,t_{s}\right )=f\!\left (\boldsymbol{x},t_{s}\right)\!, \end{equation}

for any times $t$ and $t_{s}$ and phase-space coordinates $\boldsymbol{x}$ . Thus, given the value of the distribution function $f$ at the mesh points, at any given time step, we obtain the new value at mesh point $\boldsymbol{x}_{i}$ using the fact that

(3.27) \begin{equation} f^{*}=f\!\left (\boldsymbol{x}_{i},t_{n+1}\right )=f\!\left (\boldsymbol{X}\!\left (t_{n};\boldsymbol{x}_{i},t_{n+1}\right)\!,t_{n}\right)\!. \end{equation}

For each mesh point $\boldsymbol{x}_{i}$ , $f$ is computed in two steps.

1. Step1: Find the starting point of the characteristic ending at $\boldsymbol{x}_{i}$ , i.e. $\boldsymbol{X}(t_{n};\boldsymbol{x}_{i},t_{n}+\triangle t)$ ; $\triangle t$ is the time step used in the numerical scheme.

2. Step2: Compute $f(\boldsymbol{X}(t_{n};\boldsymbol{x}_{i},t_{n}+\triangle t),t_{n})$ by interpolation, $f$ being known only on the mesh points at time $t_{n}=n\triangle t$ .

For step 1, the starting point of the characteristic is obtained to second-order accuracy by

(3.28) \begin{equation} \frac {\boldsymbol{X}\!\left (t_{n+1}\right )-\boldsymbol{X}\!\left (t_{n}\right )}{\triangle t}=\frac {\boldsymbol{x}_{i}-\boldsymbol{X}\!\left (t_{n}\right )}{\triangle t}=\boldsymbol{U}\left (\boldsymbol{X}\!\left (t_{n+ 1/2}\right)\!,t_{n+ 1/2}\right)\!. \end{equation}

Writing $\boldsymbol{X}(t_{n+(1/2)})={[\boldsymbol{X}(t_{n}+\Delta t)+\boldsymbol{X}(t_{n})]}/{2}$ , there exists $\boldsymbol{\psi }_{i}$ such that the quantity $\boldsymbol{X}(t_{n})$ can be written in the form $\boldsymbol{X}(t_{n})=\boldsymbol{x}_{i}-\boldsymbol{\psi }_{i}$ and $\boldsymbol{X}(t_{n}+\triangle t)=\boldsymbol{x}_{i}$ , which allows us to write (3.28) in the following form:

(3.29) \begin{equation} \boldsymbol{\psi }_{i}=\triangle t\boldsymbol{U}\left (\boldsymbol{x}_{i}-\frac {\boldsymbol{\psi }_{i}}{2},t_{n+ 1/2}\right)\!. \end{equation}

This can be solved iteratively for the unknown $\boldsymbol{\psi }_{i}$ by writing

(3.30) \begin{equation} \boldsymbol{\psi }_{i}^{\left (k+1\right )}=\triangle t\boldsymbol{U}\left (\boldsymbol{x}_{i}-\frac {\boldsymbol{\psi }_{i}^{\left (k\right )}}{2},t_{n+ 1/2}\right)\!. \end{equation}

Once $\boldsymbol{\psi }_{i}$ are known, $f(\boldsymbol{x}_{i}-\boldsymbol{\psi }_{i},t_{n})$ is interpolated by a tensor product of cubic B-splines.

Figure 1.

Illustration of the numerical scheme used in the 1-D multi-stream code: the different shifts or advections for solving the set of Vlasov equations (bottom frame) are represented by black arrows; the top and middle frames correspond to the computation of the electric field and to the magnetic field together with the vector potential.

The overall scheme is shown in figure 1: the distribution function $f^{n+1}$ , calculated at time $t_{n+1}=(n+1)\triangle t$ from the data of $f$ at time $t_{n}$ , is duplicated; a second distribution $f'^{n+ 1/2}$ has been introduced with an initial shift of half a time step. In this way, the two distributions $f$ and $f'$ are alternately advected. To advect the function $f$ , the Poisson equation (3.18) is solved using data from $f'^{n+ 1/2}$ , which is itself advected in turn using the longitudinal electric field $E_{x}^{n+1}$ calculated at time $t_{n+1},$ from data from $f^{n+1}$ . This process, together with the iterative sequence (3.30), enables both 2-D advections to be solved to the second order in time. The transverse electric field is calculated successively at times $t_{n+ 1/4}$ and $t_{n+ 3/4}$ , alternately with the magnetic field $\boldsymbol{B}_{\perp }$ and vector potential $\boldsymbol{A}_{\perp }$ at times $t_{n+ 1/2}$ and $t_{n+1}$ .

The scheme used for a single stream can be easily extended to a number $N_\textit{stream}$ of streams, increasing the level of code parallelism since the various streams can be distributed over different MPI processes. The reconstruction of $f$ is made using the characteristics in the backward direction in time, which is equivalent to considering, at each time step, a new set of ‘particles’ located on the grid and to determining the characteristic feet and then interpolating the function using the grid points in the immediate neighbourhood.

3.2. The extension to a two-dimensional system

The Hamiltonian reduction is a procedure that ‘eliminates’ a redundant degree of freedom from a Hamiltonian system by providing an equivalent description in a smaller-dimensional phase space. In essence, it is equivalent to the reduction of the configuration space in analytical mechanics, allowed by the use of holonomic constraints. This procedure is based on the homogeneity of space in a given direction – say $z$ – i.e. a translational invariance in space along $z$ for a Hamiltonian of the kind $H=mc^{2}\gamma +e\phi$ , and it leads to the invariance of momentum according to the Noether theorem. It results in the invariance of the generalised canonical momentum $P_{cz}=p_{z}+eA_{z}(x,y,t)=\text{const.}$ through the Hamilton equation

(3.31) \begin{equation} \frac {{\rm d}P_{cz}}{{\rm d}t}=-\frac {\partial H}{\partial z}=0. \end{equation}

In this case the reduced description is characterised by a subgroup of $2N_{s}+1$ ‘streams’, each stream being described by the Hamiltonian $H_{j}$ . It is the symmetry of the global Hamiltonian (here associated to the lacking $z$ variable) that motivates the choice of the reduced variables, the so-called distribution of the generalised canonical momentum $C_{j}$ . Thus the general form of the distribution function may be represented in the following form (where $\boldsymbol{x}=(x,y,0)$ ):

(3.32) \begin{equation} f\!\left (\boldsymbol{x,p},t\right )=\sum _{j=-N_{s}}^{N_{s}}f_{j}\left (\boldsymbol{x},\boldsymbol{p}_{\parallel },t\right )\delta \left [p_{z}-\left (C_{j}-eA_{z}\left (x,y,t\right )\right )\right]\!. \end{equation}

This Hamiltonian reduction technique is not merely an approximation of the physical system under study; it provides a description of the system’s dynamics based on a class of exact solutions that arises from a particular symmetry of the system. This technique reduces the dimensionality of the phase space due to the conservation of the generalised canonical momentum. Furthermore in the numerical integration, the existence of an analytical solution enhances the convergence of the global solution by providing a better description of the current density $J_{z}$ in the $z$ direction, as it is constructed from the multi-stream representation of the distribution function given in (3.32). Finally, the model can be applied to any type of distribution function, including relativistic distribution cases, thus allowing for the reintroduction of a notion of temperature anisotropy, which is generally impossible with a Maxwell–Jüttner distribution.

3.2.1. Basic equations

By writing the Hamiltonian in the standard form

(3.33) \begin{equation} H=mc^{2}\sqrt {1+\frac {\left (\boldsymbol{P}_{c}-e\boldsymbol{A}\right ).\left (\boldsymbol{P}_{c}-e\boldsymbol{A}\right )}{m^{2}c^{2}}}+e\phi \left (x,y,t\right)\!, \end{equation}

where $\boldsymbol{P}_{c}=\boldsymbol{p}+e\boldsymbol{A}$ is the canonical momentum, i.e. the conjugate variable to the configuration space $\boldsymbol{x}=(x,y,0)$ , it is possible to reduce the dimensionality of the phase space. In the Coulomb gauge ( $\mathrm{div}\boldsymbol{A}=0$ ), we have $\boldsymbol{A}=(A_{x},A_{y},A_{z})=\boldsymbol{A}_{\parallel }+\boldsymbol{A}_{\perp }$ , with $\boldsymbol{A}_{\parallel }=(A_{x},A_{y},0)$ , $\boldsymbol{A}_{\perp }=(0,0,A_{z})$ . Since $\boldsymbol{A}$ does not depend on the $z$ variable, we have $(\partial A_{x}/\partial x)+{(\partial A_{y}}/{\partial y)}=0$ in the Coulomb gauge. The magnetic field $\boldsymbol{B}=\mathrm{rot}\boldsymbol{A}$ satisfies the following set of equations:

(3.34) \begin{equation} B_{x}=\frac {\partial A_{z}}{\partial y},\quad B_{y}=-\frac {\partial A_{z}}{\partial x}\quad \mathrm{and}\quad B_{z}=\frac {\partial A_{y}}{\partial x}-\frac {\partial A_{x}}{\partial y}, \end{equation}

allowing one to separate the different contributions of the transverse magnetic (TM) $(E_{x},E_{y},B_{z})$ and transverse electric (TE) $(E_{z},B_{x},B_{y})$ modes. By denoting the in-plane (by symbol $\parallel$ ) and the $z$ -component (as the perpendicular or $\perp$ ) contributions of the canonical momentum, we have $\boldsymbol{P}_{c}=\boldsymbol{P}_{c\parallel }+\boldsymbol{P}_{c\perp }$ with $\boldsymbol{P}_{c\perp }=P_{cz}\boldsymbol{e}_{z}$ and $\boldsymbol{P}_{c\parallel }=P_{cx}\boldsymbol{e}_{x}+P_{cy}\boldsymbol{e}_{y}$ , the Hamilton equations read as

(3.35) \begin{equation} \frac {{\rm d}\boldsymbol{x}}{{\rm d}t}=\frac {\partial H}{\partial \boldsymbol{P}_{c}}\mathrm{\quad and}\quad \frac {{\rm d}\boldsymbol{P}_{c}}{{\rm d}t}=-\frac {\partial H}{\partial \boldsymbol{x}}. \end{equation}

Along the $z$ direction, the perpendicular direction here, we have

(3.36) \begin{equation} P_{cz}=p_{z}+eA_{z}\left (x,y,t\right )={\rm const.}\equiv C_{j}, \end{equation}

since $\mathrm{d}P_{cz}/\mathrm{d}t\equiv \dot {P_{cz}}=-\partial H/\partial z=0$ . Therefore, without loss of generality, we can consider a plasma where particles are divided into $N_\textit{stream}=2N_{s}+1$ ‘streams’, each stream $j$ (for $j\in [-N_{s},N_{s}]$ ) having the initial perpendicular momentum $C_{j}$ . We can now define for each particle population $j$ , a reduced Hamiltonian denoted $H_{j}$ which satisfies the definition

(3.37) \begin{equation} H_{j}=mc^{2}\gamma _{j}+e\phi \left (x,y,t\right)\!, \end{equation}

where the Lorentz factor $\gamma _{j}$ , associated with the particle stream $j$ , is defined by

(3.38) \begin{equation} \gamma _{j}=\sqrt {1+\frac {\left (\boldsymbol{P}_{c\parallel }-e\boldsymbol{A}_{\parallel }\right ).\left (\boldsymbol{P}_{c\parallel }-e\boldsymbol{A}_{\parallel }\right )}{m^{2}c^{2}}+\frac {\left (C_{j}-eA_{z}\left (x,y,t\right )\right )^{2}}{m^{2}c^{2}}}. \end{equation}

Denoting $\boldsymbol{p}_{\parallel }=(p_{x},p_{y},0)$ and $\boldsymbol{x}=(x,y,0)$ , for each particle stream $j$ , the corresponding Hamilton equations read as

(3.39) \begin{equation} \frac {{\rm d}\boldsymbol{x}}{{\rm d}t}=\frac {\partial H_{j}}{\partial \boldsymbol{P}_{c\parallel }}=\frac {\boldsymbol{P}_{c\parallel }-e\boldsymbol{A}_{\parallel }}{m\gamma _{j}}=\frac {\boldsymbol{p}_{\parallel }}{m\gamma _{j}}, \end{equation}

(3.40) \begin{align} \frac {{\rm d}\boldsymbol{P}_{c\parallel }}{{\rm d}t}&=-\boldsymbol{\nabla }H_{j}=-\frac {1}{2m\gamma _{j}}\boldsymbol{\nabla }\left (\boldsymbol{P}_{c\parallel }-e\boldsymbol{A}_{\parallel }\right )^{2}-e\boldsymbol{\nabla }\phi -\frac {1}{2m\gamma _{j}}\boldsymbol{\nabla }\left (C_{j}-eA_{z}\right )^{2}\nonumber \\ &=\frac {e}{m\gamma _{j}}\left [\boldsymbol{p}_{\parallel }.\boldsymbol{\nabla }\boldsymbol{A}_{\parallel }+\boldsymbol{p}_{\parallel }\times \left (\boldsymbol{\nabla }\times \boldsymbol{A}_{\parallel }\right )\right ]-e\boldsymbol{\nabla }\phi -\frac {1}{2m\gamma _{j}}\boldsymbol{\nabla }\left (C_{j}-eA_{z}\right )^{2}, \end{align}

where $\boldsymbol{\nabla }=\partial _{x}\boldsymbol{e}_{x}+\partial _{y}\boldsymbol{e}_{y}$ . By noting that $\mathrm{d}\boldsymbol{P}_{c\parallel }/\mathrm{d}t={\rm d}\boldsymbol{p}_{\parallel }/dt+e\partial \boldsymbol{A}_{\parallel }/\partial t+e\boldsymbol{p}_{\parallel }.\boldsymbol{\nabla }\boldsymbol{A}_{\parallel }/$ $m\gamma _{j}$ , (3.40) reads as

(3.41) \begin{equation} \frac {{\rm d}\boldsymbol{p}_{\parallel }}{{\rm d}t}+\frac {e\partial \boldsymbol{A}_{\parallel }}{\partial t}=\frac {e}{m\gamma _{j}}\boldsymbol{p}_{\parallel }\times B_{z}\boldsymbol{e}_{z}-\frac {1}{2m\gamma _{j}}\boldsymbol{\nabla }\left (C_{j}-eA_{z}\right )^{2}-e\boldsymbol{\nabla }\phi . \end{equation}

Finally, the Hamilton equations become

(3.42) \begin{equation} \frac {{\rm d}\boldsymbol{x}}{{\rm d}t}=\frac {\boldsymbol{p}_{\parallel }}{m\gamma _{j}}, \end{equation}

(3.43) \begin{equation} \frac {{\rm d}\boldsymbol{p}_{\parallel }}{{\rm d}t}=e\boldsymbol{E}_{\parallel }+\frac {e}{m\gamma _{j}}\boldsymbol{p}_{\parallel }\times B_{z}\boldsymbol{e}_{z}-\frac {1}{2m\gamma _{j}}\boldsymbol{\nabla }\left (C_{j}-eA_{z}\left (x,y,t\right )\right )^{2}, \end{equation}

where we have used the definition $\boldsymbol{E}_{\parallel }=-\boldsymbol{\nabla }\phi -\partial \boldsymbol{A}_{\parallel }/\partial t$ . We can now define, for each particle population $j$ , a distribution function denoted $f_{j}(\boldsymbol{x},\boldsymbol{p}_{\parallel },t)$ which satisfies the following relativistic Vlasov equation:

(3.44) \begin{align} \frac {{\rm d}f_{j}}{{\rm d}t}=\frac {\partial f_{j}}{\partial t}+\frac {{\rm d}\boldsymbol{x}}{{\rm d}t}\boldsymbol{\cdot}\boldsymbol{\nabla }f_{j}+\frac {{\rm d}\boldsymbol{p}_{\parallel }}{{\rm d}t}\boldsymbol{\cdot}\frac {\partial f_{j}}{\partial \boldsymbol{p}_{\parallel }}=0, \end{align}

or equivalently,

(3.45) \begin{align} &\frac {\partial f_{j}}{\partial t}+\frac {p_{x}}{m\gamma _{j}}\frac {\partial f_{j}}{\partial x}+\frac {p_{y}}{m\gamma _{j}}\frac {\partial f_{j}}{\partial y}+\left [eE_{x}+{\frac {ep_{y}B_{z}}{m\gamma _{j}}}-\frac {1}{2m\gamma _{j}}\partial _{x}\left (C_{j}-eA_{z}\right )^{2}\right ]\frac {\partial f_{j}}{\partial p_{x}}\nonumber \\ &\qquad +\left [eE_{y}-{\frac {ep_{x}B_{z}}{m\gamma _{j}}}-\frac {1}{2m\gamma _{j}}\partial _{\mathrm{y}}\left (C_{j}-eA_{z}\right )^{2}\right ]\frac {\partial f_{j}}{\partial p_{y}}=0. \end{align}

We have used the invariance property $P_{cz}=p_{z}-eA_{z}(x,y,t)=\text{const.}=C_{j}$ in the perpendicular direction to the plane $(x,y)$ . Equation (3.45) is coupled to

(3.46) \begin{equation} \frac {\partial A_{z}}{\partial t}=-E_{z}. \end{equation}

For each population of particles corresponding to the ‘bunch’ $j$ , the source terms used in the Maxwell equations are defined by the following quantities:

(3.47) \begin{equation} \boldsymbol{J}_{\parallel ,j}=e\int {\rm d}^{2}p\frac {\boldsymbol{p}_{\parallel }}{m\gamma _{j}}f_{j} ,\quad J_{z,j}=e\left (C_{j}-eA_{z}\right )\int {\rm d}^{2}p\frac {f_{j}}{m\gamma _{j}}. \end{equation}

Finally, the global source terms are

(3.48) \begin{equation} \boldsymbol{J}_{\parallel }=\sum _{j=-N_{s}}^{+N_{s}}\boldsymbol{J}_{\parallel ,j}\quad \mathrm{and}\quad J_{z}=\sum _{j=-N_{s}}^{+N_{s}}J_{z,j}. \end{equation}

At this stage, several points are worth highlighting:

1. (i) The relation (3.46) is not an approximation, since the electric potential does not depend on the $z$ variable.

2. (ii) In (3.45) the Lorentz factor $\gamma _{j}$ becomes dependent both on the particle momentum and on the space variables, rendering obsolete the use of a standard time-splitting technique by separating the global solution of the Vlasov equation into four 1-D equations.

3. (iii) The Hamiltonian reduction technique can also be used in a Lagrangian scheme (PIC codes).

3.2.2. The numerical scheme with two dimensions in configuration space

In applied mathematics, the Strang splitting method (Strang Reference Strang1968) is a numerical technique used to solve differential equations that can be decomposed into a sum of differential operators. This method can be extended to tackle multi-dimensional partial differential equations by reducing them to a sum of lower-dimensional problems. The time-splitting technique was first introduced by Cheng & Knorr (Reference Cheng and Knorr1976) for solving the Vlasov–Poisson system. More recently, following the work of Morrison (Reference Morrison1980), Crouseilles, Einkemmer & Faou (Reference Crouseilles, Einkemmer and Faou2015) proposed a new splitting method based on the decomposition of the Hamiltonian functional in the following form:

(3.49) \begin{align} \mathcal{H}\left [\boldsymbol{\chi }\right ]&=\frac {1}{2}\epsilon _{0}\int \left |\boldsymbol{E}\right |^{2}{\rm d}^{3}x+\frac {1}{2}\epsilon _{0}c^{2}\int \left |\boldsymbol{B}\right |^{2}{\rm d}^{3}x+\frac {1}{2}m\int \left |v\right |^{2}\!f\!\left (\boldsymbol{x},\boldsymbol{v},t\right ){\rm d}^{3}x{\rm d}^{3}v\nonumber\\&=\mathcal{H}_{E}+\mathcal{H}_{B}+\mathcal{H}_{f}. \end{align}

This formulation relates to the Vlasov–Maxwell system in the non-relativistic regime, where $\mathcal{H}_{E},\mathcal{H}_{B}$ and $\mathcal{H}_{f}$ represent the contributions of electric, magnetic and kinetic energy, respectively. Here, $\boldsymbol{\chi }=(f,\boldsymbol{E},\boldsymbol{B})$ is a vector dependent on the distribution function $f$ and on the electromagnetic components $\boldsymbol{E}$ and $\boldsymbol{B}$ . In the Hamiltonian formalism, the different components $\chi _{i}$ satisfy the equation $\partial \chi _{i}/\partial t=\{ \chi _{i},\mathcal{H}\}$ , where $\{ \boldsymbol{\cdot}\}$ is the Poisson bracket. From the relation $\mathcal{H}=\mathcal{H}_{E}+\mathcal{H}_{B}+\mathcal{H}_{f}$ one can construct an approximate solution in the form

(3.50) \begin{equation} \boldsymbol{\chi }\left (t\right )=\exp \left (\frac {\mathcal{H}_{E}t}{2}\right )\exp \left (\frac {\mathcal{H}_{B}t}{2}\right )\exp \left (\mathcal{H}_{f}t\right )\exp \left (\frac {\mathcal{H}_{B}t}{2}\right )\exp \left (\frac {\mathcal{H}_{E}t}{2}\right )\boldsymbol{\chi }\left (0\right)\!. \end{equation}

This expression represents the solution obtained through the Strang splitting scheme, which is a second-order-accurate method. The primary challenge when attempting to cast the Vlasov–Maxwell equations into Hamiltonian form lies in identifying a suitable set of canonical variables. This task is complicated by the fact that the components $\chi _{i}$ are non-canonical and that the Poisson bracket takes on a complex form. The multi-stream technique provides a more straightforward approach as it utilises canonical variables.

In the multi-stream approach, a similar Strang splitting technique can be applied to solve the system consisting of $2N_{s}+1$ reduced Vlasov equations (3.45). Thus, the Vlasov equations (3.45) can be split in the following sequence, where each resolution step preserves the exact conservation of $f_{j}$ :

(3.51) \begin{equation} \mathrm{step}\:\mathcal{X}\left (t\right )\left [f_{j}\right ]:\;\frac {\partial f_{j}}{\partial t}+\frac {p_{x}}{m\gamma _{j}}\frac {\partial f_{j}}{\partial x}-\frac {1}{2m\gamma _{j}}\partial _{x}\left (C_{j}-eA_{z}\right )^{2}\frac {\partial f_{j}}{\partial p_{x}}=0, \end{equation}

(3.52) \begin{equation} \mathrm{step}\:\mathcal{Y}\left (t\right )\left [f_{j}\right ]:\;\frac {\partial f_{j}}{\partial t}+\frac {p_{y}}{m\gamma _{j}}\frac {\partial f_{j}}{\partial y}-\frac {1}{2m\gamma _{j}}\partial _{y}\left (C_{j}-eA_{z}\right )^{2}\frac {\partial f_{j}}{\partial p_{y}}=0, \end{equation}

together with

(3.53) \begin{equation} \mathrm{step}\:\mathcal{R}\left (t\right )\left [f_{j}\right ]:\quad \frac {\partial f_{j}}{\partial t}+e\left (E_{x}+\frac {p_{y}B_{z}}{m\gamma _{j}}\right )\frac {\partial f_{j}}{\partial p_{x}}+e\left (E_{y}-\frac {p_{x}B_{z}}{m\gamma _{j}}\right )\frac {\partial f_{j}}{\partial p_{y}}=0, \end{equation}

where the Lorentz factor is given by

(3.54) \begin{equation} \gamma _{j}=\sqrt {1+\frac {p_{x}^{2}}{m^{2}c^{2}}+\frac {p_{y}^{2}}{m^{2}c^{2}}+\frac {\left (C_{j}-eA_{z}\left (x,y,t\right )\right )^{2}}{m^{2}c^{2}}}. \end{equation}

Taken separately, the previous steps, linked to the $\mathcal{X},\mathcal{Y}$ and $\mathcal{R}$ operators, can be exactly integrated along their respective characteristics, using 2-D B-spline interpolation techniques. Thus each step, denoted by the $\mathcal{X},\mathcal{Y}$ and $\mathcal{R}$ operators, allows one to exactly conserve the mass since the corresponding advections are made globally in the considered phase-space variables (see § 3.2.3). Another advantage of the method is the possibility of treating large 2-D (global) advections without limitation on the time step. Indeed, in comparison with the semi-Lagrangian scheme used in the fully kinetic code VLEM2D3V, the advections employed in the multi-stream code are now global and require the global inversion of a tridiagonal matrix. In the VLEM2D3V code, the spatial advections are local and necessitate a Courant–Friedrichs–Lewy condition on the time step, which is not the case here.

At each step $\mathcal{X}$ or $\mathcal{Y}$ , a transposition is applied between advected and non-advected variables (a domain decomposition may then be applied for ‘non-advected’ variables). However, in the full kinetic VLEM2D3V version, local spline interpolations are used, so transpositions are not necessary.

Thus a formal solution of (3.51)–(3.53) can be obtained in the form of a sequence of successive ‘shifts’ of advections $\mathcal{X}(\Delta t/2),\mathcal{Y}(\Delta t/2)$ made over half a time step $\triangle t/2$ :

(3.55) \begin{equation} \mathrm{step}\:\mathcal{X}\left (\Delta t/2\right ):\left \{ \begin{array}{l} \mathrm{transpose}\;f^{n}\left (\overline {x},\overline {y},p_{x},p_{y}\right )\Rightarrow f_{t}^{n}\left (x,\overline {y},p_{x},\overline {p_{y}}\right )\!,\\[5pt] \mathrm{compute}\:f_{t}^{*}\left (x,\overline {y},p_{x},\overline {p_{y}}\right )={{f_{t}^{n}}}\left (x-\dfrac {p_{x}}{m\gamma _{j}}\dfrac{\Delta t}{2},\overline {y},p_{x}\right.\\[10pt] \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\,\,\,\left.+\dfrac {\Delta t}{4m\gamma _{j}}\partial _{x}\left (C_{j}-eA_{z}\right )^{2},\overline {p_{y}}\right )\!,\\[10pt] \mathrm{transpose}\;f_{t}^{*}\left (x,\overline {y},p_{x},\overline {p_{y}}\right )\Rightarrow f^{*}\left (\overline {x},\overline {y},p_{x},p_{y}\right)\!. \end{array}\right . \end{equation}

At the end of this step, performed by advancing the system over half a time step $\triangle t/2$ , a global 2-D advection in the $\tilde {x}=(x,p_{x})$ phase space is required. A matrix transposition is then performed between the variables $(x,p_{x})$ and the spatial coordinates $(x,y)$ so as to enable calculation of the electromagnetic field components at the end of each sequence; here the electromagnetic field uses a domain decomposition in configuration space. The 2-D advection in $(x,p_{x})$ is performed via a global inversion of the matrix with respect to the $(x,p_{x})$ coordinates, the domain decomposition being performed in the perpendicular directions $\tilde {y}=(y,p_{y})$ . The bars introduced in the notation refer to the domain decomposition used for parallelism. This procedure is reproduced for all $\mathcal{X}(\Delta t/2)$ , $\mathcal{Y}(\Delta t/2)$ shifts, and the $\mathcal{R}(\Delta t)$ rotation to restore the domain decomposition:

(3.56) \begin{equation} \mathrm{step}\:\mathcal{Y}\left (\Delta t/2\right ):\left \{ \begin{array}{l} \mathrm{transpose}\;{{{f}^{{*}}}}\left (\overline {x},\overline {y},p_{x},p_{y}\right )\Rightarrow {{f}_{{t_{bis}}}^{{*}}}\left (\overline {x},y,\overline {p_{x}},p_{y}\right )\!,\\[8pt]\mathrm{compute}\:{{f_{t_{bis}}^{**}}}\left (\overline {x},y,\overline {p_{x}},p_{y}\right )={{f}_{{t_{bis}}}^{{*}}}\left (\overline {x},y-\dfrac {p_{y}}{m\gamma _{j}}\dfrac {\Delta t}{2},\overline {p_{x}},p_{y}\right.\\[10pt]\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\left.+\dfrac {\Delta t}{4m\gamma _{j}}\partial _{y}\left (C_{j}-eA_{z}\right )^{2}\right )\!,\\[10pt] \mathrm{transpose}\;{{f}_{{t_{bis}}}^{{**}}}\left (\overline {x},y,\overline {p_{x}},p_{y}\right )\Rightarrow f^{**}\left (\overline {x},\overline {y},p_{x},p_{y}\right)\!. \end{array}\right . \end{equation}

This is followed by the $\mathcal{R}$ rotation integration over a full time step $\Delta t$ :

(3.57) \begin{equation} \mathrm{step}\:\mathcal{R}\left (\Delta t\right ):\mathrm{compute}\;f^{***}\left (\overline {x},\overline {y},\boldsymbol{p}\right )=f^{**}\left (\boldsymbol{P}^{n+1}\left (\overline {x},\overline {y},\boldsymbol{p}\right )\right )\!, \end{equation}

where a Boris (Reference Boris1970) method is used to solve the step $\mathcal{R}(\Delta t)$ by separating the electric and magnetic contribution of the electromagnetic field, and by rewriting the characteristics of (3.53) in the form

(3.58) \begin{align} \boldsymbol{v}^{+}-\boldsymbol{v}^{-}=\frac {e}{2m}\left (\boldsymbol{v}^{+}+\boldsymbol{v}^{-}\right )\times B_{z}\boldsymbol{e}_{z}, \end{align}

where $\boldsymbol{v}^{n}=\boldsymbol{p}^{n}/(m\gamma _{j}^{n})=\boldsymbol{v}^{-}-e\boldsymbol{E}_{\parallel }\Delta t/(2m)$ and $\boldsymbol{v}^{n+1}\equiv \boldsymbol{p}^{n+1}/(m\gamma _{j}^{n+1})=\boldsymbol{v}^{+}+e\boldsymbol{E}_{\parallel }{}\Delta t/(2m)$ . The shifts of the distribution function $f_{j}$ , over half a time step, along the directions $x{-}p_{x}$ and $y{-}p_{y}$ , are denoted with $\mathcal{X}(\Delta t/2)$ , $\mathcal{Y}(\Delta t/2)$ , respectively, and the ‘rotation’ in the momentum space over a full time step is indicated as $\mathcal{R}(\Delta t)$ .

Several remarks must be made:

1. (i) A straightforward way to distribute the computational work involved across processors is to assign a part of the $x{-}y$ grid to each processor. Thus, the code is initially parallelised by using a regular spatial domain decomposition in the $x$ and $y$ directions, which is used when solving the Maxwell equations and the rotation step $\mathcal{R}$ for the Vlasov equation. Transpositions are then required to adapt this domain decomposition for each shift in $x{-}p_{x}$ (denoted $\mathcal{X}$ ) and $y{-}p_{y}$ (denoted $\mathcal{Y}$ ).

2. (ii) Transpositions of the global matrix are required before and after each $\mathcal{X}(\Delta t/2)$ or $\mathcal{Y}(\Delta t/2)$ shift. The transposed quantities are denoted by $\overline {x},\overline {y},\overline {p_{x}}$ and $\overline {p_{y}}$ in (3.55)–(3.57), i.e. the phase-space variables for which a domain decomposition has been applied.

3. (iii) Information shared between the threads of neighbouring subdomains combines OpenMP and MPI. The numerical scheme parallelises computed-bounds loop with OpenMP on a local processus, sharing communications between threads, and inter-processor communications are handled by MPI.

4. (iv) Another level of parallelisation is also introduced in the multi-stream model by distributing the different ‘streams’ on different MPI processes. According to the MPI terminology, non-blocking communications were used, point-to-point communications can be performed with non blocking semantics as well, having resources free for computation.

5. (v) The global scheme $f_{j}^{n+1}=\mathcal{X}(\Delta t/2)\circ \mathcal{Y}(\Delta t/2)\circ \mathcal{\mathcal{R}}(\Delta t)\circ \mathcal{Y}(\Delta t/2)\circ \mathcal{X}(\Delta t/2)f_{j}^{n}$ , illustrated in figure 2, exhibits second-order accuracy in time (see Appendix A). Thus the semi-Lagrangian model used in the multi-stream approach can be seen as a generalisation of the well-known Cheng & Knorr (Reference Cheng and Knorr1976) technique, based on the use of the time-splitting scheme and non-local cubic spline interpolation of the distribution function. The operator $\circ$ acts as a convolution product of two operators $\mathcal{A}$ and $\mathcal{B}$ , which applies to the distribution function $f(\boldsymbol{x},\boldsymbol{p},t)$ between two given times, defined by $f(t_{2})[\boldsymbol{x,p}]=\mathcal{A}\circ \mathcal{B}f(t_{1})=\mathcal{A}(\mathcal{B}f(t_{1})[\boldsymbol{x,p}])$ ; the operators $\mathcal{A}$ and $\mathcal{B}$ being chosen from the operators $\mathcal{X}$ , $\mathcal{Y}$ and the rotation $\mathcal{R}$ . The second-order accuracy results from the fact that the time stencil is centred, which provides a symmetric self-adjoint numerical integrator.

Figure 2.

Illustration of the numerical scheme used in the multi-stream code in two dimensions in configuration space: the different shifts or advections for solving the set of Vlasov equations (bottom frame) are represented by black arrows; the top and middle frames correspond to the computation of the electric field and to the magnetic field together with the vector potential.

3.2.3. Conservation of the global scheme

The Vlasov equations (3.45) can be rewritten in a global advective form:

(3.59) \begin{equation} \frac {\partial f_{j}}{\partial t}+\boldsymbol{U}\left (\boldsymbol{X},t\right )\boldsymbol{\cdot}\frac {\partial f_{j}}{\partial \boldsymbol{X}}=0, \end{equation}

where $\boldsymbol{X}=(x,y,p_{x},p_{y})$ denotes the generalised phase-space coordinate and

(3.60) \begin{align} \boldsymbol{U}=\left (\begin{array}{c} \dfrac {p_{x}}{m\gamma _{j}}\\[14pt] \dfrac {p_{y}}{m\gamma _{j}}\\[14pt] eE_{x}+e\dfrac {p_{y}B_{z}}{m\gamma _{j}}-\dfrac {1}{2m\gamma _{j}}\partial _{x}\left (C_{j}-eA_{z}\right )^{2}\\[14pt] eE_{y}-e\dfrac {p_{x}B_{z}}{m\gamma _{j}}-\dfrac {1}{2m\gamma _{j}}\partial _{y}\left (C_{j}-eA_{z}\right )^{2} \end{array}\right ) \end{align}

the generalised advection field. Since $\boldsymbol{U}$ is divergence free, we can put (3.59) in a conservative form:

(3.61) \begin{equation} \frac {\partial f_{j}}{\partial t}+\frac {\partial }{\partial \boldsymbol{X}}\boldsymbol{\cdot}\left [\boldsymbol{U}\left (\boldsymbol{X},t\right )f_{j}\right ]=0. \end{equation}

By now separating each 2-D advection, (3.61) reads

(3.62) \begin{equation} \frac {\partial f_{j}}{\partial t}+\frac {\partial }{\partial \boldsymbol{\boldsymbol{X}_{x}}}\boldsymbol{\cdot}\left [\boldsymbol{U}_{x}\left (\boldsymbol{\boldsymbol{X}_{}},t\right )f_{j}\right ]+\frac {\partial }{\partial \boldsymbol{\boldsymbol{X}_{y}}}\boldsymbol{\cdot}\left [\boldsymbol{\boldsymbol{U}_{y}}\left (\boldsymbol{\boldsymbol{X}_{}},t\right )f_{j}\right ]+\frac {\partial }{\partial \boldsymbol{\boldsymbol{X}_{p}}}\boldsymbol{\cdot}\left [\boldsymbol{\mathbf{\boldsymbol{U}_{p}}}\left (\boldsymbol{X}_{},t\right )f_{j}\right ]=0, \end{equation}

where

(3.63) \begin{equation} \boldsymbol{\boldsymbol{U}_{x}}=\left (\begin{array}{c} \dfrac {p_{x}}{m\gamma _{j}}\\[14pt] 0\\[3pt] -\dfrac {1}{2m\gamma _{j}}\partial _{x}\left (C_{j}-eA_{z}\right )^{2}\\[14pt] 0 \end{array}\right )\quad \mathrm{and}\quad \boldsymbol{X_{x}}=\left (\begin{array}{c} x\\[3pt] 0\\[3pt] p_{x}\\[3pt] 0 \end{array}\right ), \end{equation}

(3.64) \begin{equation} \boldsymbol{\boldsymbol{U}_{y}}=\left (\begin{array}{c} 0\\[3pt] \dfrac {p_{y}}{m\gamma _{j}}\\[14pt] 0\\[3pt] -\dfrac {1}{2m\gamma _{j}}\partial _{y}\left (C_{j}-eA_{z}\right )^{2} \end{array}\right )\quad \mathrm{and}\quad \boldsymbol{X_{y}}=\left (\begin{array}{c} 0\\[3pt] y\\[3pt] 0\\[3pt] p_{y} \end{array}\right ), \end{equation}

and finally

(3.65) \begin{equation} \boldsymbol{U_{p}}=\left (\begin{array}{c} 0\\[3pt] 0\\[3pt] e\left (E_{x}+\dfrac {p_{y}B_{z}}{m\gamma _{j}}\right )\\[12pt] e\left (E_{y}-\dfrac {p_{x}B_{z}}{m\gamma _{j}}\right ) \end{array}\right )\quad\mathrm{and}\quad\boldsymbol{X_{p}}=\left (\begin{array}{c} 0\\[3pt] 0\\[2pt] p_{x}\\[2pt] p_{y} \end{array}\right ). \end{equation}

It is well known that we can solve separately (3.62) with keeping second-order accuracy (see Marchuk Reference Marchuk1982) by solving the following sequence of conservative equations:

(3.66) \begin{equation} \frac {\partial f_{j}}{\partial t}+\frac {\partial }{\partial \boldsymbol{X_{x}}}\boldsymbol{\cdot}\left [\boldsymbol{U_{x}}\left (\boldsymbol{X_{}},t\right )f_{j}\right ]=0, \end{equation}

(3.67) \begin{equation} \frac {\partial f_{j}}{\partial t}+\frac {\partial }{\partial \boldsymbol{X_{y}}}\boldsymbol{\cdot}\left [\boldsymbol{U_{y}}\left (\boldsymbol{X_{}},t\right )f_{j}\right ]=0, \end{equation}

(3.68) \begin{equation} \frac {\partial f_{j}}{\partial t}+\frac {\partial }{\partial \boldsymbol{X_{p}}}\boldsymbol{\cdot}\left [\boldsymbol{U_{p}}\left (\boldsymbol{X_{}},t\right )f_{j}\right ]=0, \end{equation}

in a symmetric time-centred way, i.e. (3.66) and (3.67) are integrated over half a time step, while (3.68) is solved over a full time step. We repeat the sequence of resolution of (3.67) and (3.66) over half a time step. As each advection field $\boldsymbol{U_{x}},\boldsymbol{U_{y}}$ and $\boldsymbol{U_{p}}$ is divergence free, we can recover the sequence of advective equations given in (3.51)–(3.53).

3.2.4. Resolution of the Maxwell equations

Finally, the reduced Vlasov equations given by (3.45) for $-N_{s}\leqslant j\leqslant +N_{s}$ (with a total number of $2N_{s}+1$ ‘streams’) are self-consistently coupled to the Maxwell equations for the two sets of TM and TE components of the electromagnetic field:

(3.69) \begin{equation} \frac {\partial E_{x}}{\partial t}=c^{2}\frac {\partial B_{z}}{\partial y}-\frac {J_{x}}{\varepsilon _{0}}, \end{equation}

(3.70) \begin{equation} \frac {\partial E_{y}}{\partial t}=-c^{2}\frac {\partial B_{z}}{\partial x}-\frac {J_{y}}{\varepsilon _{0}}, \end{equation}

(3.71) \begin{equation} \frac {\partial B_{z}}{\partial t}=\frac {\partial E_{x}}{\partial y}-\frac {\partial E_{y}}{\partial x} \end{equation}

and

(3.72) \begin{equation} \frac {\partial E_{z}}{\partial t}=c^{2}\left (\frac {\partial B_{y}}{\partial x}-\frac {\partial B_{x}}{\partial y}\right )-\frac {J_{z}}{\varepsilon _{0}}, \end{equation}

(3.73) \begin{equation} \frac {\partial B_{x}}{\partial t}=-\frac {\partial E_{z}}{\partial y}, \end{equation}

(3.74) \begin{equation} \frac {\partial B_{y}}{\partial t}=\frac {\partial E_{z}}{\partial x}. \end{equation}

We have used here the well-known Yee (Reference Yee1966) scheme to solve the Maxwell equations where the TE and TM modes are decoupled. In order to ensure high speed and accuracy, it is convenient to use a domain decomposition in the configuration $(x,y)$ grid and a leapfrog scheme for the integration of the electromagnetic field in time. The leapfrog scheme is thus applied to the space coordinates too, which means that electric and magnetic fields are shifted in time by $\Delta t/4$ and some components of them are also shifted by $\Delta x/2$ and $\Delta y/2$ , as illustrated in the top and middle frames in figure 2. Thus, the electric components are defined at times $t_{n-1/4}$ , $t_{n+1/4}$ and $t_{n+3/4}$ (as a usual subcycling technique used in PIC codes), while the magnetic and $A_{z}$ components are defined at times $t_{n}$ , $t_{n+1/2}$ and $t_{n+1}$ , and the current density $\boldsymbol{J}$ is computed from (3.48) at times $t_{n}$ and $t_{n+1/2}$ . As a result, the discretised expressions of (3.69)–(3.71), in Cartesian coordinates, can be written as

(3.75) \begin{equation} E_{x\,i+\frac {1}{2},j}^{n+\frac {1}{4}}=E_{x,i+\frac {1}{2},j}^{n-\frac {1}{4}}+\frac {c^{2}\Delta t}{2\Delta y}\left (B_{z\,i+\frac {1}{2},j+\frac {1}{2}}^{n}-B_{z\,i+\frac {1}{2},j-\frac {1}{2}}^{n}\right )-\frac {\Delta t}{2\varepsilon _{0}}J_{x\,i+\frac {1}{2},j}^{n}, \end{equation}

(3.76) \begin{equation} E_{y\,i,j+\frac {1}{2}}^{n+\frac {1}{4}}=E_{y,i,j+\frac {1}{2}}^{n-\frac {1}{4}}+\frac {c^{2}\Delta t}{2\Delta x}\left (B_{z\,i+\frac {1}{2},j+\frac {1}{2}}^{n}-B_{z\,i-\frac {1}{2},j+\frac {1}{2}}^{n}\right )-\frac {\Delta t}{2\varepsilon _{0}}J_{y\,i,j+\frac {1}{2}}^{n}, \end{equation}

(3.77) \begin{equation} B_{z\,i+\frac {1}{2},j+\frac {1}{2}}^{n+\frac {1}{2}}=B_{z\,i+\frac {1}{2},j+\frac {1}{2}}^{n}-\frac {\Delta t}{2\Delta x}\left (E_{y\,i+1,j+\frac {1}{2}}^{n+\frac {1}{4}}-E_{y\,i,j+\frac {1}{2}}^{n+\frac {1}{4}}\right )+\frac {\Delta t}{2\Delta y}\left (E_{x\,i+\frac {1}{2},j+1}^{n+\frac {1}{4}}-E_{x\,i+\frac {1}{2},j}^{n+\frac {1}{4}}\right)\!. \end{equation}

The discretised versions of (3.72)–(3.74) are

(3.78) \begin{equation} B_{x\,i,j+\frac {1}{2}}^{n+\frac {1}{2}}=B_{x\,i,j+\frac {1}{2}}^{n}-\frac {\Delta t}{2\Delta y}\left (E_{z\,i,j+1}^{n+\frac {1}{4}}-E_{z\,i,j}^{n+\frac {1}{4}}\right )\!, \end{equation}

(3.79) \begin{equation} B_{y\,i+\frac {1}{2},j}^{n+\frac {1}{2}}=B_{y\,i+\frac {1}{2},j}^{n}+\frac {\Delta t}{2\Delta x}\left (E_{z\,i+1,j}^{n+\frac {1}{4}}-E_{z\,i,j}^{n+\frac {1}{4}}\right )\!, \end{equation}

(3.80) \begin{equation} E_{z\,i,j}^{n+\frac {1}{4}}=E_{z,i,j}^{n-\frac {1}{4}}+\frac {c^{2}\Delta t}{2\Delta x}\left (B_{y\,i+\frac {1}{2},j}^{n}-B_{y\,i-\frac {1}{2},j}^{n}\right )-\frac {c^{2}\Delta t}{2\Delta y}\left (B_{x\,i,j+\frac {1}{2}}^{n}-B_{x\,i,j-\frac {1}{2}}^{n}\right )-\frac {\Delta t}{2\varepsilon _{0}}J_{z\,i,j}^{n}, \end{equation}

together with

(3.81) \begin{equation} J_{z\,i,j}^{n}=e\sum _{l=-N_{s}}^{+N_{s}}\left (C_{l}-eA_{z\,i,j}^{n}\right )\int \frac {f_{l}^{n}}{m\gamma _{l}}{\rm d}p_{x}{\rm d}p_{y}\quad \mathrm{and}\quad A_{z\,i,j}^{n+\frac {1}{2}}=A_{z,i,j}^{n}-\frac {\Delta t}{2}E_{z\,i,j}^{n+\frac {1}{4}}. \end{equation}

The sequence for the integration of the Maxwell equations is repeated two times, in order to compute the values of the advections terms $(p_{x}/(m\gamma _{j}),\partial _{x}(C_{j}-eA_{z})^{2}/(2m\gamma _{j}))$ and $(p_{y}/(m\gamma _{j}),\partial _{y}(C_{j}-eA_{z})^{2}/(2m\gamma _{j}))$ in (3.51), (3.52) at times $t_{n+1/4}$ and $t_{n+3/4}$ required for solving the $x-p_{x}$ and $y-p_{y}$ 2-D advections, respectively. Finally the 2-D advection in $p_{x}-p_{y}$ can also be performed after computing the advection field $\boldsymbol{E}_{\parallel }^{*}=(\boldsymbol{E}_{\parallel }^{n+1/4}+\boldsymbol{E}_{\parallel }^{n+3/4})/2$ and using the data of $B_{z}^{n+1/2}$ . As indicated in figure 2, the current densities $\boldsymbol{J}_{\parallel }^{n}$ and $\boldsymbol{J}_{\parallel }^{n+1/2}$ are computed at time $t_{n}=n\Delta t$ and at the end of the first sequence $\mathcal{X}(\Delta t/2)\circ \mathcal{Y}(\Delta t/2)$ , respectively.

3.3. Relevance of the multi-stream approach to systems beyond classical Vlasov plasmas

We conclude this section by noting that the Hamiltonian reduction implemented with the multi-stream description may find broader applications than classical plasmas. In principle, it involves any Hamiltonian system described by Vlasov- or Liouville-like equations, and displaying the appropriate symmetry (i.e. the presence of a cyclic coordinate). Two examples that we discuss here are the non-relativistic Vlasov–Maxwell system including spin effects and the covariant formulation of the relativistic Vlasov model.

First, the scalar non-relativistic spin Vlasov–Maxwell system is described in the extended phase space $(\boldsymbol{x},\boldsymbol{p},\boldsymbol{s})$ for an electron distribution function $\mathcal{\mathit{\mathcal{G}}}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{s},t)$ . As shown in Marklund & Morrison (Reference Marklund and Morrison2011), the previous Hamiltonian formulation can be viewed as a classical pre-quantisation property where a new distribution is introduced in the form $\mathcal{F}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{\hat {s}},t)=\mathcal{G}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{s},t)\delta (|\boldsymbol{s}|-(\hbar/2))$ , where it obeys the same kind of Vlasov equation but with a reduced dimensionality in phase space. This is the consequence of the existence of Casimir invariants $\int {\rm d}^{9}z\mathcal{C}(\mathcal{G},s^{2})={\rm const.}$ and the integral reduces from integration over the global phase space element ${\rm d}^{9}z={\rm d}^{3}x{\rm d}^{3}p{\rm d}^{3}s$ to ${\rm d}^{3}x{\rm d}^{3}p{\rm d}^{2}s$ with ${\rm d}^{2}s=\sin \theta _{s}{\rm d}\theta _{s}{\rm d}\varphi _{s}$ and $\boldsymbol{\hat {s}}$ is the unit spin vector. Second, the Vlasov equation can also be extended to a more general treatment, i.e. in a covariant description expressed in terms of the phase-space coordinates $(x^{0}=ct,\boldsymbol{x};p^{0}=\mathcal{E}/c,\boldsymbol{p})$ (Marsden et al. Reference Marsden, Montgomery, Morrison and Thompson1986; Liboff Reference Liboff2003; Cary & Brizard Reference Cary and Brizard2009). Furthermore, a variational formulation for the Vlasov–Maxwell system has been presented in Brizard (Reference Brizard2000) that uses constrained variations for the Vlasov distribution in an eight-dimensional phase space.

Here we adopt the general notation $\boldsymbol{x}$ for the three-dimensional configuration space. Here $x^{\mu }=(x^{0},\boldsymbol{x})$ denotes the space–time position of a particle and $p^{\mu }=(p^{0},\boldsymbol{p})$ its momentum–energy coordinates. In the case of a system composed of charged particles submitted to an electromagnetic field $F^{\mu \nu }$ , where $F^{\mu }=(e/m)p_{\nu }F^{\mu \nu }$ , the equation of motion for the four-momentum $p^{\mu }$ is

(3.82) \begin{equation} \frac {{\rm d}p^{\mu }}{{\rm d}t}=\frac {eF^{\mu \nu }}{m}p_{\nu }\equiv eF^{\mu \nu }u_{\nu }, \end{equation}

where $u^{\mu }=(u^{0}=\gamma c,\boldsymbol{u}=\gamma \boldsymbol{v})={\rm d}x^{\mu }/{\rm d}\tau$ is the four-velocity and ${\rm d}/{\rm d}\tau =\gamma {\rm d}/{\rm d}t$ is the derivative with the proper time. We use the Minkowski space–time metric $g^{\mu \nu }=\mathrm{diag}(-1,+1,+1,+1)$ . We have $p^{\mu }=mu^{\mu }$ and $F^{\mu \nu }=\partial ^{\mu }A^{\nu }-\partial ^{\nu }A^{\mu }$ denotes the Faraday tensor. The contravariant space–time derivatives are $\partial ^{\mu }=g^{\mu \nu }\partial _{\nu }=(-\partial /\partial x^{0},\boldsymbol{\nabla })$ and $A^{\mu }=(A^{0}=\phi /c,\boldsymbol{A})$ is the four-potential. The relativistic Vlasov equation is recovered as usual and reads

(3.83) \begin{equation} \frac {{\rm d}\mathcal{F}}{{\rm d}\tau }=u^{\mu }\frac {\partial \mathcal{F}}{\partial x^{\mu }}+F^{\mu }\left (x^{\mu },p^{\mu }\right )\frac {\partial \mathcal{F}}{\partial p^{\mu }}=0. \end{equation}

Since the number of particles in the system is assumed to be conserved, the eight-current in $\mu$ space must also be conserved and the continuity equation reads

(3.84) \begin{equation} \frac {\partial }{\partial x^{\mu }}\left [\frac {{\rm d}x^{\mu }}{{\rm d}\tau }\mathcal{F}\left (x^{\mu },p^{\mu }\right )\right ]+\frac {\partial }{\partial p^{\mu }}\left (\frac {{\rm d}p^{\mu }}{{\rm d}\tau }\mathcal{F}\right )=0. \end{equation}

Using the $(x^{\mu },p^{\mu })$ notation, with the usual definition of the momentum–energy four-vector for $p^{\mu }$ , the electron distribution function $\mathcal{F}=\mathcal{F}(x^{\mu },p^{\mu })$ satisfies the covariant relativistic Vlasov equation (3.83) which reads

(3.85) \begin{equation} p^{\mu }\frac {\partial \mathcal{F}}{\partial x^{\mu }}+e\frac {\partial \mathcal{F}}{\partial p^{\mu }}F^{\mu \beta }p_{\beta }=0. \end{equation}

In an eight-dimensional phase space, this scalar equation takes the explicit form

(3.86) \begin{equation} \frac {\partial \mathcal{F}}{\partial t}+\frac {\boldsymbol{p}}{m\gamma }\boldsymbol{\cdot}\frac {\partial \mathcal{F}}{\partial \boldsymbol{x}}+e\left (\boldsymbol{E}+\frac {\boldsymbol{p}}{m\gamma }\times \boldsymbol{B}\right )\boldsymbol{\cdot}\frac {\partial \mathcal{F}}{\partial \boldsymbol{p}}+e\left (\frac {\boldsymbol{p}}{m\gamma }\boldsymbol{\cdot}\boldsymbol{E}\right )\frac {\partial \mathcal{F}}{\partial p^{0}}=0. \end{equation}

The covariant relativistic equation (3.86) formally has an additional degree of freedom compared with the ‘classical’ expression of the Vlasov equation (2.1) and the motion of a charged particle is located on the surface $\mathcal{H}={(p^{\mu }p_{\mu })}/{(2m)}=-mc^{2}/2$ , that derives from the condition $p_{\mu }p^{\mu }=-m^{2}c^{2}$ , $\mathcal{H}$ being a Lorentz scalar (see Brizard & Chan (Reference Brizard and Chan1999) for more details). By considering a solution in the form

(3.87) \begin{equation} \mathcal{F}\left (x^{\mu },p^{\mu }\right )=\delta \left (p^{0}-m\gamma c\right )\frac {f\!\left (\boldsymbol{x},\boldsymbol{p},t\right )}{\gamma }, \end{equation}

where $f$ is the usual distribution written in the six-dimensional phase space $(\boldsymbol{x,p})$ and $\gamma$ the Lorentz factor, and after a $p^{0}$ -integration over (3.86), we recover the usual expression of the relativistic Vlasov equation given by (2.1).

It can be seen that the formal solution provided by (3.87) constitutes a class of exact solutions describing the physics of collisionless plasmas in a relativistic regime (Vlasov–Maxwell model) under the constraint $p_{\mu }p^{\mu }=-m^{2}c^{2}$ , which reflects the Einstein relation between the energy of an electron $\mathcal{E}$ and its momentum $\boldsymbol{p}$ , i.e. $p^{2}-({(\mathcal{E}}/{c)})^{2}=-m^{2}c^{2}$ .

We thus see that both the non-relativistic spin and the (relativistic) covariant formulation in the Vlasov description correspond to Hamiltonian reduction techniques based on the introduction of a Dirac distribution function which corresponds to the metric constraint $p^{\mu }p_{\mu }=-m^{2}c^{2}$ in the covariant formulation of the Vlasov equation or two spin degrees of freedom with a fixed $|\boldsymbol{s}|=\hbar/2$ . Thus such a Hamilton reduction procedure can be applied again in the presence of translation symmetry using the invariance of canonical momentum in a transverse direction.

4.1. Filamentation aspects in Lagrangian and Eulerian schemes

In order to understand the advantages and disadvantages of the numerical technique associated with the resolution of the reduced Vlasov–Maxwell system governed by (3.16) and (3.23), coupled in a self-consistent way with the Maxwell equations (3.20)–(3.22), we consider the case of the interaction of a high-intensity laser wave interacting with a plasma.

From a numerical point of view, the modelling of the Vlasov–Maxwell system relies mainly on two numerical techniques. The first, and most widely used, is the PIC method.

Particle codes involve a Lagrangian formulation: the PIC scheme can be regarded as a discretisation of phase space in terms of a superposition of moving elements, usually referred to as super-particles. In PIC codes, super-particle trajectories are computed from the electromagnetic fields, prescribed on a fixed grid whose typical size is of order of the Debye length. At the end of the time step, the charge of each ‘super-particle’ is redistributed among the neighbouring grid points, allowing one to solve the Maxwell equations. This process used in PIC codes (the redistribution of the charge of super-particles) involves a smoothing of the information which efficiently decreases the individual effects introduced by the grid. This makes the PIC code a very well adapted scheme to follow the kinematic filamentation of the distribution function in velocity space, one of the fundamental properties of the Vlasov equation. In the filamentation process, the information transfer from small to large wavenumbers follows from the energy transfer between scales and 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.

In Vlasov simulations that are performed on a (Eulerian) fixed grid in phase space, the grid inevitably becomes too coarse as the fine graining develops. Such a process leads to an entropy increase. Thus, this ‘filamentation’ problem is well known in Eulerian Vlasov simulations, while the PIC approach is not affected by this problem and the filamentation is usually very well described in a Lagrangian framework. However, the kinematic filamentation of $f$ can also interact with the filamentation induced by certain physical instabilities, leading to modifications in the dynamics of $f$ through a heating process in the presence of developed turbulence (Ghizzo, Del Sarto & Betar Reference Ghizzo, Del Sarto and Betar2023; Ghizzo & Del Sarto Reference Ghizzo and Del Sarto2023). We shall see that the multi-stream model can improve the numerical description of the filamentation in the semi-Lagrangian case, by introducing an exact solution of the distribution in the form of a Dirac comb, while preserving the ‘noiseless’ character of the code. In this way, the multi-stream approach could offer certain advantages for certain physical solutions where the filamentation of $f$ can play a primordial role, as can be the case in CFIs or in Weibel-type instabilities, both associated with the emergence and amplification of a magnetic field in the system.

On the other hand, the Eulerian or Lagrangian techniques use a statistical approach where the graininess parameter $g=1/(n_{0}\lambda _{D}^{3})$ tends towards zero. This characteristic is one of the undeniable advantages of the Eulerian or semi-Lagrangian methods, and is often referred to as a noiseless code (in the statistical sense). The re-introduction of a ‘numerical grain’ effect in the PIC code, associated with the finite size of the super-particles, leads to an amplification of the ‘thermal’ noise of the particles; this property is often associated with the increase in numerical ‘noise’ encountered in PIC codes. This process is obviously a major drawback in PIC codes, and a number of numerical techniques have been developed to reduce this noise without, however, succeeding in eliminating it completely. This problem makes it difficult to describe the tails of distribution functions in the PIC approach, whereas Eulerian and semi-Lagrangian techniques are less prone to this problem, provided sufficient phase-space sampling is used. Thus, the use of a multi-stream model can considerably reduce this problem through judicious selection of ‘streams’, by positioning two streams in these very low-density regions. This is an advantage that can in principle be used both in the semi-Lagrangian code and in a PIC approach.

4.2. The case of a single stream in relativistic laser–plasma interaction and kinematic filamentation aspects

Before making a systematic comparison between the reduced multi-stream model and the fully kinetic VLEM code, it is worth returning to the single-stream case, which models a cold plasma in the perpendicular direction. This example, although simplified, highlights some of the fundamental properties of the Vlasov equation, such as the filamentation of $f$ in the velocity space, and how the information is processed in the usual techniques for numerically solving the Vlasov equation. In the case of a single-stream model, (3.16) is used to solve the reduced Vlasov equation in phase space. This model constitutes the elementary ‘brick’ which is used in any Hamiltonian reduction approach in the case of a higher dimensionality of the phase space. The corresponding numerical scheme has been presented in § 3.1.3.

We focus then on the relativistic regime of laser–plasma interaction. A well-known solution is given by Akhiezer & Polovin (Reference Akhiezer and Polovin1956) for a circularly polarised electromagnetic wave in a homogeneous plasma, corresponding to $\mathbf{A}_{\bot }^{2}={\rm const.}$ , and the linear dispersion relation of a propagating electromagnetic wave of frequency $\omega _{0}$ and wavevector $k_{0}$ , is given by $\omega _{0}^{2}={(\omega _{p}^{2}}/{\gamma _{0})}+k_{0}^{2}c^{2}$ , where the Lorentz factor $\gamma _{0}$ is given by $\gamma _{0}^{2}=\sqrt {1+a_{0}^{2}}$ .

The simulation that we describe now has been performed using a ratio of the electron mean density to the critical density of $n_{0}/n_{c}=1$ ( $n_{c}$ being usually defined as $m\omega _{0}^{2}\varepsilon _{0}/e^{2}$ ) and a quiver momentum value of $a_{0}=eE_{0}/m\omega _{0}c=\sqrt {3}$ , which corresponds to a laser intensity of $I\lambda _{0}^{2}=8.2\times 10^{18}$ $\mathrm{\mathrm{W}\,cm^{-2}\, \mathrm{\mu}m^{2}}$ (where $\lambda _{0}$ is the pump wavelength). The initial electron temperature is $k_{\mathrm{B}}T_{e}=3\,\mathrm{keV}$ . The initial large-amplitude pump wave $(\omega _{0},k_{0})$ is described by the electric field component $\mathbf{E_{\perp }}=(0,E_{0}\cos k_{0}x,E_{0}\sin k_{0}x)$ . Choosing $k_{0}c/\omega _{p}=1/\sqrt {2}$ , we obtain $\omega _{0}=\omega _{p}$ , i.e. a ratio of $n_{0}/\gamma _{0}n_{c}=0.5$ , i.e. well below the relative plasma density. Here, we have used a phase-space sampling of $N_{x}N_{p_{x}}$ of $512\times 768$ points and a time step of $\triangle t\omega _{p}=0.01.$

The corresponding plasma wave predicted by theory is then very close to 2 $\Delta k$ . Here, we have solved the dispersion relation, i.e.

(4.1) \begin{equation} D_{+}D_{-}=\frac {\omega _{p}^{2}a_{0}^{2}}{2\gamma _{0}^{3}}\left (\frac {k^{2}c^{2}}{D_{\mathrm{p}}}-1\right )\left (D_{+}+D_{-}\right )\!, \end{equation}

where $D_{\mathrm{p}}$ and $D_{\pm }$ correspond, respectively, to the dispersion relation of the electron plasma wave and of the electromagnetic waves in the relativistic plasma. We have the usual relations:

(4.2) \begin{equation} D_{p}=\omega ^{2}-\frac {\omega _{p}^{2}}{\gamma _{0}};\quad D_{\pm }=\left (\omega \pm \omega _{0}\right )^{2}-\left (k\pm k_{0}\right )^{2}c^{2}-\frac {\omega _{p}^{2}}{\gamma _{0}}. \end{equation}

By solving the dispersion relation in the case of a cold plasma, the growth rate of the relativistic parametric instability of the most unstable mode is located at $kc/\omega _{p}=1.40$ and corresponds to a maximum growth rate of $\eta _{\mathrm{th}}/\omega _{p}\simeq 0.409$ . In figure 3 we have plotted, on top, the time evolution of the most unstable plasma mode (here mode 2) on a logarithmic scale: the curve indicates a growth rate of around $\eta _{\mathrm{num}}/\omega _{p}\simeq 0.406\pm 0.003$ in good agreement with the expected value predicted by the linear theory (for a cold plasma). But the more striking point is that we did not introduce any initial density perturbation in our simulation. The instability starts just from the round-off error, which for a 64-bit machine is of order $10^{-15}$ . This result clearly demonstrates that the low-noise character of the multi-stream code allows a very powerful and extremely precise study of the growth of instabilities over a large number of decades.

Figure 3.

Top: time evolution of the most unstable plasma mode (here the mode 2) on a logarithmic scale obtained by the multi-stream code with one ‘stream’. The curve exhibits a growth rate close to $\eta _{\mathrm{num}}/\omega _{p}\simeq 0.406$ in good agreement with the expected value predicted by the linear theory (for a cold plasma). Bottom: the phase-space representation of the electron distribution function at two different instants. The arrow indicates the first time close to the saturation where the distribution is plotted (at left).

The last point concerns the problem of the filamentation in $x{-}p_{x}$ phase space. To demonstrate the efficiency of our algorithm, we look now at the particle dynamics in phase space. Figure 3 (bottom frames) shows the behaviour of the electron distribution function in the $x{-}p_{x}$ phase space at two times taken in the nonlinear regime of the instability. Here, the formation of smaller and smaller filaments in phase space results from the physical process induced by the relativistic parametric instability, which is well recovered in the multi-stream approach.

4.3. The case of the current filamentation instability

The following example demonstrates the possibility of describing the CFI (Fried Reference Fried1959) with a limited number of streams. Here, two streams are sufficient for an accurate description of the nonlinear saturation process of a CFI instability produced by two counter-streaming electron (physical) beams.

Potential applications include laser–plasma interaction (Silva et al. Reference Silva, Fonseca, Tonge, Mori and Dawson2002; Tzoufras et al. Reference Tzoufras, Ren, Tsung, Tonge, Mori, Fior, Fonseca and Silva2006; Okada & Ogawa Reference Okada and Ogawa2007), astrophysics (Tautz & Schlickeiser, Reference Tautz and Schlickeiser2005a ,Reference Tautz and Schlickeiser b ; Bret et al. Reference Bret, Stocken, Fiuza, Ruyer, Gremillet, Narayan and Silva2013; Petersen, Glenzer & Fiuza Reference Petersen, Glenzer and Fiuza2021) and gamma-ray burst sources (Medvedev & Loeb Reference Medvedev and Loeb1999; Medvedev Reference Medvedev2006).

The initial distribution condition is composed of two Maxwellians with beam momenta centred at $C_{1}/mc=-0.9$ and $C_{2}/mc=0.9$ . Here, we focus attention on purely transverse initial perturbation (on the magnetic $B_{z}$ component) with wavevector of type $\boldsymbol{k}_{0}=k_{0}\boldsymbol{e}_{x}$ , i.e. along the $x$ direction (here the longitudinal direction), perpendicular to the two counter-streaming electron beams. We have chosen a symmetric case corresponding to beam densities of $n_{01}=n_{02}=0.5n_{0}$ . Only the fundamental mode (for the magnetic field) $k_{0}c/\omega _{p}=1$ is excited. The corresponding initial magnetic field is then given by ${eB_{z}}/{m\omega _{p}}={(eB_{0}}/{m\omega_{p})}\sin ({(k_{0}c}/{\omega _{p})}\, {(x\omega _{p}}/{c)})$ in normalised units, with an amplitude perturbation of $eB_{0}/m\omega _{p}=10^{-4}$ . Electron plasma temperature is $k_{\mathrm{B}}T_{e}=2\,\mathrm{keV}$ , in both the $p_{x}$ and $p_{y}$ directions in the VLEM1D2V version, while in the reduced model only the $p_{x}$ direction has a non-zero temperature. We keep the same temperature in the reduced model. The phase-space sampling used here in the VLEM1D2V version is $N_{x}N_{p_{x}}N_{p_{y}}$ $=256\times 257^{2}$ , i.e. $1.69\times 10^{7}$ grid points or ‘particles’. The time step used in both simulations is $\Delta t\omega _{p}=0.003$ . The phase-space sampling, used in the multi-stream simulation, is somewhat higher, with $N_{x}N_{p_{x}}$ corresponding to $513^{2}$ grid points.

Figure 4.

Time evolution of the normalised magnetic field component ${(eB_{z}}/{m\omega _{p})}(x={(L_{x}}/{2)},t)$ , in logarithmic scale. As expected, the magnetic field is amplified and the linear growth rate is found to be in good agreement with the theoretical value predicted by the linear dispersion relation of CFI.

First, in order to obtain an accurate estimation of the growth rate of the instability in the linear phase, we have solved, for the corresponding physical parameters of simulation, the linear dispersion relation obtained directly from the multi-stream model which reads as in the fluid approximation

(4.3) \begin{align} &\left (1-\sum _{j=1}^{N_{s}}\frac {\omega _{pj}^{2}}{\omega ^{2}\varGamma _{0j}}\right )\left (-\omega ^{2}+k^{2}c^{2}+\sum _{j=1}^{N_{s}}\frac {\omega _{pj}^{2}}{\varGamma _{0j}^{3}}+\frac {k^{2}c^{2}}{\omega ^{2}}\sum _{j=1}^{N_{s}}\frac {\omega _{pj}^{2}}{\varGamma _{0j}^{3}}\frac {C_{j}^{2}}{m^{2}c^{2}}\right )\nonumber\\&\quad\,\,\quad=-{\frac {k^{2}c^{2}}{\omega ^{4}}}\left (\sum _{j=1}^{N_{s}}\frac {\omega _{pj}^{2}}{\varGamma _{0j}^{2}}\frac {C_{j}}{mc}\right )^{2}, \end{align}

where $\varGamma _{0j}=\sqrt {1+{(C_{j}^{2}}/{m^{2}c^{2})}}$ . In figure 4, we present the time evolution of the component ${(eB_{z}}/{m\omega _{p})}(x={(L_{x}}/{2)},t)$ of the magnetic field, measured at half of the box length, in a logarithmic scale. As expected, the magnetic mode is unstable and grows with a linear growth rate of $\eta _{\mathrm{num}}/\omega _{p}\simeq 0.448\pm 0.003$ , in good agreement with the theoretical value of $\eta _{\mathrm{th}}/\omega _{p}=0.450$ , for a value of $kc/\omega _{p}=1$ , obtained by solving (4.3). Indeed the magnetic energy shows the same behaviour in the full kinetic VLEM1D2V model, exhibiting the same growth rate (not shown here).

Figure 5 shows the behaviour of electrons, respectively, in the $x{-}p_{y}$ and $x{-}p_{x}$ phase space in the top and middle frames using the full kinetic VLEM1D2V code. The distribution functions have been averaged over the lacking dimension, i.e. here in the $p_{x}$ and $p_{y}$ direction, respectively.

Figure 5.

Representation of the distribution function in the $x{-}p_{y}$ phase space (top) and in the $x{-}p_{x}$ phase space (middle) using the full kinetic VLEM1D2V code. In the bottom frame, the corresponding $x{-}p_{x}$ phase-space representation of the sum of the streams in the case of CFI.

The corresponding numerical results obtained from the reduced multi-stream code are shown in the bottom frame in figure 5, at the same time $t\omega _{p}=27$ , during the saturation regime. In order to make a direct comparison in the $x{-}p_{x}$ phase-space plane, we have made the sum of the different bunches $\sum _{j=1,2}f_{j}$ . Clearly, the wave–particle dynamics (trapping and acceleration) is correctly described in the reduced model, even with just two ‘streams’: the detailed mechanism of plasma wave breaking is identical in both simulations.

We clearly see that the behaviour of $f$ in the transverse phase space, i.e. in $x{-}p_{y}$ , shows that the different particle streams responsible for CFI keep a temperature quite close to the temperature given with the initial condition, despite the appearance of an inhomogeneous structure in $\widetilde {f}(x,p_{y},t)=\int fdp_{x}$ . It can be also noted that the linear stage of CFI (so as the classical temperature-driven WI) does not depend on the heat flux term (Sarrat, Del Sarto & Ghizzo Reference Sarrat, Del Sarto and Ghizzo2016). These features of the CFI may thus explain why only two streams ultimately prove sufficient to recover, with a high degree of accuracy, the characteristics of CFI both in the linear regime and in the nonlinear regime of the instability.

4.4. The case of the Weibel instability

We now briefly illustrate the physical mechanism for the thermal anisotropy-driven WI (Weibel Reference Weibel1959; Morse & Nielson Reference Morse and Nielson1971), in a 1-D configuration space. Weibel instabilities are encountered in laser–plasma interaction in the relativistic regime (Silva et al. Reference Silva, Fonseca, Tonge, Mori and Dawson2002; Okada & Ogawa Reference Okada and Ogawa2007; Petersen et al. Reference Petersen, Glenzer and Fiuza2021), in astrophysics (Califano, Pegoraro & Bulanov Reference Califano, Pegoraro and Bulanov1997; Fonseca et al. Reference Fonseca, Silva, Tonge, Mori and Dawson2003; Schaefer-Rolffs, Lerche & Schlickeiser Reference Schaefer-Rolffs, Lerche and Schlickeiser2006; Schaefer-Rolffs & Tautz Reference Schaefer-Rolffs and Tautz2008; Lazar et al. Reference Lazar, Schlickeiser, Wielebinski and Poedts2009; Grassi et al. Reference Grassi, Grech, Amiranoff, Pegoraro, Macchi and Riconda2017; Schoeffler et al. Reference Schoeffler, Loureiro, Fonseca and Silva2016; Schlickeiser Reference Schlickeiser2004; Schaefer-Rolffs & Schlickeiser Reference Schaefer-Rolffs and Schlickeiser2005) and collisionless shocks in astrophysics (Medvedev & Loeb Reference Medvedev and Loeb1999; Medvedev Reference Medvedev2006; Park et al. Reference Park2015).

We first carry out a VLEM1D2V simulation of WI, using as initial condition a Maxwellian distribution function with a temperature anisotropy corresponding to $k_{\mathrm{B}}T_{x}=1\,\mathrm{keV}$ and $k_{\mathrm{B}}T_{y}=50\,\mathrm{keV}$ along the longitudinal $p_{x}$ and perpendicular $p_{y}$ components, respectively. We first consider a single unstable mode $k_{0}$ . The numerical space domain, in dimensionless units, is given by $L_{x}=2\pi /k_{0}$ . Here, we choose $k_{0}c/\omega _{p}=1.75$ . We perturb the system by a magnetic field term $\delta B_{z}=B_{0}\sin k_{0}x$ with $eB_{0}/m\omega _{p}=10^{-4}$ as initial amplitude. The phase-space sampling for the full kinetic VLEM1D2V code is $N_{x}N_{p_{x}}N_{p_{y}}$ equal to $256^{3}$ , and we choose a time step of $\triangle t\omega _{p}=0.005$ .

The introduction of a small perturbation of the magnetic field leads to the separation of the currents in space, thereby producing the amplification of the magnetic field. Thus, both currents and magnetic fields increase exponentially. Indeed it is these alternating currents that drive the instability and lead to the complex ‘Y’ shape of the distribution function met in the $x{-}p_{y}$ phase space, already observed in Vlasov simulations in Ghizzo (Reference Ghizzo2013b ).

Figure 6.

Illustration of the phase-space representation of the distribution function in the $x{-}p_{x}$ plane (top) and in the $x{-}p_{y}$ plane (bottom), obtained from the simulation of WI carried out using the VLEM1D2V code. Note the appearance in the dynamic of $f$ of two kinds of ‘anti-symmetric’ coherent structures in the bottom panel. The plots have been realised by averaging the distribution on the lacking variable (on $p_{y}$ and $p_{x}$ for the top and bottom panels, respectively).

In figure 6, the plots of the electron distribution function in the $x{-}p_{x}$ and $x{-}p_{y}$ phase space are shown in the top frame and in the bottom frame, respectively, at time $t\omega _{p}=71$ , i.e. in the saturation phase of WI. Although the distribution exhibits a complex shape in the $x{-}p_{y}$ phase space, we see that the multi-stream code is able to recover the dynamics of $f$ in the longitudinal $x{-}p_{x}$ phase space, even with a small number of ‘numerical streams’ (five here). It is clear from the plot of $f$ in $x{-}p_{y}$ that the dynamics of the distribution results from a separation of the electron population into two kinds of ‘bunches’ of particles linked to the separation of the currents. These particle bunches seem to be linked directly to the property of invariance of the canonical momentum in the perpendicular direction. The results of the nonlinear VLEM1D2V simulations show not only that the nonlinear saturation is governed by strong magnetic trapping as expected, but also that the concept of ‘stream’ is important in WI.

Because we want to compare the case of the WI driven by a temperature anisotropy with the multi-stream code, we need to know precisely how we may reintroduce the concept of ‘temperature’ along the ‘lacking’ $p_{y}$ variable. This task is made easier by considering the equivalence in the sense of the moments of the distribution function or, in other words, by looking at the equivalence in the fluid momentum sense between the reduced multi-stream distribution and a continuous distribution function (see also Ghizzo et al. (Reference Ghizzo, Sarrat and Del Sarto2017) and Del Sarto et al. (Reference Del Sarto, Ghizzo and Sarrat2024) for the non-relativistic case).

Let us now consider a Maxwellian distribution $F_{M}(\boldsymbol{p})$ of thermal velocity $v_{\mathrm{th},y}$ along the $p_{y}$ coordinate. By defining the averaged quantity $h(p_{y})=\int F_{M}(p_{x},p_{y}){\rm d}p_{x}$ , we can write the $h$ function in form $h(p_{y})=\sum _{j=-N_{s}}^{+N_{s}}\alpha _{j}\delta (p_{y}-C_{j})$ , i.e. in the form of $2N_{s}+1$ ‘streams’, where the coefficients $\alpha _{j}$ and $C_{j}$ allow one to rebuild the shape of the distribution in the $p_{y}$ direction and to recover the concept of temperature $T_{y}$ (Inglebert et al. Reference Inglebert, Ghizzo, Reveille, Del Sarto, Bertrand and Califano2011; Ghizzo & Bertrand Reference Ghizzo and Bertrand2013). We may define the $2n$ moment of $h(p_{y})$ as

(4.4) \begin{equation} \int _{-\infty }^{+\infty }p_{y}^{2n}h(p_{y}){\rm d}p_{y}=\sum _{j=-N_{s}}^{+N_{s}}\alpha _{j}C_{j}^{2n}. \end{equation}

By assuming that the $j{\rm th}$ stream verifies the symmetry properties $C_{-j}=-C_{j}$ and $\alpha _{-j}=\alpha _{j}$ for $j=1,\ldots ,N_{s}$ , together with the choice $C_{0}=0$ , (4.4) leads in the case of a Maxwellian distribution to

(4.5) \begin{equation} {{2\sum _{j=1}^{+N_{s}}\alpha _{j}C_{j}^{2n}=\left (2n-1\right )!!\left (k_{B}T_{y}m\right )^{n}}}. \end{equation}

Equation (4.5), together with the normalisation condition $\sum _{j=-N_{s}}^{+N_{s}}\alpha _{j}=1$ , allows us to reconstruct the distribution in the form of $2N_{s}+1$ streams whose distribution along the $p_{y}$ direction can be used to describe any given type of $T_{y}$ temperature distribution. In the general case, (4.5) has the form of a Vandermonde system, which can become ill-conditioned for large values of the number of streams. Following Inglebert et al. (Reference Inglebert, Ghizzo, Reveille, Bertrand and Califano2012b ), a more convenient way consists of a regular sampling of the $p_{y}$ axis, for an equispaced set of $C_{j}$ values. For instance for $N_{s}=2$ (i.e. five streams), we have to solve the following reduced system:

(4.6) \begin{equation} {\alpha _{0}+2\alpha _{1}+2\alpha _{2}=1}, \end{equation}

(4.7) \begin{equation} 2\alpha _{1}C_{1}^{2}+2\alpha _{2}C_{2}^{2}=m^{2}v_{\mathrm{th},y}^{2}, \end{equation}

(4.8) \begin{equation} 2\alpha _{1}C_{1}^{4}+2\alpha _{2}C_{2}^{4}=3m^{4}v_{\mathrm{th},y}^{4}, \end{equation}

where $v_{\mathrm{th},y}^{2}=k_{\mathrm{B}}T_{y}/m$ . By imposing $C_{0}=0$ , $C_{1}=mv_{\mathrm{th},y}$ and $C_{2}=2mv_{\mathrm{th},y}$ , we obtain $\alpha _{0}=1/2$ , $\alpha _{1}=1/6$ and finally $\alpha _{2}=1/12$ .

A comparison between both VLEM1D2V and the multi-stream codes (using five ‘streams’) has already been presented in figure 4 of Inglebert et al. (Reference Inglebert, Ghizzo, Reveille, Bertrand and Califano2012b ). For the multi-stream code, we have reproduced two simulations using the same phase-space sampling of $N_{x}N_{p_{x}}$ of $512\times 1024$ grid points (with the same time step of $\triangle t\omega _{p}=0.0025$ ), by changing the number of streams. Numerical simulations using the multi-stream code have been carried out using identical parameters of $k_{0}$ , $k_{\mathrm{B}}T_{x}=1\,\mathrm{keV}$ and $k_{\mathrm{B}}T_{y}=50\,\mathrm{keV}$ as previously used in the full kinetic treatment.

The temporal evolution of the magnetic field component $eB_{z}/m\omega _{p}(x={(L_{x}}/{2)},t)$ , measured at half the box length, is shown in figure 7 in the case of just three streams (top) and then seven streams (bottom). The curves show no major difference in the linear regime of the instability, and both vary with the expected growth rate of $\eta _{\mathrm{num}}/\omega _{p}\simeq 0.40\pm 0.003$ , in perfect agreement with the theoretical predictions and expected results of Inglebert et al. (Reference Inglebert, Ghizzo, Reveille, Bertrand and Califano2012b ), demonstrating the model convergence with a fairly modest number of 5–7 streams. Remarkably, the linear phase is perfectly described with just three streams. However, differences appear in the saturation regime, notably in the oscillatory behaviour of the magnetic field, characterised by field oscillation at the bounce frequency induced here by magnetic trapping. This frequency is better described in the multi-stream model, which requires the selection of at least five streams.

Figure 7.

For the study of WI, plot of the magnetic field $eB_{z}(x={(L_{x}}/{2)},t)/m\omega _{p}$ in a logarithmic scale versus time, in the case of the multi-stream code: three streams (top) and seven streams (bottom). The linear WI regime is perfectly described in both cases, although some differences persist in the WI saturation regime. The bounce frequency is well described from five streams upwards. The case of seven streams is similar to that of five streams, showing that convergence is obtained for a low number of streams.

Last but not least, another point deserving some comments concerns the importance of the invariance of the canonical momentum in the complex shape of $f$ , for instance observed in figure 6 (bottom panel). First, a ‘reduced’ Hamiltonian analysis is fundamental in connection with the behaviour of the distribution function. Figure 8 highlights the importance of the trapping process for certain values of the canonical momentum. The top panel shows the behaviour of the distribution function, obtained from the simulation based on the full kinetic approach (the VLEM1D2V code). At the top, the distribution function $f(x,p_{x},p_{y}=p_{y,\mathrm{th}})$ , i.e. plotted for a value of the $p_{y}$ variable close to the electron thermal momentum (with $v_{y,\mathrm{th}}\simeq \sqrt {k_{\mathrm{B}}T_{y}/m}\sim 0.312c$ and $p_{y,\mathrm{th}}=m\gamma _{\mathrm{th}}v_{y,\mathrm{th}}$ , $\gamma _{\mathrm{th}}$ being the Lorentz factor corresponding to the thermal velocity), is presented at time $t\omega _{p}=71$ . The middle and bottom frames in figure 8 correspond, respectively, to plots of the five-stream and seven-stream models, for the stream corresponding to the same value of the momentum $p_{y}=C_{j}-eA_{y}$ . This clearly shows that the complex behaviour of the distribution function in the $x{-}p_{y}$ plane is intrinsically linked to the overall dynamics of the distribution, and in part to the magnetic trapping process of $f_{j}.$ We can see that the main characteristics of the distribution function are perfectly recovered for a relatively small number of streams.

However, some differences are observed in the representation of the distribution function in the $x{-}p_{x}$ phase space. The representation of the central trapping region (O-point) is slightly distorted in the case of the multi-stream model (middle and bottom plots in figure 8), compared with the plot obtained with the fully kinetic code (at the top). Regarding the X-point, at the centre of the plots in figure 8, an increase of ‘noise’ is observed near the X-point, which can lead to an oscillatory behaviour of the distribution function and the appearance of negative values in the distribution function (a numerical artefact, quantitatively negligible, produced by the semi-Lagrangian method). However, these negative values are observed in regions where the density of the distribution function remains very low (of the order of $10^{-5}$ in absolute values), and their absolute value decreases when seven streams are used (see the bottom plot).

The last plot shows the formation of ‘arms’ induced by the rapid rotation of the trapping structure around itself. The trapping mechanism is primarily magnetic in nature. The rotation of the central vortex (trapping structure) induces the formation of two rotating arms that become increasingly thin as they rotate, reminiscent of the filamentation process of the distribution function in velocity space. Filamentation leads to the creation of a finer structure in velocity space and then in phase space. While the filamentation process of $f$ is generally well described in Lagrangian codes such as PIC, it can pose a significant challenge in an Eulerian approach. The multi-stream technique allows for a significant reduction in computation time (of the order of $N_{p_{y}}/N_\textit{stream}$ , i.e. the ratio between the sampling used for the description in the $p_{y}$ momentum space and the number of streams, here close to $256/7\sim 36$ ). This reduction in the computation time can be used to increase the sampling of the distribution function in phase space (the sampling is multiplied by four in the multi-stream model), while maintaining a reasonable CPU time. Therefore, the use of a semi-Lagrangian approach, coupled with the Hamiltonian reduction technique, offers a modelling tool capable of finely tracking this filamentation process.

Figure 8.

Representation of the distribution function in $x{-}p_{x}$ phase space in the case of WI. Top: plot obtained from the VLEM1D2V code for a population of particles positioned at $p_{y}\sim -2p_{\mathrm{th},y}$ ; middle: plot of the localised stream at $p_{y}=C_{-2}-eA_{y}=-2p_{\mathrm{th},y}-eA_{y}$ from the multi-stream code with five streams; bottom: plot of the equivalent stream obtained from the multi-stream code with seven streams. Temperatures are $k_{\mathrm{B}}T_{x}=1\,\mathrm{keV}$ and $k_{\mathrm{B}}T_{y}=50\,\mathrm{keV}$ , respectively, in the $p_{x}$ and $p_{y}$ directions.

Thus the reduced model can provide a deeper physical insight, depending on the kind of model simplification that has been performed. By reducing the general solution to a subclass of exact solutions of the Vlasov–Maxwell system based on the invariance of the generalised canonical momentum, the driving idea is to reduce the complexity of the physical problem, not only from the mathematical and numerical, but also from the conceptual point of view, so that the fundamental ingredients and mechanisms at play can be better recognised. Thus, the multi-stream model is well suited to describe the full nonlinear dynamics of a plasma in the presence of Weibel-type instabilities: the conservation of the canonical momentum appears to play a fundamental role in the saturation scenario of the instability due to magnetic trapping. Additional features related to nonlinear charge-separation effects associated with the linear growth of WIs have been highlighted through the multi-stream representation (see Del Sarto et al. (Reference Del Sarto, Ghizzo and Sarrat2024) for more details).

4.5. The case of the oblique current filamentation instability

The example here discussed corresponds to the case of so-called ‘oblique’ CFI modes. Choosing the $\boldsymbol{k}$ wavevector in the $(x,y)$ plane, with two electron (physical) beams always positioned along the $p_{y}$ direction, leads to a kind of coupling between CFI, described in § 4.3, and the two-stream instability (TSI), of a purely electrostatic nature, which requires modelling in two spatial dimensions. In addition, this case provides a test case for evaluating the effectiveness of the various 2-D advections, $\mathcal{X}$ , $\mathcal{Y}$ and $\mathcal{R}$ which are made in the respective $x{-}p_{x}$ , $y{-}p_{y}$ and finally $p_{x}{-}p_{y}$ spaces.

4.5.1. The electrostatic–electromagnetic nature of the initial perturbation

Guided by the linear analysis of the dispersion relation of oblique Weibel-type modes reported in Ghizzo et al. (Reference Ghizzo, Del Sarto and Sarrat2020a ) and Ghizzo & Del Sarto (Reference Ghizzo and Del Sarto2020b ), we start with a simulation carried out with the VLEM2D3V code for the study of the OI in the weak relativistic regime. The system is initialised by two symmetric electron beams of normalised velocities $\beta _{1}=v_{y1}/c\simeq -0.67$ and $\beta _{2}=v_{y2}/c\simeq 0.67$ , counter-propagating in the $y$ direction. Charge and current neutrality are initially ensured by imposing $n_{1}+n_{2}=n_{0}$ and $n_{1}\beta _{1}+n_{2}\beta _{2}=0$ for all simulations corresponding to the non-propagative branch of the OI: the box lengths are $L_{x}=4.088d_{e}$ and $L_{y}=2.513d_{e}$ (where $d_{e}=c/\omega _{p}$ ). A sampling of the distribution function in the configuration space is chosen to be $N_{x}N_{y}=256^{2}$ , i.e. the simulation uses a grid spacing of $\Delta x=0.0160d_{e}$ and $\Delta y=0.00985d_{e}$ . The most unstable mode (determined via the linear dispersion relation in Ghizzo et al. (Reference Ghizzo, Del Sarto and Sarrat2020a)) is excited with $k_{x}d_{e}\simeq 3.05$ and $k_{y}d_{e}=5$ , which correspond to Fourier modes $(n_{x},n_{y})=(2,2)$ . The initial distribution function is composed of two drifted Maxwellian distributions $f(\boldsymbol{x},\boldsymbol{p},t=0)=\sum _{j=1,2}n_{j}F_{0j}(\boldsymbol{p})$ , where

(4.9) \begin{equation} F_{0j}={{\frac {1}{\left (2\pi m^{2}c^{2}\beta _{\mathrm{th}}^{2}\right )^{3/2}}}}\exp \left (-\frac {p_{x}^{2}+\left (p_{y}-p_{j}\right )^{2}+p_{z}^{2}}{2m^{2}c^{2}\beta _{\mathrm{th}}^{2}}\right)\!, \end{equation}

and where $\beta _{\mathrm{th}}=\sqrt {k_{\mathrm{B}}T_{\mathrm{eq}}/(mc^{2})}$ is the thermal velocity normalised to the light velocity $c$ , $p_{j}=m\varGamma _{0j}\beta _{j}c$ and $\varGamma _{0j}=1/\sqrt {1-\beta _{j}^{2}}$ the corresponding Lorentz factor. Here we choose, for each propagating electron beam, the same temperature of $k_{\mathrm{B}}T_{\mathrm{eq}}=6\,\mathrm{keV}$ . The symmetry of the temperature of the beams is a further requirement in order to afford the decoupling of electrostatic and electromagnetic modes (see Sarrat et al. (Reference Sarrat, Del Sarto and Ghizzo2016) for more details). The distribution function is perturbed by purely transverse magnetic fluctuations introduced in the $z$ component of the magnetic field in the form of a sine perturbation $\delta B_{z}=\delta B_{0}\sin (n_{x}\Delta k_{x}x+n_{y}\Delta k_{y}x)$ , allowing one to excite the most unstable mode of wavevector components $(k_{x},k_{y},0)$ with a perturbation amplitude of $e\delta B_{0}/m\omega _{p}=0.005$ . A sampling of $N_{p_{x}}N_{p_{y}}N_{p_{z}}$ of $64^{3}$ points was used in the momentum space and the time step is $\Delta t\omega _{p}=0.008$ .

The multi-stream model is integrated by solving the reduced Vlasov equations given in (3.45), with five streams, coupled in a self-consistent way with (3.46) and with the Maxwell equations. For the initial distribution function, we consider a similar case to that previously studied in the VLEM2D3V code. For modelling a temperature in the $p_{z}$ direction, we again use (4.6). In the case of the multi-stream model, a small perturbation is applied on both transverse magnetic components in the following form:

(4.10) \begin{align} E_{x}&=-\frac {k_{y}\delta B_{0}\omega _{p}}{k_{x}^{2}+k_{y}^{2}}\sin \left (k_{x}x+k_{y}y\right )\!,\quad E_{y}=\frac {k_{x}\delta B_{0}\omega _{p}}{k_{x}^{2}+k_{y}^{2}}\sin \left (k_{x}x+k_{y}y\right )\!,\nonumber\\&\quad B_{z}=\delta B_{0}\sin \left (k_{x}x+k_{y}y\right )\!. \end{align}

The choice of this perturbation is due to the fact that, for oblique propagation with respect to the beams, the CFI, quasi-static if excited alone with a perpendicular propagation, couples with the TSI by becoming propagative: both $E_{x}$ and $E_{y}$ components are thus involved, although only the component $B_{z}$ of the magnetic field grows in time.

The magnetic field generated in the purely electromagnetic CFI is quasi-static, as it is spatially localised and typically does not propagate. When the initial perturbation is strictly perpendicular to the direction of the antiparallel electron beams and the beams are symmetric, the instability remains purely electromagnetic. In this case, the linear growth can be entirely decoupled from electrostatic fluctuations. However, if the wavevector of the perturbation becomes oblique to the orientation of the beams, the CFI couples with the electrostatic TSI and leads to the so-called OI. This coupling introduces a time-dependent characteristic to the electrostatic mode. As a result, the combined CFI–TSI mode imparts a propagative nature to the magnetic field perturbation that gives rise to a resonant character of the instability. Note that, in the considered geometry, the only electromagnetic component of the OI is related to the perturbation on $A_{z}$ . Therefore, the electromagnetic perturbation on the component $E_{z}$ is completely decoupled from the electrostatic fluctuations associated with $E_{x}$ and $E_{y}$ .

Thus, the manner in which the initial perturbation is introduced into the system can lead to a transient regime that varies in duration in the multi-stream code. The length of the transition phase is highly dependent on the accuracy with which the unstable eigenmode is represented. While a simple perturbation of the magnetic field components is sufficient to initiate the instability in the VLEM2D3V code, a coherent perturbation must be introduced in the multi-stream code, affecting both the electromagnetic fields and the vector potential, as well as the distribution function. This difference is related to the construction of the multi-stream code, which somewhat ‘decouples’ the dynamics of the distribution function $f$ obtained from the advections, as described by (3.51)–(3.52), and the vector potential calculated from equation (3.46), from the behaviour of $f$ described by the global rotation in (3.53).

Concerning the transverse electric components, a similar perturbation has been introduced in the form

(4.11) \begin{align} E_{z}&=\delta E_{0}\sin \left (k_{x}x+k_{y}y\right ),\quad B_{x}=\frac {k_{y}\delta E_{0}}{\omega _{0}}\sin \left (k_{x}x+k_{y}y\right )\!,\nonumber\\&\quad B_{y}=-\frac {k_{x}\delta E_{0}}{\omega _{0}}\sin \left (k_{x}x+k_{y}y\right )\!, \end{align}

where $\omega _{0}^{2}=\omega _{p}^{2}+(k_{x}^{2}+k_{y}^{2})c^{2}$ . A corresponding perturbation is added in the $z$ component of the vector potential:

(4.12) \begin{equation} A_{z}=-\frac {\delta E_{0}}{\omega _{0}}\cos \left (k_{x}x+k_{y}y\right )+\frac {\delta E_{0}}{\omega _{0}}\times 10^{-3}\sum _{n_{x},n_{y}=1}^{15}\sin (n_{x}\Delta k_{x}x+n_{y}\Delta k_{y}y)\!, \end{equation}

together with a small-amplitude perturbation to excite a broad spectrum of unstable modes, characterised by (4.12). Concerning the distribution function, both the perturbation in density and momentum are introduced in a self-consistent way with the perturbations given on the electromagnetic field components. Thus we have

(4.13) \begin{equation} f\!\left (\boldsymbol{x},\boldsymbol{p},t=0\right )=\left [1+\delta n\cos \left (k_{x}x+k_{y}y\right )\right ]\sum _{j=1,2}n_{j}F_{0j}\!\left (p_{x}-P_{x},p_{y}-\beta _{j}c\varGamma _{0j}-P_{y}\right)\!, \end{equation}

where the perturbed (fluid) momenta are $P_{x}={(ek_{y}\delta B_{0}}/{(k_{x}^{2}+k_{y}^{2}))}\cos (k_{x}x+k_{y}y)$ , $P_{y}={(ek_{x}\delta B_{0}}/{(k_{x}^{2}+k_{y}^{2}))}\cos (k_{x}x+k_{y}y)$ and $\delta n/n_{0}=0.005$ , $e\delta B_{0}/m\omega _{p}=0.005$ and $e\delta E_{0}/m\omega _{p}c=0.05$ . In (4.12), the quantities $\Delta k_{x}=2\pi /L_{x}$ and $\Delta k_{y}=2\pi /L_{y}$ are the fundamental Fourier modes in the $x$ and $y$ directions, respectively. We choose their normalised values to be $\Delta k_x d_e \simeq 1.53$ and $\Delta k_y d_e = 2.50$ . The perturbation is applied to the mode $(n_x,n_y)=(2,2)$ , which corresponds to the wavevector components $k_x d_e \simeq 3.05$ and $k_y d_e = 5$ . The phase space sampling is $N_{x}N_{y}=256^{2}$ and $N_{p_{x}}N_{p_{y}}=128^{2}$ , the time step used in simulations is $\Delta t\omega _{p}=0.005$ .

4.5.2. Nonlinear saturation processes

The CFI is the main basic plasma process which generates a strong magnetic field $B_{z}$ . In its oblique version (OI), studied here, a coupling with the electrostatic TSI is also possible, allowing also for the growth of an electrostatic component of the electric field. This configuration corresponds to a 2-D system, in which the kinetic energy density $\epsilon _{K}=({1}/{L_{x}L_{y})}\int {\rm d}^{2}x\int {\rm d}^{3}p\:mc^{2}(\gamma -1)f(x,p,t)$ is converted into magnetic energy density $\epsilon _{m}={(1}/{L_{x}L_{y})}\int {\rm d}^{2}x{(1)}/{(2)}\varepsilon _{0}c^{2}|\boldsymbol{B}|^{2}$ and electric energy density $\epsilon _{e}={(1}/{L_{x}L_{y})}\int {\rm d}^{2}x{(1)}/{(2)}\varepsilon _{0}|\boldsymbol{E}|^{2}$ . From a numerical point of view, the various energy densities are normalised to $n_{0}mc^{2}$ .

Figure 9.

Time evolution of the $z$ component of the magnetic energy $\epsilon _{m,z}$ , in a logarithmic scale, in the case of an OI with a hybrid character, both electromagnetic (CFI) and electrostatic (TSI), for a system consisting of two electron counter-streaming beams. The black curve corresponds to the VLEM2D3V code, while the blue curve was obtained from the multi-stream (ms) code (with five streams). Physical parameters are identical in both simulations.

In figure 9, we have superimposed the time evolution of the magnetic energy in a logarithmic scale obtained from the VLEM2D3V code (in black) and that obtained from the multi-stream code (in blue). The curves show the same growth rate, close to $\eta _{\mathrm{num}}/\omega _{p}=0.405\pm 0.003$ , in agreement with the theoretical maximum growth rate $\eta _{th}/\omega _{p}\sim 0.40$ predicted by the dispersion relation obtained by Ghizzo et al. (Reference Ghizzo, Del Sarto and Sarrat2020a ). The curves show a time delay in the start-up of the instability (close to the time difference $t_{\mathrm{delay}}\omega _{p}\sim 60$ ), with the OI in the VLEM2D3V code starting later. This time delay is attributed to a difference in the perturbations introduced in the two different numerical codes.

In both numerical approaches (fully kinetic with the VLEM code and reduced with the multi-stream code), the Vlasov equation is solved using a semi-Lagrangian scheme, which exhibits a very low level of numerical noise. The start-up of the instability is highly sensitive to the type of perturbation initially introduced. If the component of the vector potential $A_{z}$ is not initially perturbed, the code remains perfectly stable over very long time scales.

Thus, the multi-stream code requires the introduction of a coherent perturbation not only in the components of the electromagnetic fields $(\boldsymbol{E,B})$ and the vector potential $\boldsymbol{A}$ , but also in the distribution function, which accounts for perturbations in both density and velocities induced by the initial electromagnetic field. Such a perturbation is necessary to initiate the WI, as the code begins in a noiseless state with the electrostatic and electromagnetic components of the electric field being decoupled from the outset. In contrast, the VLEM2D3V code does not require such modifications to the initial conditions; only a perturbation in the magnetic field component is necessary. However, we subsequently introduced the same modifications in the VLEM2D3V code to facilitate the comparison shown in figures 15–19.

Figure 10 shows the temporal evolution of the magnetic energy $\epsilon _{m,z}$ calculated from the $B_{z}$ component (top) as well as the evolution of the average density (bottom). A perturbation is introduced in the components of the electromagnetic field, similar to the simulation case shown in figure 9.

Figure 10.

An example of the OI in which the coupling between CFI and TSI does not take place. The initial perturbation is introduced only in the fields but not in the distribution function: because of the way the algorithm advances quantities, this induces a much longer transient before the OI eigenmode is formed and starts growing. Therefore, for a long time interval the system remains stable. Top: the time evolution of the $z$ component of the magnetic energy $\epsilon _{m,z}$ . Bottom: the time evolution of the mean density which exhibits a good conservation.

Note that this kind of perturbation in the multi-stream code does not allow a fast onset of the OI, because of the reasons provided in § 4.5.1, related to the way the numerical integration is performed. In this simulation case, the particle density is very well conserved (the relative density variation is here limited to $10^{-8}$ ). The simulation case depicted in figure 10 shows the excellent stability of the code, since the initial metastable equilibrium condition is preserved up to about $10^{5}$ iterations. The variation ( $\sim\!10^{-8}$ ) of the relative density begins at about $t\omega _{p}\simeq 200$ due to the onset of the instability from a long transient time (the given perturbation allows, nevertheless, the onset of a small-amplitude OI). The density conservation could be in principle improved by increasing the refinement of the phase-space mesh. However, its imperfect conservation cannot be avoided, since it is mostly due to the open boundary conditions in the momentum space: this, combined with the broadening of the distribution function (visible, for the case of an OI, on the $p_{x}$ component in the bottom frame of figure 12), unavoidably makes the particles on the tail of the distribution ‘escape’ from the simulation box.

When, instead, the initial perturbation is introduced in the multi-stream model on both the electromagnetic field components and the initial distribution function (to ensure in-phase start-up of the electromagnetic field and vector potential), the OI eigenmode develops over a much shorter transient time. However, the instability growth rate remains unchanged, and the magnetic field fluctuation level and magnetic bounce frequency remain similar in both numerical approaches, even in the instability saturation regime. It should also be noted that the polarisation of the electromagnetic wave remains quasi-linear, with a dominant TM mode component.

Figure 11.

Top: time evolution of the mean density (total mass) obtained in the case of the multi-stream (ms) code (with five streams) using the numerical scheme shown in figure 2, based on the alternating use of 2-D advections and a time-splitting scheme. Bbottom: the corresponding total energy versus time. Note that the average relative density is kept at $(1.005-1.000)/1=0.5\,\%$ for 32 000 iterations, i.e. a simulation time of $t\omega _{p}\simeq 160$ .

Figure 12.

Representation of the distribution function in the $x{-}p_{x}$ phase space. Top: plot of the mean distribution $\widetilde {f}(x,p_{x})=\int f(x,y=L_{y}/2,\boldsymbol{p},t){\rm d}p_{y}\,{\rm d}p_{z}$ , obtained from the full kinetic VLEM2D3V code. The colour plots in the middle and bottom frames correspond to the distribution of the central beam (positioned in $p_z \sim C_0=0$ ) and the last stream obtained from the multi-stream code (with five streams). Note that the last stream (bottom), corresponding to a momentum of $p_{z}\sim 2p_{z,\mathrm{th}}$ (twice the thermal momentum), is made up of an electron population of lower density. Thus, the multi-stream (ms) code is perfectly suited to a fine description of wave–particle interaction, including in the tail regions of the distribution function.

Figure 11 shows the conservation of mass and energy in a simulation involving 2-D advections. The relative total mass is conserved within $(1.005-1.000)/1=0.5\,\%$ when applying advection methods $\mathcal{X}$ $(x,p_{x})$ , $\mathcal{Y}$ $(y,p_{y})$ and $\mathcal{R}$ $(p_{x},p_{y})$ . In contrast, total energy conservation is afforded within a relative error of $(0.3580-0.3568)/0.3568\simeq 0.3\,\%$ over 32 000 iterations, i.e. a simulation time of $t\omega _{p}\simeq 160$ (for a time step of $\Delta t\omega _{p}=0.005$ ). However, mass conservation is less effective compared with the results shown in figure 10. This can be attributed to the instability causing significant heating of the distribution function, which broadens the distribution in momentum variables $p_{x}$ and $p_{y}$ as observed in figure 12. A much better conservation of density would require an increase in the size of the simulation box in momentum space, thereby enhancing the sampling of the distribution function with respect to the momentum coordinates. While this is entirely feasible, it necessitates an increase in CPU computation time. An alternative can be to employ an adaptive mesh refinement technique by incorporating grid refinement strategies as in Antoine et al. (Reference Antoine, Ghizzo, Deriaz and Del Sarto2025), so as to reduce computational time. In figure 12 plots of the distribution function $\widetilde {f}(x,p_{x})=\int f(x,y=L_{y}/2,\boldsymbol{p},t){\rm d}p_{y}\,{\rm d}p_{z}$ are presented showing when the saturation regime is well established.

At the top of figure 12, the plot of the distribution function corresponds to the result obtained with the fully kinetic VLEM2D3V code, which accounts for the complete dynamics of $f$ in momentum space (i.e. there are no ‘streams’ here). This colour plot corresponds to the time $t\omega _{p}=400$ , when the saturation of the OI is well established. The middle panel highlights a similar dynamics obtained with the multi-stream code at an equivalent time (i.e. when the instability has reached the saturation regime). This plot highlights the dynamics of $f$ for the central stream, specifically for the population of electrons that belong to the core of the distribution (with a momentum $p_z\sim 0$ ). In the saturation regime, the dynamics of $f$ is characterised by the emergence of two magnetic trapping structures that rotate one around the other. As a result, the distribution function experiences significant broadening in $p_{x}$ , indicative of a heating process.

The multi-stream model provides a detailed description of wave–particle interactions, even in regions of phase space with very low densities. An example is illustrated in the bottom plot of figure 12, which shows the behaviour of the distribution function for a momentum value of $p_{z}\sim 2p_{\mathrm{th},z}$ (twice the thermal momentum), corresponding to the distribution function of the fifth ‘stream’ (with an initial density much lower compared with the central stream). Such an analysis would require a very large number of ‘particles’ in a PIC code and/or very high sampling of $f$ in momentum space in the case of a direct approach using the VLEM2D3V code. It is observed that the trapping process, which is of a mixed nature – both electrostatic (induced by TSI) and magnetic (induced by CFI) – is more significant in regions of very low densities.

The bottom frame in figure 12, shows the corresponding behaviour of the distribution function $\widetilde {f}_{j=2}(x,p_{x})=\int f_{j=2}(x,y=L_{y}/2,\boldsymbol{p},t){\rm d}p_{y}$ , averaged over the $p_{y}$ component, obtained from the multi-stream code, with five streams. It can be noted that the actual density of the distribution function, at a location of $p_{z}\sim 2p_{z,\mathrm{th}}$ in momentum space, must be corrected by the normalisation factor $\alpha _{2}={1}/{12}$ (recalling that $\alpha _{0}={1}/{2}$ and $\alpha _{-1}=\alpha _{1}={1}/{6}$ and $\alpha _{-2}=\alpha _{2}$ ).

In the case where the number of streams in the code is increased, for example to $2N_{s}+1=11$ , this normalisation factor $\alpha _{5}$ becomes of the order of $10^{-5}$ for the eleventh stream, which corresponds approximately to a position of the last stream at $p_{z}=5p_{z,\mathrm{th}}$ , meaning it is in the tail of the Gaussian distribution. This demonstrates that it is possible to describe, at a lower CPU cost, the wave–particle interactions occurring in regions of very low density, provided that one of the streams is positioned in the tail of the Gaussian distribution.

The plots exhibit a structure characterised by the appearance of two vortices, due to strong magnetic trapping, and a heating process in the longitudinal direction (according to momentum along $p_{x}$ ). Note, in the lower plot in figure 12, that the heating process is slightly more pronounced for a population of electrons chosen at the value of $p_{z}\sim 2p_{z,\mathrm{th}}$ . These results demonstrate that the multi-stream model can be used to recover the electron dynamics during OI at a lower numerical cost.

The reconstruction of the global distribution function from the five streams does not allow for a precise determination of the nature of wave–particle interactions. For example, the mechanism of saturation of CFIs is often induced by the formation of magnetic trapping structures. Such structures can appear in the tails of the distribution function at very low densities: they represent information that is often ‘hidden’ in the global distribution function and is difficult to observe directly. Thus, the early PIC simulations by Morse & Nielson (Reference Morse and Nielson1971) highlighted magnetic trapping structures through the implementation of a type of diagnostics that utilises the invariance of the canonical momentum, which allows for the selection of the ’population’ of super-particles involved in trapping. An example of how magnetic and electrostatic trapping differently affect streams of particles corresponding to different values of the momentum was shown in Del Sarto et al. (Reference Del Sarto, Ghizzo and Sarrat2024) (see e.g. figure 15 therein). In general, the multi-stream technique allows for an easier characterisation of the role played by different types of particle trapping, thanks to having partitioned the distribution into ‘streams’. The actual reconstruction of the streams is used in the code only for determining source terms.

4.6. Multi-stream simulation using five streams: coupling CFI–magnetic reconnection

The last example, which we discuss in this subsection, demonstrates the capabilities of the multi-stream code to model a case of coupling between OI and WI, in which a temperature anisotropy is combined with two electron counter-propagating beams. This last example shows the possibility of the emergence of a MR process, previously observed by Califano et al. (Reference Califano, Attico, Pegoraro, Bertin and Bulanov2001), driven by the plasma self-organisation.

The topological configuration, used here, where the beams are propagating inside the $(x,y)$ plane and where a MR process takes place, allows us to study the coupling of OI with MR in presence of a secondary instability, as for instance WI, driven by a temperature anisotropy with $T_{\perp }=T_{z}\geqslant T_{\parallel }=T_{x}=T_{y}$ . The relative importance of energy conversion and topological change as defining properties of reconnection and the extent to which they are linked are key issues in MR.

Two simulations have been performed using the full kinetic VLEM2D3V and the multi-stream code (with five streams) with the same physical parameters as before: the in-plane temperature is taken identical in all cases and maintained at $k_{\mathrm{B}}T_{\parallel }=6\,\mathrm{keV}$ , with a perturbation introduced into both the components of the electromagnetic field and the initial distribution function (perturbed in both density and velocity). We take a lower level of perturbation in amplitude, i.e. $\delta n/n_{0}=0.005$ , $e\delta B_{0}/m\omega _{p}=0.005$ and $e\delta E_{0}/m\omega _{p}c=0.005$ .

The simulation, carried out with VLEM2D3V, uses numerical parameters analogous to the results presented in § 4.5, except for the sampling of $N_{p_{x}}N_{p_{y}}N_{p_{z}}=32^{3}$ points for representing the distribution function in momentum space, and for values of the momentum vector included in the intervals $|p_{x}|,|p_{y}|\leqslant 3mc$ and $|p_{z}|\leqslant 5mc$ . A second simulation, for identical physical parameters, has been performed with the multi-stream code, with five streams, with a higher sampling of $N_{p_{x}}N_{p_{y}}=128^{2}$ points and with a time step of $\triangle t\omega _{p}=0.008$ .

Identical perturbations have been implemented in both codes, as well as a magnetic perturbation term. Here, both OI and WI are initially excited at the same time, leading to the generation and amplification of an electromagnetic field with two types of modes: a TM mode, induced by OI, which leads to the amplification of an electromagnetic wave governed by the $(E_{x},E_{y},B_{z})$ components, in agreement with the results of § 4.5, and a TE mode, induced by the coupling with WI, of $(E_{z},B_{x},B_{y})$ components implicated in the MR process.

Figure 13 shows the dynamics of the magnetic field lines in the plane $(x,y)$ , for the $B_{x}{-}B_{y}$ components of the magnetic field, as the result of the coupling between OI and WI. At the top, the reconnection processes in the $(x,y)$ plane are clearly visible at time $t\omega _{p}=48$ in the plot of the magnetic field lines obtained from the simulation performed by the VLEM2D3V code. The bottom panel shows the result obtained from the reduced model at an equivalent time: we observe a very good agreement between the two models, showing the emergence of an ultrafast process of reconnection of the field lines in the $(x,y)$ plane, induced here by the coupling between the OI and WI.

Figure 13.

Illustration of the lines of the components $(B_{x},B_{y})$ of the magnetic field in the $(x,y)$ plane. The curves are plotted during a first MR process in the nonlinear saturation regime in presence of a temperature anisotropy; the streams have a temperature $T_{z}$ larger than the parallel temperature. The simulation has been performed using the VLEM2D3V code (top) and the multi-stream model with five streams (bottom).

Figure 14 shows the dynamics of the magnetic field lines at two later instants, obtained from the multi-stream code. The density of the field lines is the same as that used in figure 13. The top panel shows the formation of two chains of magnetic islands at time $t\omega _{p}=155$ , in the saturation regime, located at $x\omega _{p}/c\sim 1$ and $x\omega _{p}/c\sim 3$ , each exhibiting around 15 reconnection sites. The dynamics of these chains, however, remains complex and displays an oscillatory behaviour over time. The chains of islands show a certain stability over time, although there are modulations in the size of the magnetic islands. The magnetic components $(B_{x},B_{y})$ are amplified through the WI, introduced by an anisotropy in the plasma temperature in the direction perpendicular to the reconnection plane (in $T_{z}$ here). Note that both chains of magnetic islands are almost homogeneous in the $y$ direction.

Figure 14.

Nonlinear evolution of the system, showing magnetic reconnection/annihilation processes, obtained with the multi-stream model. The system is identical to that observed in figure 13 in the nonlinear saturation regime at later times, in the presence of temperature anisotropy. About 15 reconnection sites are observed, forming two chains of magnetic islands (top). The self-organisation of the system leads to processes of merging the islands until the appearance of two mesoscopic islands; the streams have a temperature $T_{z}$ larger than the parallel temperature. The simulation has been performed using the multi-stream model with five streams.

The bottom panel in figure 14 shows a process of self-organisation of the magnetic field that leads to the emergence of two large-scale magnetic vortices, occupying the size of the box. This results from the coupling between the CFI (current-driven instability) associated with the two electron counter-propagating beams and the WI induced by the temperature anisotropy.

Figure 15.

Corresponding behaviours of the $z$ component of the magnetic energy versus time for the different numerical approaches: the full kinetic VLEM2D3V code (in black) together with the corresponding multi-stream code (in blue). Top: evolution versus time. Bottom: the same diagnostic is now presented in a logarithmic scale. We have used five streams in the multi-stream (ms) code to describe the OI–WI interaction.

A global view of the temporal dynamics of the magnetic field is presented in figure 15, which shows the $z$ component of the magnetic energy versus time (top) for the two simulations performed using the VLEM2D3V code (in black) and the multi-stream code (in blue). A slight time delay between the energy amplification in the linear regime is clearly observed. The growth rate remains identical in both simulations, as expected. A plot in a logarithmic scale is added in the bottom frame.

The maximum amplitude of the magnetic energy is, however, slightly higher in the case of the multi-stream model. In both simulations, the dominant component of the magnetic field $eB_{z}/m\omega _{p}$ is driven by CFI.

The oscillatory behaviour in the multi-stream simulation, observed for instance in figure 15 (top panel), and carried out with five streams and a temperature of $T_{z}=300\,\mathrm{keV}$ , has a numerical ‘bounce’ frequency of the order of $\omega _{\mathrm{num}}\simeq 0.46\omega _{p}$ in the nonlinear saturation regime, i.e. close to the magnetic bouncing frequency, estimated from the relation (see Davidson et al. Reference Davidson, Hammer, Haber and Wagner1972; Ghizzo Reference Ghizzo2013b )

(4.17) \begin{equation} \frac {\omega _{b}}{\omega _{p}}=\sqrt {\frac {k_{0}c}{\omega _{p}}\frac {p_{1}}{m\varGamma _{01}}\frac {\omega _{ce}}{\omega _{p}}}\simeq \sqrt {1.536\times \frac {0.9}{1.3453^{2}}\times 0.35}\simeq 0.51. \end{equation}

In relation (4.17), we have used a value of $0.35$ for the maximum of the magnetic field component $eB_{z}/m\omega _{p}$ (to calculate the ratio of the electron cyclotron frequency to the plasma frequency $\omega _{ce}/\omega _{p}=eB_{z,\mathrm{max}}/m\varGamma _{01}\omega _{p}$ ) where the Lorentz factor $\varGamma _{01}$ is introduced to take into account the relativistic correction; the wavevector $k_{0}c/\omega _{p}=1.536$ corresponding to the CFI on the mode $(n_{x},n_{y})=(1,0)$ .

We obtain a lower estimate than the theoretical value given by formula (4.17), indicating that magnetic trapping is not the only process involved in the nonlinear saturation of the instability. Thus, electrostatic trapping (induced by the ‘TSI-like component’ of the OI) plays an important role in the saturation of the OI. Another process may also be involved, related to the generation of temperature anisotropy: since $T_{z}\gg T_{\parallel }$ , specifically a secondary WI can develop, which leads to the amplification of TM modes and MR process. The spontaneous generation of a secondary WI, developing in the nonlinear stage of primary relativistic CFI-type modes, has been already observed by Tomita and Ohira (Reference Tomita and Ohira2016).

In the presence of a WI, the dynamics of the distribution function in the phase space can become particularly complex (as observed for example in figure 6 in § 4.4). In order to study this aspect, it is interesting to analyse in more detail the behaviour of the electron distribution, in the VLEM2D3V simulation, during the start of the plasma reorganisation process that leads to the reconnection mechanism.

Figure 16.

Representation of the distribution function $\widetilde {f}(p_{x},p_{z})=\int {\rm d}p_{y}f(x={(L_{x}}/{2)},{}y={(L_{y}}/{2)},p_{x},p_{y},p_{z})$ in the $p_{x}{-}p_{z}$ momentum space, obtained from the VLEM2D3V code, in the case of the coupling of WI and OI, for a relativistic transverse temperature of $k_{\mathrm{B}}T_{\perp }=k_{\mathrm{B}}T_{z}=300\,\mathrm{keV}$ .

We present in figure 16 the behaviour of the distribution in the momentum space $p_{x}{-}p_{z}$ , at two different times. In the top panel of figure 16, the initial distribution is represented, showing the distribution’s anisotropy in temperature along the ‘hot’ direction $p_{z}$ . In the bottom panel of figure 16 the distribution is plotted at time $t\omega _{p}=48$ , i.e. at the beginning of the saturation process, when the magnetic energy begins to decrease in figure 15. A heating process is observed in the longitudinal direction along $p_{x}$ , associated with the decrease in the magnetic field. As we see in figure 17, this mechanism is indeed associated with an energy transfer process towards the $p_{x}$ component, although the multi-stream model slightly overestimates it. This energy transfer occurs via the $xx$ component of the second-order anisotropic fluid moment:

(4.18) \begin{equation} \varPi _{ij}=\int \frac {p_{i}p_{j}}{m\gamma (x,y,p_{x},p_{y})}f(x,y,p_{x},p_{y}){\rm d}p_{x}\,{\rm d}p_{y}. \end{equation}

This tensor encompasses the contribution of both the kinetic energy and of the internal energy of the system.

Note that the tensor $\varPi _{ij}$ is not an extension of the classical (non-relativistic) pressure tensor $\boldsymbol{\pi }=nm(\langle \boldsymbol{v}\boldsymbol{v}\rangle -\boldsymbol{uu})$ in the relativistic regime. Tensor $\boldsymbol{\pi }$ was previously introduced with (4.14)–(4.16). Nevertheless, the tensor $\varPi _{ij}$ can provide an indication of the presence of anisotropy (off-diagonal term), allowing for the identification of potential couplings between various instabilities.

In MR, wave–particle interactions give rise to an energy transfer from the magnetic energy to kinetic energy. This is supported by the plot of the spatial average of the tensor component $\varPi _{ik}$ , which in the multi-stream model is given by

(4.19) \begin{equation} \langle \varPi _{ik}\rangle =\frac {1}{L_{x}L_{y}}\sum _{j=-N_{s}}^{N_{s}}\int {\rm d}x{\rm d}y\int {\rm d}p_{x}{\rm d}p_{y}\,\frac {p_{i}p_{k}}{m\gamma _{j}(x,y,p_{x},p_{y})}f_{j}, \end{equation}

versus time, which is shown in figures 17–19. It must be pointed out that a similar definition can be used to define the quantity $\varPi _{zz}$ just by replacing $p_{z}$ with $p_{z}\equiv C_{j}-eA_{z}(x,y,t)$ . The corresponding formula in the full kinetic treatment is

(4.20) \begin{equation} \langle \varPi _{ij}\rangle =\frac {1}{L_{x}L_{y}}\int {\rm d}x{\rm d}y\int {\rm d}^{3}p\,\frac {p_{i}p_{j}}{m\gamma (\boldsymbol{p})}f(x,y,\boldsymbol{p},t). \end{equation}

Figure 17.

Plot of the $\langle \varPi _{xx}\rangle$ component versus time in both numerical approaches, the multi-stream model with five stream (in blue), together with the full kinetic VLEM2D3V code (in black). We have used an initial distribution with a perpendicular temperature $T_{z}=300\,\mathrm{keV}$ leading to MR via the secondary WI.

From figure 18, which presents the dynamics of the $\langle \varPi _{yy}\rangle$ component of the tensor, we observe a significant decrease in the $\langle \varPi _{yy}\rangle$ component associated with the growth of $\langle \varPi _{xx}\rangle$ . It should be noted that the transverse components $\langle \varPi _{xy}\rangle$ , $\langle \varPi _{xz}\rangle$ and $\langle \varPi _{yz}\rangle$ remain negligible compared with the components of the diagonal of the tensor, namely the components $\langle \varPi _{xx}\rangle$ , $\langle \varPi _{yy}\rangle$ and $\langle \varPi _{zz}\rangle$ . However, in figure 18 we note a small difference in the values of the $\langle \varPi _{yy}\rangle$ component, at the initial time $t\omega _{p}=0$ , obtained with the fully kinetic (VLEM2D3V) and with the reduced (multi-stream) versions, the reduced model giving a slight overestimation of the $\langle \varPi _{yy}\rangle$ component.

Figure 18.

Plot of the $\langle \varPi _{yy}\rangle$ component versus time in both numerical approaches, the multi-stream model with five streams (in blue) together with the full kinetic VLEM2D3V code (in black). We have used an initial distribution with a perpendicular temperature $T_{z}=300\,\mathrm{keV}$ leading to MR via the secondary WI.

A global reorganisation takes place in which collisionless wave–particle interactions transfer the energy stored in the beams and in the magnetic field into the kinetic energy of the plasma. While the beginning of the dynamics, driven by this propagative OI mode, is very similar to the dynamics met in § 4.5, a transition towards 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. Such a cascade-like structure can be observed in figure 13 in the structure of the chain of magnetic vortices in the magnetic field lines. This mechanism can induce deep modifications in the energy transfer and can lead to a strong stochastic heating of the distribution function in the $p_{x}$ direction.

Finally, a second ‘heating’ process is also observed in the $p_{z}$ direction as well (in the time interval $25\leqslant t\omega _{p}\leqslant 50)$ , as indicated by the observation of the temporal variation of the $\langle \varPi _{zz}\rangle$ component of the tensor in figure 19. The blue curve represents the evolution of the quantity $\langle \varPi _{zz}\rangle$ as a function of time, estimated from formula (4.19) in the multi-stream code. We have also plotted the same quantity calculated from the distribution function data obtained from the VLEM2D3V code.

The full kinetic VLEM2D3V code provides a detailed view of the dynamics of the $\langle \varPi _{zz}\rangle$ component over time. After the first phase, linked to the transfer of energy $\epsilon _{m}\rightarrow \epsilon _{K}$ (magnetic into kinetic energy), we observe a kinetic cooling process in the dynamics of $\langle \varPi _{zz}\rangle$ , which is much more pronounced in VLEM2D3V, suggesting that a larger number of streams is now required to describe this process in the multi-stream code.

Figure 19.

Plot of the $\langle \varPi _{zz}\rangle$ component versus time in both numerical approaches, the multi-stream model with five streams (in blue) together with the full kinetic VLEM2D3V code (in black). We have used an initial distribution with a perpendicular temperature $k_{\mathrm{B}}T_{z}=300\,\mathrm{keV}$ leading to MR via the secondary WI.

5.1. Full kinetic VLEM code

The Vlasov–Maxwell VLEM code is currently running on the Jean Zay computer of the IDRIS centre in a hybrid OpenMP/MPI version. Performance tests have been performed on this architecture mainly for the 2D3V version of the semi-Lagrangian VLEM code. To study the performance of the parallel algorithm, we introduced the speed-up $S_{\mathrm{MPI}}(\mathrm{\mathit{p}})$ , defined as follows:

(5.1) \begin{equation} S_{\mathrm{MPI}}\left (\mathrm{\mathit{p}}\right )=\mathrm{\frac {\mathrm{Elapsed}\:process\:time\:on\:16\:cores}{Elapsed\:process\:time\:on\:\mathit{p}\:cores}}.\quad \end{equation}

The algorithm is parallelised using a domain decomposition method in the position $(x,y)$ space. The MPI version allows the communication of data necessary for domain decomposition (in particular the transfer of boundary functions of subdomains and first derivatives for Hermite polynomials). These simulations were performed with a phase-space sampling of $N_{x}N_{y}N_{p_{x}}N_{p_{y}}N_{p_{z}}=1024\times 1024\times 16^{3}$ for 500 time iterations. Given the memory size constraints, we used $N_{\mathrm{thread}}=4$ OpenMP tasks for these simulations (it was not possible to use a single OpenMP task here due to the large sampling chosen).

In figure 20, we observe a good behaviour of the speed-up $S_{\mathrm{MPI}}(\mathrm{\mathit{p}})$ as a function of the number of MPI processes p, with an efficiency close to 84 % even with a very high number of MPI processes. This result is mainly due to the type of parallelism chosen for OpenMP, here carried out at fine grains, that is, by using PARALLEL type directives, which corresponds to a parallelisation on the loops; in this case the latency times are more important.

Figure 20.

The theoretical speed-up is plotted as a solid line. The blue circles represent the value obtained from the speed-up during the test simulations. The speed-up is evaluated from formula (5.1), i.e. the effective time spent in the processor (elapsed time). Here, we observe an efficiency of around 94 % up to 2048 MPI processes (for four OpenMP tasks). This study was carried out on the Jean Zay architecture for the VLEM2D3V code. At 4096 MPI processes, the speed-up remains high with an efficiency close to 84 %.

We have now run VLEM2D3V with several OpenMP tasks for a constant number of MPI processes (here chosen to be 1024 MPI processes), using a phase-space sampling of $N_{x}N_{y}=512^{2}$ and $N_{p_{x}}N_{p_{y}}N_{p_{z}}=16^{3}$ , for 500 temporal iterations. The speed-up $S_{\mathrm{OpenMP}}(\mathrm{\mathit{p}})$ is calculated from the following relation:

(5.2) \begin{equation} S_{\mathrm{OpenMP}}\left (\mathrm{\mathit{p}}\right )=\mathrm{\frac {\mathrm{Elapsed}\:process\:time\:on\:one\:OpenMP\,task}{Elapsed\:process\:time\:on\:\mathit{p}\:OpenMP\,tasks}}.\quad \end{equation}

The results obtained are presented in figure 21. The speed-up $S_{\mathrm{OpenMP}}(\mathrm{\mathit{p}})$ increases with the number of OpenMP tasks, $\mathit{\mathrm{\mathit{p}}}$ . The test was carried out up to a maximum number of $\mathrm{\mathit{p}}=16$ threads. Up to 8 threads, we observe a correct efficiency. We show in figure 21 that the efficiency strongly decreases although the speed-up remains more or less constant beyond four OpenMP tasks. However, in the case of this configuration with a high number of MPI processes, the efficiency remains of the order of 50 % at four OpenMP tasks (with a parameter $T_{\mathrm{CPU}}/T_{\mathrm{Elapsed}}\sim 3.574$ close to $N_{\mathrm{thread}}=4$ ) up to 40 % for $N_{\mathrm{stream}}=8$ (with a parameter $T_{\mathrm{CPU}}/T_{\mathrm{Elapsed}}\sim 6.66$ ). On the other hand, the efficiency falls to 20 % for $N_{\mathrm{thread}}=16$ .

Figure 21.

The theoretical speed-up is plotted as a solid line. The blue squares represent the value obtained from the speed-up during the test simulations using the VLEM2D3V code. The speed-up is evaluated using formula (5.2). This study was carried out on the Jean Zay calculator.

The performance gain of up to eight OpenMP tasks is encouraging. It indicates that improvements can be achieved through multi-task coupling and single-processor vectorisation, particularly by rewriting certain sections of the code related to the processing of internal loops. In principle, it is indeed possible to use a coarse-grained parallelisation, i.e. by introducing a second domain decomposition according to the OpenMP tasks. Maybe this improvement will be implemented in the future, but for the moment it is not a priority task, since it would imply a substantial rewriting of the source code. Note also that speed-up and efficiency are measured in effective run conditions using significant resources in memory size.

5.2. The multi-stream scheme

In this subsection, we present a performance analysis of the multi-stream code. The code uses a domain decomposition similar to that of VLEM, enabling the parallel solution of Maxwell’s equations using a numerical scheme analogous to that employed by VLEM2D3V. However, the advections are now mainly 2-D. This implies transpositions of the matrix, in order to use inversion algorithms for tridiagonal matrices, the interpolation technique being realised with a tensor product of 1-D advections. We have thus run the multi-stream code with several MPI processes, for a constant number of OpenMP tasks, here chosen to be four OpenMP threads, using a phase-space sampling of $N_{x}N_{y}=256^{2}, N_{p_{x}}N_{p_{y}}=32^{2}$ , for 500 time iterations. In all the simulations performed, the multi-stream code has five streams and therefore at least five MPI processes are needed. The speed-up is calculated using

(5.3) \begin{equation} S_{\text{MPI-ms}}\left (p\right )=\mathrm{\frac {\mathrm{Elapsed}\:process\:time\:on\:5\:MPI\:processes}{Elapsed\:process\:time\:on\:\mathit{p}\:MPI\:processes}}.\quad \end{equation}

Figure 22.

The theoretical speed-up is shown as a solid line. The blue circles represent the speed-up values obtained from test simulations using the multi-stream code. The speed-up $S_{\text{ MPI-ms}}$ is calculated using formula (5.3). This study was carried out on the Jean Zay computer.

The results are depicted in figure 22, showing a good level of parallelisation up to 1280 processors, despite the use of successive matrix transpositions during 2-D advection changes. It is interesting here to estimate the CPU times required for both simulations, based on a 2-D spatial configuration. For the OI, studied in § 4.5, the simulation carried out using the VLEM2D3V code required 12 successive restarts (for a time of $t\omega _{p}=500$ , i.e. 60 000 time iterations), for 256 MPI processes and 8 OpenMP tasks. This corresponds to a total elapsed time of $7.5\times 10^{8}\,{\rm s}$ for a simulation using a sampling of $N_{x}N_{y}N_{p}^{3}=256^{2}\times 64^{3}$ grid points in the global phase space, i.e. a computation time of $0.356\times 10^{-9}\,{\rm s}$ per time step, per grid point and per core.

The multi-stream code requires a time of $2.9\times 10^{8}\,{\rm s}$ for a sampling of $N_{x}N_{y}N_{p}^{2}=256^{2}\times 128^{2}$ grid points and $256\times 5$ MPI processes and 4 OpenMP tasks and $10^{5}$ iterations, i.e. $0.542\times 10^{-9}\,{\rm s}$ per time step, per grid point and per core. This is of the same order of magnitude as the elementary time, which seems justified since the multi-stream code uses a global matrix inversion technique for global B-spline interpolation. The VLEM code, on the other hand, uses a local interpolation method, which speeds up the computation time. However, given the large memory requirements ( $16\times 64^{3}$ points per MPI process), this leads to an OpenMP efficiency of $4.1$ for eight OpenMP tasks (i.e. around 50 %), whereas the multi-stream version gives an acceleration rate of 3.5 for four OpenMP tasks, i.e. an OpenMP efficiency of 87.5 %. The use of B-splines, on the other hand, delivers a gain of at least a factor of 2 over the more time-sensitive local advection technique. An equivalent simulation with $N_{x}N_{y}N_{p}^{2}=256^{2}\times 64^{2}$ grid points in phase space, and 30 000 temporal iterations, using 1024 MPI processes and four OpenMP tasks, for example, finally leads to a reduction factor of around 33 to achieve a simulation equivalent to that carried out with VLEM2D3V for equivalent accuracy.