2.1. The Vlasov–Poisson system

A fully ionised, collisionless, electrostatic electron–ion plasma, where magnetic field contributions can be neglected, can be described kinetically by the VP system. This system consists of two Vlasov equations, one for ions (i) and one for electrons (e), coupled with the Poisson equation for the electrostatic potential:

(2.1) \begin{eqnarray} \partial _t f_s + \boldsymbol{v} \boldsymbol{\cdot }\boldsymbol{\nabla } f_s + \frac {q_s}{m_s} {\boldsymbol{E}} \boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} f_s = 0, \quad s = i, e, \end{eqnarray}

(2.2) \begin{eqnarray} \epsilon _0 \Delta \phi = -\sum _{s} q_s \int {d}^3v\, f_s, \end{eqnarray}

(2.3) \begin{eqnarray} {\boldsymbol{E}}(\boldsymbol{x},t)=-\boldsymbol{\nabla }\phi (\boldsymbol{x},t), \end{eqnarray}

where $\varDelta$ is the Laplacian operator, $f_s(\boldsymbol{x},\boldsymbol{v},t)$ are the distribution functions of the ions ( $s=i$ ) and electrons ( $s=e$ ), ${\boldsymbol{E}}(\boldsymbol{x},t)$ is the electric field and $\phi (\boldsymbol{x},t)$ is the electrostatic potential; $q_i=e$ and $q_e=-e$ , with $e$ being the fundamental electric charge; $m_s$ stands for the masses of the ions and electrons; $\epsilon _0$ is the vacuum permittivity; and $\boldsymbol{x}$ and $\boldsymbol{v}$ are the spatial and velocity coordinates, respectively.

The electrostatic potential is a solution to the Poisson equation (2.2):

(2.4) \begin{eqnarray} \phi (\boldsymbol{x},t) = \epsilon _0^{-1}\sum _{s}q_s\int {d}^3x' {d}^3v'\, G(\boldsymbol{x},\boldsymbol{x}') f_s(\boldsymbol{x}',\boldsymbol{v}',t), \end{eqnarray}

where $G(\boldsymbol{x},\boldsymbol{x}')$ is the Green function for the Laplacian $\varDelta$ .

The Vlasov equation (2.1) can be formulated as non-canonical Hamilton equations (Morrison Reference Morrison1980, Reference Morrison1982, Reference Morrison1987; Li Reference Li2025), with the following Hamiltonian functional:

(2.5) \begin{eqnarray} {\mathcal{H}}_{_{VP}}=\sum _{s}\int {d}^3x\,{d}^3v \, f_s\left ( m_s\frac {v^2}{2}+ \frac {q_s}{2}\phi \right ) \!, \end{eqnarray}

where $\phi$ is given by (2.4), and the standard particle Poisson bracket (Morrison Reference Morrison1980, Reference Morrison1982; Marsden & Weinstein Reference Marsden and Weinstein1982):

(2.6) \begin{eqnarray} \{F,G\} = \sum _{s}\int {d}^3x\, {d}^3v\, \frac {f_s}{m_s}\left [\frac {\delta F}{\delta f_s},\frac {\delta G}{\delta f_s}\right ]_{x,v}\! , \end{eqnarray}

where $\delta F/\delta f_s$ is the functional derivative of $F$ with respect to $f_s$ and

(2.7) \begin{align} \big [a,b\big ]_{x,v} = \boldsymbol{\nabla } a\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} b- \boldsymbol{\nabla } b\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} a. \end{align}

Vlasov equations (2.1) follow from

(2.8) \begin{eqnarray} \partial _t f_s = \{f_s,{\mathcal{H}}_{_{VP}}\}_{_{VP}}, \quad s=i,e , \end{eqnarray}

noticing that

(2.9) \begin{align} \frac {\delta {\mathcal{H}}_{_{VP}}}{\delta f_s}= \frac {1}{2}m_s v^2+q_s \phi , \end{align}

which is the total particle energy.

Finally, it is well known that the bracket (2.6) has the following infinite family of Casimir invariants, i.e. functionals $\mathcal{C}$ satisfying $\{F,\mathcal{C}\}=0,$ $\forall F$ :

(2.10) \begin{eqnarray} \mathcal{C}_s = \int {d}^3x{d}^3v\, \mathscr{C}_s(f_s), \end{eqnarray}

where $\mathscr{C}_s$ are arbitrary, well-behaved functions of $f_{s}$ .

2.2. The Vlasov--Ampère system

Another model used to describe electrostatic evolution of plasmas with negligible or totally vanishing magnetic field is the VA system. Now, the Vlasov equations are closed with the Ampère equation for the calculation of the electric field, instead of the Poisson equation (2.2), i.e. we have the following dynamical equation for ${\boldsymbol{E}}(\boldsymbol{x},t)$ :

(2.11) \begin{eqnarray} \epsilon _0 \partial _t{\boldsymbol{E}} = -\boldsymbol{J}, \end{eqnarray}

where

(2.12) \begin{eqnarray} \boldsymbol{J}(\boldsymbol{x},t) = \sum _{s} q_s \int {d}^3v\, f_s \boldsymbol{v} \end{eqnarray}

is the electric current density. The VA system can be formulated as a non-canonical Hamiltonian system with the Hamiltonian

(2.13) \begin{eqnarray} {\mathcal{H}}_{_{VA}}=\sum _{s}\int {d}^3x\,{d}^3v \, f_s m_s \frac {v^2}{2} + \frac {1}{2\epsilon _0}\int {d}^3x\, |{\boldsymbol{E}}|^2 \end{eqnarray}

and the bracket

(2.14) \begin{eqnarray} \{F,G\}_{_{VA}} = \sum _{s}\int {d}^3x\, {d}^3v\,\bigg \{ \frac {f_s}{m_s}\left [\frac {\delta F}{\delta f_s},\frac {\delta G}{\delta f_s}\right ]_{x,v}\nonumber \\ +\frac {q_s}{\epsilon _0 m_s}\left (\frac {\delta G}{\delta f_s}\frac {\delta F}{\delta {\boldsymbol{E}}} \boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} f_s -\frac {\delta F}{\delta f_s} \frac {\delta G}{\delta {\boldsymbol{E}}}\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}}f_s \right ) \bigg \}. \end{eqnarray}

The dynamical equations can be cast in the form of the Hamilton equations:

(2.15) \begin{align} \partial _t \boldsymbol{u} = \{\boldsymbol{u},{\mathcal{H}}_{_{VA}}\}_{_{VA}}, \quad \boldsymbol{u} = (f_s,E). \end{align}

It is important to note that the bracket (2.14) does not, in general, define a Poisson bracket, as it fails to satisfy the Jacobi identity. However, it does satisfy the identity and thus becomes a true Poisson bracket if the electric field remains irrotational, a condition required for the validity of the VA system, as mentioned in the introduction. This can be shown by following the proof of the Jacobi identity for the full VM bracket presented in Morrison (Reference Morrison2013) (see Appendix A). Note that this condition is trivially satisfied in one-dimensional (1D) plasmas, so the VA system and the corresponding Hamiltonian description are always valid in one dimension. Hence, as noted in the introduction, while we adopt the generic three-dimensional formalism in this subsection to maintain uniform notation with the analysis of the VP system, our treatment of the VA system ultimately concentrates on the 1D case.

The bracket (2.14) has two infinite families of Casimirs. One is given by (2.10) and the second one is

(2.16) \begin{eqnarray} {\mathcal{C}}_{{\boldsymbol{E}}} &\,=\,& \int {d}^3x\, g(\boldsymbol{x}) \left ( \epsilon _0 \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{E}} -\sum _{s}q_s \int {d}^3v\, f_{s} \right )\!, \end{eqnarray}

where $g(\boldsymbol{x})$ is an arbitrary function of $\boldsymbol{x}$ . Thus, in the VA model, the Poisson equation arises as a consequence of the conservation of ${\mathcal{C}}_{{\boldsymbol{E}}}$ , since ${\rm d}{\mathcal{C}}_{{\boldsymbol{E}}}/\text{d}t=0$ implies

(2.17) \begin{align} \partial _t \left ( \epsilon _0 \boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{E}} -\sum _{s}q_s \int {d}^3v\, f_{s} \right )=0. \end{align}

Thus, if the Gauss law is initially satisfied, it must be satisfied for all times.

3.1. The reformulated Vlasov–Poisson system

Quasineutrality, i.e. $n_i(\boldsymbol{x},t) = n_e(\boldsymbol{x},t)$ , arises naturally when considering large length scales, much larger than the electron Debye length, where electric fields are effectively screened by the particle response. This can be seen by introducing non-dimensional quantities:

(3.1) \begin{eqnarray} \tilde {x} = \frac {x}{L}, \quad \tilde {v}=\frac {v}{v_{th,e}}, \quad \tilde {t} = \frac {t}{L/v_{th,e}}, \quad f_{s} = \frac {f}{n_0/ v_{th,e}^3}, \quad \tilde {E} = \frac {E}{k_B T_e/eL}, \end{eqnarray}

where $L$ and $n_0$ are the characteristic length scale and density, respectively, and

(3.2) \begin{align} v_{th,e} = \sqrt {\frac {k_B T_e}{m_e}} \end{align}

is the electron thermal velocity. By this normalisation the VP system is written in the following non-dimensional form:

(3.3) \begin{eqnarray} \partial _t f_s + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla } f_{s} + \mu _{s} {\boldsymbol{E}}\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} f_{s}=0, \end{eqnarray}

(3.4) \begin{eqnarray} -\lambda ^2 \Delta \phi = n_i-n_e, \end{eqnarray}

where $\mu _i=\mu = m_e/m_i$ , $\mu _e = -1$ and $\lambda =\lambda _{D,e}/L$ , with

(3.5) \begin{align} \lambda _{D,e} = \sqrt {\frac {\epsilon _0 k_B T_e}{n_0 e^2}} \end{align}

being the electron Debye length. For length scales $L \gg \lambda _{D,e}$ , i.e. $\lambda \ll 1$ , Poisson equation (3.4) implies $n_i=n_e$ , where

(3.6) \begin{eqnarray} n_s = \int {d}^3v\, f_s, \quad s =i,e, \end{eqnarray}

are the particle densities of the two particle species.

A series of papers (Crispel et al. Reference Crispel, Degond and Vignal2005a , Reference Crispel, Degond and Vignalb ; Degond et al. Reference Degond, Deluzet and Navoret2006; Belaouar et al. Reference Belaouar, Crouseilles, Degond and Sonnendrücker2009; Degond et al. Reference Degond, Deluzet, Navoret, Sun and Vignal2010) considered the asymptotic limit $\lambda \rightarrow 0$ , where the Poisson equation can no longer close the system as it merely yields $n_i=n_e$ . To close the Vlasov equation, those works consider a reformulated VP system where the Poisson equation is replaced by another partial differential equation for the determination of the electrostatic potential, obtained by taking zeroth- and first-order moments of the Vlasov equation. Manipulating appropriately the moment equations, the following equation arises:

(3.7) \begin{eqnarray} -\boldsymbol{\nabla }\boldsymbol{\cdot }\left [ (n_e+\mu n_i+\lambda ^2\partial _{tt})\boldsymbol{\nabla } \phi \right ] = \boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{\colon } (\boldsymbol{P}_i-\boldsymbol{P}_e) , \end{eqnarray}

which allows the calculation of $\phi$ even in the limit $\lambda \rightarrow 0$ . In (3.7) the tensors $\boldsymbol{P}_s$ are defined as

(3.8) \begin{eqnarray} \boldsymbol{P}_s=\int {d}^3v\, \boldsymbol{v}\boldsymbol{v} f_s. \end{eqnarray}

3.2. The Dirac method of constraints

Here, we follow a different approach, exploiting the Hamiltonian formulation of the model to impose the quasineutrality condition as a Dirac constraint following the Dirac algorithm for the construction of a generalised Poisson bracket (Dirac Reference Dirac1950, Reference Dirac1958; Sundermeyer Reference Sundermeyer1982; Morrison et al. Reference Morrison, Andreussi and Pegoraro2020). The quasineutrality condition

(3.9) \begin{eqnarray} \int {d}^3v \, f_i = \int {d}^3v \, f_e \end{eqnarray}

can be seen as a constraint:

(3.10) \begin{eqnarray} \varPhi _1(\boldsymbol{x}) = \int {d}^3v\, \left [f_i(\boldsymbol{x},\boldsymbol{v},t)-f_e(\boldsymbol{x},\boldsymbol{v},t)\right ]=0. \end{eqnarray}

A consistency conditions is that the constraint function $\varPhi _1(\boldsymbol{x})$ should be preserved by the dynamics, so $\{\varPhi _1,{\mathcal{H}}\}\approx 0$ , where ‘ $\approx$ ’ denotes weak equality, i.e. the equality is satisfied on submanifolds of the phase space identified by the constraint. To impose this consistency condition, we need to write the constraint function $\varPhi _1$ as a phase-space functional, i.e.

(3.11) \begin{eqnarray} \varPhi _1(\boldsymbol{x})= \int {d}^3x' {d}^3v\, \delta (\boldsymbol{x}-\boldsymbol{x}')\left [f_i(\boldsymbol{x}',\boldsymbol{v},t)-f_e(\boldsymbol{x}',\boldsymbol{v},t)\right ]\!, \end{eqnarray}

so that the functional derivatives of $\varPhi _1$ are

(3.12) \begin{eqnarray} \frac {\delta \varPhi _1}{\delta f_i} = \delta (\boldsymbol{x}-\boldsymbol{x}'), \quad \frac {\delta \varPhi _1}{\delta f_e}=-\delta (\boldsymbol{x}-\boldsymbol{x}'). \end{eqnarray}

3.2.1. The Vlasov--Poisson system

In view of (3.12) and (2.6) we find

(3.13) \begin{eqnarray} \{\varPhi _1,{\mathcal{H}}\}_{_{VP}} = -\boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\, \boldsymbol{v}(f_i-f_e) = -\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}. \end{eqnarray}

Hence, the requirement $\{\varPhi _1,{\mathcal{H}}\}_{_{VP}}\approx 0$ results in a secondary constraint:

(3.14) \begin{eqnarray} \varPhi _2(\boldsymbol{x}) = \int {d}^3x'\, \delta (\boldsymbol{x}-\boldsymbol{x}') \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} = \int {d}^3x'\, \delta (\boldsymbol{x}-\boldsymbol{x}') \boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\,\boldsymbol{v}(f_i-f_e)=0,\qquad \end{eqnarray}

a condition that also arises in the analysis of Burby et al. (Reference Burby, Kaltsas, Morrison, Tassi and Throumoulopoulos2025). The consistency condition $\{\varPhi _2,{\mathcal{H}}\}_{_{VP}}\approx 0$ can verified that holds true in view of the two unconstrained Vlasov equations (2.1).

We can verify that $\varPhi _1$ and $\varPhi _2$ are second-class constraints, i.e. the quantity $\{\varPhi _1,\varPhi _2\}_{VP}$ is non-vanishing. So, we can form a non-singular, antisymmetric constraint matrix of the form

(3.15) \begin{eqnarray} C=\begin{pmatrix} C_{11} \quad& &\quad C_{12}\\ C_{21}\quad &&\quad C_{22} \end{pmatrix}\!, \end{eqnarray}

where $C_{jk}(\boldsymbol{x},\boldsymbol{x}') = \{\varPhi _j(\boldsymbol{x}),\varPhi _k(\boldsymbol{x}')\}_{VP}$ with $j,k=1,2$ . The details of the calculation of the elements $C_{jk}$ of the constraint matrix are presented in Appendix B. These elements are

(3.16) \begin{eqnarray} & C_{11}=0, \quad C_{12}(\boldsymbol{x},\boldsymbol{x}')=-C_{21}(\boldsymbol{x},\boldsymbol{x}') = {\mathcal{L}} \delta (\boldsymbol{x}'-\boldsymbol{x}),& \nonumber \\& C_{22}(\boldsymbol{x},\boldsymbol{x}')= \boldsymbol{\nabla }\boldsymbol{\cdot }\left [(\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\nabla } \delta (\boldsymbol{x}'-\boldsymbol{x})+\boldsymbol{M} \Delta \delta (\boldsymbol{x}'-\boldsymbol{x})+\boldsymbol{\nabla }\delta (\boldsymbol{x}'-\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{M}\right ]\!, & \qquad \end{eqnarray}

where

(3.17) \begin{eqnarray} {\mathcal{L}} := \boldsymbol{\nabla }\boldsymbol{\cdot }\left [ \left (\frac {n_i}{m_i}+\frac {n_e}{m_e}\right )\boldsymbol{\nabla } \right ] \end{eqnarray}

is a self-adjoint elliptic operator and

(3.18) \begin{eqnarray} \boldsymbol{M}(\boldsymbol{x},t) = \int {d}^3v\, \boldsymbol{v}\left (\frac {f_i}{m_i}+\frac {f_e}{m_e}\right )\!. \end{eqnarray}

In order to eliminate the second-class constraints, Dirac defined a new bracket algebra on the phase space such that the bracket of any phase- space function with a constraint vanishes, and thus the constraints become Casimir invariants of the new bracket. The Dirac bracket of two functionals $F$ and $G$ on the phase space is defined as follows:

(3.19) \begin{eqnarray} \{F,G\}^* = \{F,G\}- \int \int {d}^3x {d}^3x' \{F,\varPhi _j(\boldsymbol{x})\}C^{-1}_{jk}(\boldsymbol{x},\boldsymbol{x}')\{\varPhi _{k}(\boldsymbol{x}'),G\}, \end{eqnarray}

where $\{F,G\}$ is the standard Poisson bracket and $C^{-1}_{jk}$ are the inverse matrix elements calculated by

(3.20) \begin{eqnarray} \int {d}^3x'' C_{j\ell }(\boldsymbol{x},\boldsymbol{x}'')C^{-1}_{\ell k}(\boldsymbol{x}'',\boldsymbol{x}') = \delta _{jk} \delta (\boldsymbol{x}-\boldsymbol{x}'). \end{eqnarray}

After some straightforward manipulations we find

(3.21) \begin{align} &C_{22}^{-1}=0, \quad C_{12}^{-1}(\boldsymbol{x},\boldsymbol{x}^\prime ) = -C_{21}^{-1}(\boldsymbol{x},\boldsymbol{x}^\prime )=- {\mathcal{L}} ^{-1}\delta (\boldsymbol{x}-\boldsymbol{x}^\prime ), \nonumber \\ &C_{11}^{-1}(\boldsymbol{x},\boldsymbol{x}^{\prime} ) = {\mathcal{L}} ^{-1} {\mathcal{L}} ^{\prime-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\left [(\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\nabla }\delta (\boldsymbol{x}^{\prime} -\boldsymbol{x})+\boldsymbol{M}\Delta \delta (\boldsymbol{x}^{\prime} -\boldsymbol{x})+\boldsymbol{\nabla }\delta (\boldsymbol{x}^{\prime} -\boldsymbol{x})\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{M}\right ]\!,\end{align}

where $ {\mathcal{L}} ^{-1}$ is the inverse of $ {\mathcal{L}}$ . To derive these elements, we have used that if $ {\mathcal{L}}$ is a self-adjoint and invertible operator on $L^2$ (so that its inverse $\mathcal{L}^{-1}$ exists on $L^2$ ), then its inverse $ {\mathcal{L}} ^{-1}$ is also self-adjoint.

Substituting the elements ${\mathcal{C}}^{-1}_{jk}$ from (3.21) into the general form of the field-theoretic Dirac bracket (3.19), and after some tedious manipulations, we arrive at the following bracket:

(3.22) \begin{eqnarray} \{F,G\}^\star = \{F,G\}_{_{VP}}-\int {d}^3x\, \Big \{ \boldsymbol{\nabla }\left ( {\mathcal{L}} ^{-1} \varUpsilon _F\right ) \boldsymbol{M} \boldsymbol{:} \boldsymbol{\nabla }\boldsymbol{\nabla }\left ( {\mathcal{L}} ^{-1} \varUpsilon _G\right ) \nonumber \\ - \boldsymbol{\nabla }\left ( {\mathcal{L}} ^{-1} \varUpsilon _G\right ) \boldsymbol{M} \boldsymbol{:} \boldsymbol{\nabla }\boldsymbol{\nabla }\left ( {\mathcal{L}} ^{-1} \varUpsilon _F\right )+\varOmega _G {\mathcal{L}} ^{-1} \varUpsilon _F-\varUpsilon _G {\mathcal{L}} ^{-1}\varOmega _F \Big \}, \end{eqnarray}

where

(3.23) \begin{eqnarray} \varUpsilon _X &:=\, & \boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\, \left [ \frac {\delta X}{\delta f_i}\boldsymbol{\nabla }_{\boldsymbol{v}}\left ( \frac {f_i}{m_i}\right ) - \frac {\delta X}{\delta f_e}\boldsymbol{\nabla }_{\boldsymbol{v}}\left ( \frac {f_e}{m_e}\right )\right ]\!, \end{eqnarray}

(3.24) \begin{eqnarray} \varOmega _X &:=\,& \boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\, \boldsymbol{v}\left (\left [\frac {\delta X}{\delta f_i},\frac {f_i}{m_i}\right ]_{x,v}-\left [\frac {\delta X}{\delta f_e},\frac {f_e}{m_e}\right ]_{x,v}\right )\!, \quad X=F,G \end{eqnarray}

and $\{F,G\}_{_{VP}}$ is the standard Poisson bracket (2.6). One can show that ${\mathcal{C}}_s$ , $\varPhi _1$ and $\varPhi _2$ given by (2.10), (3.11) and (3.14), respectively, are Casimir invariants of the bracket (3.22). Hence, the second-class constraints have been incorporated into the dynamics of the system, and their conservation is guaranteed by the structure of the Dirac bracket (3.22), which replaces the original Poisson bracket (2.6).

3.2.2. The Vlasov--Ampère system

For the VA system we require that the constraint $\varPhi _1$ satisfies the consistency condition $\{\varPhi _1,{\mathcal{H}}\}_{_{VA}}\approx 0$ , where the Poisson bracket is now given by (2.14). We can easily show that, by this condition, the same secondary constraint (3.14) arises, and we can employ the same steps as in the VP case to derive the corresponding Dirac bracket in the VA case. After repeating the procedure we find that the Dirac bracket for the VA system is formally identical to (3.22) with the difference that $\{F,G\}_{VP}$ is replaced by $\{F,G\}_{VA}$ given by (2.14) and the quantities $\varOmega _{X}$ are now given by

(3.25) \begin{eqnarray} \varOmega _X &\ =\ & \boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\, \left [\boldsymbol{v}\left (\left [\frac {\delta X}{\delta f_i},\frac {f_i}{m_i}\right ]_{x,v}-\left [\frac {\delta X}{\delta f_e},\frac {f_e}{m_e}\right ]_{x,v}\right )+\frac {e}{\epsilon _0}\frac {\delta X}{\delta {\boldsymbol{E}}} \left (\frac {f_i}{m_i}+\frac {f_e}{m_e}\right )\right ]\!, \nonumber \\ X\ &=&\ F,G. \end{eqnarray}

Casimir invariants of the new bracket are ${\mathcal{C}}_s$ , the two second-class constraints $\varPhi _1$ and $\varPhi _2$ given by (2.10), (3.11) and (3.14), respectively, and

(3.26) \begin{eqnarray} {\mathcal{C}}_{{\boldsymbol{E}}} = \int {d}^3x\, g(\boldsymbol{x})\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{E}}, \end{eqnarray}

which is consistent with (2.16), since the quasineutrality condition has been imposed.

3.3. Constrained dynamics

For the VP system, the constrained dynamical equations that respect the preservation of the constraints $\varPhi _1$ and $\varPhi _2$ given by (3.11) and (3.14), respectively, are recovered from the Hamilton–Dirac equations:

(3.27) \begin{eqnarray} \partial _t f_s=\{f_s,{\mathcal{H}}_{_{VP}}\}^\star , \end{eqnarray}

where ${\mathcal{H}}_{_{VP}}$ is given by (2.5) and $\{F,G\}^\star$ is given by (3.22) with $\varUpsilon _X$ and $\varOmega _X$ specified in (3.23) and (3.24), respectively. For the VA system, the corresponding constrained equations are

(3.28) \begin{eqnarray} \partial _t f_s = \{f_s,{\mathcal{H}}_{_{VA}}\}^\star , \quad s =i,e \end{eqnarray}

(3.29) \begin{eqnarray} \partial _t {\boldsymbol{E}} = \{{\boldsymbol{E}},{\mathcal{H}}_{_{VA}}\}^\star , \end{eqnarray}

with ${\mathcal{H}}_{_{VA}}$ given by (2.13) and $\{F,G\}^\star$ by (3.22), where $\varUpsilon _X$ and $\varOmega _X$ are given by (3.23) and (3.25), respectively. Both (3.27) and (3.28) yield the following system of constrained Vlasov equations:

(3.30) \begin{eqnarray} \partial _t f_s &\,=\,& - \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla } f_s - \frac {q_s}{m_s}{\boldsymbol{E}}\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} f_s \pm \frac {1}{m_s}\bigg \{ \boldsymbol{\nabla }\alpha \boldsymbol{\cdot }\boldsymbol{\nabla } f_s -\left [(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\nabla } \alpha - \boldsymbol{\nabla } \beta + \boldsymbol{\nabla } \gamma \right ]\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} f_{s} \nonumber \\ &&+\,q_s\boldsymbol{\nabla }_{\boldsymbol{v}} f_s\boldsymbol{\cdot }\boldsymbol{\nabla } {\mathcal{L}} ^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\,{\boldsymbol{E}}\left (\frac {f_i}{m_i}+\frac {f_e}{m_e}\right ) \bigg \}, \end{eqnarray}

where $+$ corresponds to the ion equation and $-$ to the electron equation. In (3.30) the quantities $\alpha$ , $\beta$ , $\gamma$ are

(3.31) \begin{eqnarray} \alpha = {\mathcal{L}} ^{-1} \boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\, \boldsymbol{v}(f_i-f_e) = e^{-1} {\mathcal{L}} ^{-1} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}, \end{eqnarray}

(3.32) \begin{eqnarray} \beta = {\mathcal{L}} ^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\left [(\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\nabla }\alpha +\boldsymbol{M}\Delta \alpha +(\boldsymbol{\nabla }\alpha \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{M}\right ]\!, \end{eqnarray}

(3.33) \begin{eqnarray} \gamma = {\mathcal{L}} ^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\, \boldsymbol{v}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }(f_i-f_e)= {\mathcal{L}} ^{-1}\boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{:}(\boldsymbol{P}_i-\boldsymbol{P}_e) \,; \end{eqnarray}

hence, they can be computed upon solving the following elliptic equations:

(3.34) \begin{align} {\mathcal{L}} \alpha &= e^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J} , \end{align}

(3.35) \begin{align} {\mathcal{L}} \beta & = \boldsymbol{\nabla }\boldsymbol{\cdot }\left [(\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\nabla }\alpha +\boldsymbol{M}\Delta \alpha +(\boldsymbol{\nabla }\alpha \boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{M}\right ]\!,\end{align}

(3.36) \begin{align} {\mathcal{L}} \gamma &= \boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{:}(\boldsymbol{P}_i-\boldsymbol{P}_e).\end{align}

The electric field is given by ${\boldsymbol{E}}=-\boldsymbol{\nabla }\phi$ for the general VP system, while for the VA system, (3.29) yields

(3.37) \begin{eqnarray} \epsilon _0 \partial _t {\boldsymbol{E}} =- \boldsymbol{J} +e\left (\frac {n_i}{m_i}+\frac {n_e}{m_e}\right )\boldsymbol{\nabla }\alpha . \end{eqnarray}

Taking the divergence of this equation, yields $\partial _t (\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{E}}) = 0$ , in agreement with the Casimir invariant (3.26), since (3.34) implies that the two terms on the right-hand side of (3.37) form a divergence-free vector. Obviously, in one dimension, the condition $\boldsymbol{\nabla }\times {\boldsymbol{E}} = 0$ holds automatically, whereas in higher dimensions (3.37) implies that this condition is satisfied provided the current density $\boldsymbol{J}$ obeys $ \boldsymbol{\nabla } \times \boldsymbol{J} = e \, \boldsymbol{\nabla } \left (({n_i}/{m_i}) + ({n_e}/{m_e})\right ) \times \boldsymbol{\nabla } \alpha$ . However, as previously discussed, in dimensions greater than one the VA system must be replaced by the full VM system.

Now, let us consider the VP case, and further notice that upon substituting ${\boldsymbol{E}} = -\boldsymbol{\nabla } \phi$ into (3.30), the last term becomes

(3.38) \begin{eqnarray} &\pm\,& \frac {q_s}{m_s} (\boldsymbol{\nabla }_{\boldsymbol{v}} f_s)\boldsymbol{\cdot }\boldsymbol{\nabla } {\mathcal{L}} ^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\int {d}^3v\,{\boldsymbol{E}}\left (\frac {f_i}{m_i}+\frac {f_e}{m_e}\right )\nonumber\\ &=\,&\mp \frac {e} {m_s}\boldsymbol{\nabla }_{\boldsymbol{v}}f_s\boldsymbol{\cdot }\boldsymbol{\nabla } {\mathcal{L}} ^{-1} {\mathcal{L}} \phi = \mp \frac {q_s}{m_s} \boldsymbol{\nabla } \phi \boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}}f_s. \end{eqnarray}

Substituting this last equation in (3.30) we see that the electric field term is effectively eliminated from the constrained Vlasov equations, and therefore these two equations comprise a closed system, since the new advection fields

(3.39) \begin{eqnarray} \boldsymbol{\xi } :=\boldsymbol{\nabla }\alpha , \quad \boldsymbol{\zeta }:= \boldsymbol{\nabla } \beta , \quad \boldsymbol{\eta }:=\boldsymbol{\nabla }\gamma \end{eqnarray}

depend merely on $f_i$ and $f_e$ .

Note also that, along the same lines, we can find an analogous result for the 1D VA case, since the electric field, being a continuous function, can be written as the gradient of a differentiable function. Therefore, the final system of QN Vlasov equations, in both the VP and the 1D VA case, takes the form

(3.40) \begin{eqnarray} \partial _t f_s + \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla } f_s - \frac {q_s}{em_s}\left \{ \boldsymbol{\xi } \boldsymbol{\cdot }\boldsymbol{\nabla } f_s -\left [(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\xi } - \boldsymbol{\zeta } + \boldsymbol{\eta } \right ]\boldsymbol{\cdot }\boldsymbol{\nabla }_{\boldsymbol{v}} f_{s} \right \}=0, \quad s=i,e.\qquad\quad \end{eqnarray}

For the advection of the distribution functions $f_i$ and $f_e$ , we are interested in calculating the fields $\boldsymbol{\xi }$ , $\boldsymbol{\zeta }$ and $\boldsymbol{\eta }$ , which according to (3.34)–(3.36) are

(3.41) \begin{align}\boldsymbol{\nabla } \boldsymbol{\cdot }\left [(n_e+\mu n_i) \boldsymbol{\xi }\right ] &\ \ =\ \ \frac {m_e}{e}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{J}, \end{align}

(3.42) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left [(n_e+\mu n_i)\boldsymbol{\zeta }\right ] &\ \ =\ \ m_e \boldsymbol{\nabla }\boldsymbol{\cdot }\left [(\boldsymbol{M}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\xi }+\boldsymbol{M}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\xi }+(\boldsymbol{\xi }\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{M}\right ]\!, \end{align}

(3.43) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\left [(n_e+\mu n_i) \boldsymbol{\eta }\right ] &\ \ =\ \ m_e \boldsymbol{\nabla }\boldsymbol{\nabla }\boldsymbol{:}(\boldsymbol{P}_i-\boldsymbol{P}_e).\end{align}

Note that the indices $i$ and $e$ in the particle densities are retained, even though we have imposed $n_i=n_e$ as a Dirac constraint. This is because the Dirac algorithm incorporates the constraint into the dynamics by modifying the bracket structure, such that the constraint appears as a Casimir invariant of the form (3.11). As a result, the constrained dynamics does not explicitly enforce $n_i=n_e$ , but rather enforces the invariance of $\varPhi _1$ . This implies that the difference $n_i-n_e$ is conserved over time. Therefore, only if the system is initialised with $n_i-n_e=0$ (i.e. it satisfies quasineutrality), we can write $n_i=n_e=n$ .

In this framework, the Poisson equation of the VP system and the Ampère equation of the 1D VA system become irrelevant and superfluous, since the electric field does not participate in the advection of the distribution functions. The new advection fields are

(3.44) \begin{align} \boldsymbol{v} - \frac {q_s}{em_s}\boldsymbol{\xi } \end{align}

in physical (configuration) space and

(3.45) \begin{align} \frac {q_s}{em_s}[(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\xi }-\boldsymbol{\zeta }+\boldsymbol{\eta }] \end{align}

in velocity space.

4.1. The non-dimensional 1D-1V system

Let us now consider the QN Vlasov system in the 1D spatial and velocity coordinate setting (1D-1V), where the distribution function depends on a single spatial coordinate $x$ and a single velocity coordinate $v$ . We also introduce the following normalised quantities:

(4.1) \begin{eqnarray} \tilde {x} = \frac {x}{\lambda _{D,e}}, \quad \tilde {v} = \frac {v}{v_{th,e}}, \quad \tilde {t} = \omega _{p,e} t, \quad \tilde {f} = \frac {f}{n_0/v_{th,e}}, \nonumber \\ \quad \tilde {\xi } = \frac {\xi }{m_e v_{th,e}}, \quad \tilde {\zeta } = \frac {\zeta }{m_e v_{th,e}^2/\lambda _{D,e}}, \quad \tilde {\eta } = \frac {\eta }{m_e v_{th,e}^2/\lambda _{D,e}}, \end{eqnarray}

with

(4.2) \begin{align} \omega _{p,e}=\sqrt {\frac {n_0e^2}{\epsilon _0m_e}} \end{align}

being the electron plasma frequency, in order to write the 1D-1V counterparts of (3.40) in non-dimensional form:

(4.3) \begin{eqnarray} \partial _tf_s+v \partial _x f_s -\mu _s \left [ \xi \partial _x f_s-(v\partial _x\xi -\zeta +\eta ) \partial _vf_s\right ]=0, \quad s=i,e, \end{eqnarray}

and $\mu _i = \mu =m_e/m_i$ , $\mu _e =-1$ . The functions $\xi$ , $\zeta$ and $\eta$ satisfy

(4.4) \begin{eqnarray} &&\partial _x\left [(n_e+\mu n_i)\xi \right ] = \partial _x J, \end{eqnarray}

(4.5) \begin{eqnarray} && \partial _x\left [(n_e+\mu n_i)\zeta \right ] = \partial _x\left (2 M\partial _x\xi +\xi \partial _x M\right )\!, \end{eqnarray}

(4.6) \begin{eqnarray} && \partial _x\left [(n_e+\mu n_i)\eta \right ] = \partial _{xx} (P_i-P_e), \end{eqnarray}

where $J$ , $M$ and $P_s$ are the scalar 1D counterparts of $\boldsymbol{J}$ , $\boldsymbol{M}$ and $\boldsymbol{P}_s$ given by (2.12), (3.18) and (3.8), respectively. Integrating (4.4)–(4.6) with respect to $x$ and assuming that $n_e+\mu n_i \gt 0$ $\forall x$ , we find

(4.7) \begin{eqnarray} \xi = \frac {J+c_1}{n_e+\mu n_i} , \quad \zeta =\frac {2 M\partial _x\xi +\xi \partial _x M +c_2}{n_e+\mu n_i} , \quad \eta = \frac {\partial _x (P_i-P_e) +c_3}{n_e+\mu n_i}, \end{eqnarray}

where $c_1,c_2,c_3$ are integration constants. These constants must be independent of time in order to preserve the conservation properties of the Vlasov system. Their values may be fixed by enforcing suitable boundary conditions; however, for periodic boundary conditions these constants remain undetermined, so without loss of generality, we set $c_1=c_2=c_3=0$ , to avoid introducing potential non-physical shifts in the distribution function.

4.2. Solution with the semi-Lagrangian method

We solve (4.3) using a semi-Lagrangian method on a structured Eulerian grid $(x_j, v_k) = (j \Delta x, k \Delta v)$ , with $j = 0, 1, \dotsc , N_x$ and $k = 0, 1, \dotsc , N_v$ . The time integration is performed using a Strang splitting (Strang Reference Strang1968), similar to the classical Cheng–Knorr scheme (Cheng & Knorr Reference Cheng and Knorr1976; Sonnendrücker et al. Reference Sonnendrücker, Roche, Bertrand and Ghizzo1999).

To apply the splitting, we decompose the Vlasov equation for each species into a sequence of simpler subproblems. Each corresponds to a term of the operator $ {\mathcal{D}} _s$ in $\partial _t f_s + {\mathcal{D}} _s f_s=0$ that governs the transport of the distribution function $f_s$ . This operator is given by

(4.8) \begin{align} {\mathcal{D}} _s = {\mathcal{D}} _s^x + {\mathcal{D}} _s^{v}, \end{align}

with

(4.9) \begin{align} {\mathcal{D}} _s^x:=(v - \mu _s \xi )\partial _x, \quad {\mathcal{D}} _s^{v}:=\mu _s (\eta - \zeta +v\partial _x \xi ) \partial _v . \end{align}

This decomposition yields the following subsystems describing different physical effects:

(4.10) \begin{align} &\text{Advection in } x{:} \quad \partial _t f_s + {\mathcal{D}} _s^x f_s= 0, \end{align}

(4.11) \begin{align} &\text{Advection and shearing in } v{:} \quad \partial _t f_s + {\mathcal{D}} _s^{v} f_s= 0. \end{align}

In the semi-Lagrangian method each subsystem is solved along its corresponding characteristic equations:

(4.12) \begin{align} \frac {\text{d}X}{\text{d}t} &= (V-\mu _s\xi (X)), \end{align}

(4.13) \begin{align} \frac {\text{d}V}{\text{d}t} &= \mu _s(\eta (X)-\zeta (X)+V\partial _x \xi (X)). \end{align}

Discretising time as $t^n = n \Delta t$ , a full time step $\Delta t$ in the advection of $f_s$ according to (4.10) and (4.11) is approximated by the following Strang splitting:

(4.14) \begin{eqnarray} {\rm e}^{\Delta t ( {\mathcal{D}} _s^x + {\mathcal{D}} _s^{v})} \approx {\rm e}^{(\Delta t/2) {\mathcal{D}} _s^{x}} {\rm e}^{\Delta t {\mathcal{D}} _s^v} {\rm e}^{(\Delta t/2) {\mathcal{D}} _s^{x}}. \end{eqnarray}

Thus, both (4.12) and (4.13) are advanced independently, each as part of the composition of flows defined in (4.14).

In the semi-Lagrangian method we essentially determine the starting point of the characteristic curve $(X^n, V^n)$ ending at the grid point $(X^{n+1}, V^{n+1})$ , and set $f\!(X^{n+1}, V^{n+1}, t^{n+1}) = f\!(X^n, V^n, t^n)$ , since $f$ is conserved along characteristics. As $(X^n, V^n)$ generally does not lie on the grid, we use cubic interpolation from nearby points to evaluate $f\!(X^n, V^n)$ . Hence, in order to find the starting point $(X^n, V^n)$ we need to solve the characteristic (4.12)–(4.13) backwards in time. Notice that the characteristic (4.12) is a nonlinear ordinary differential equation for $X(t)$ and generally a nonlinear solver should be invoked, or a very small time step should be used. On the other hand, employing the method of integrating factor we can find an exact solution to (4.13):

(4.15) \begin{align} V(t + \Delta t) = V\!(t)\, {\rm e}^{\mu _s\, \partial _x \xi \, \Delta t}+ \frac {\eta -\zeta }{\partial _x \xi }\left ({\rm e}^{\mu _s\partial _x \xi \Delta t}-1\right )\!, \end{align}

so in the velocity step, we can compute the starting velocity as

(4.16) \begin{align} V^n = V^{n+1}{\rm e}^{-\mu _s \partial _x \xi \,\Delta t} - \frac {\eta -\zeta }{\partial _x \xi } \left (1-{\rm e}^{-\mu _s \partial _x \xi \,\Delta t}\right )\! . \end{align}

One can see that

(4.17) \begin{eqnarray} \lim _{\partial _x\xi \rightarrow 0} V^n = V^{n+1}- \mu _s(\eta -\zeta ) \Delta t. \end{eqnarray}

Hence, if $\partial _x\xi =0$ , the starting velocity is simply $V^n=V^{n+1}-\mu _s (\eta -\zeta ) \Delta t$ .

Thus, in practice, the Strang splitting (4.14) is employed by performing the following steps:

1. (i) Half-step advection in $x$ :

   (4.18) \begin{align} f_s'= f_s(x-(v-\mu _s\xi )\Delta t/2,v), \end{align}
   where $\xi$ is computed by $f_s$ .

2. (ii) Full-step advection and shearing in $v$ :

   (4.19) \begin{align} f_s''(x,v) = f_s'\left (x,v{\rm e}^{-\mu _s \partial _x \xi '\,\Delta t} - \frac {\eta '-\zeta '}{\partial _x \xi '} \left (1-{\rm e}^{-\mu _s \partial _x \xi ' \,\Delta t}\right ) \right )\!, \end{align}
   where $\xi '$ , $\zeta '$ and $\eta '$ are computed from $f_s'$ .

3. (iii) Half-step advection in $x$ :

   (4.20) \begin{align} f_s'''= f_s''\left (x-(v-\mu _s\xi '')\Delta t/2,v\right )\!, \end{align}
   where $\xi ''$ is computed by $f_s''$ .

To smooth out spurious oscillations in the particle density profile which may develop in regions where fine-scale structures emerge, we found it effective to apply a Savitzky–Golay filter (Schafer Reference Schafer2011) of polynomial order 2 to the distribution functions. This filter minimises the least-squares error, fitting a second-order polynomial to successive frames of high-frequency data. Thus, this procedure suppresses high-frequency oscillations without overly damping physical features.

The interpolation and filtering steps break energy conservation and the symplectic structure of the continuous dynamics, leading to the small but non-zero drifts of invariants observed in the next subsection. In addition, the standard Cheng–Knorr- type semi-Lagrangian splitting is employed here in an operational rather than a geometrical manner. Thus, although it is Hamiltonian for the ordinary VP system, it is not designed to preserve the Hamiltonian or the underlying Dirac bracket structure of the constrained QN system. To improve both conservation and long-time stability, an important direction for future work is the development of numerical schemes that respect the modified Hamiltonian structure induced by the Dirac bracket. This includes constructing energy-conserving or, more broadly, structure-preserving strategies within the semi-Lagrangian framework (e.g. Crouseilles et al. Reference Crouseilles, Mehrenberger and Sonnendrücker2010, Reference Crouseilles, Einkemmer and Faou2015; Liu et al. Reference Liu, Cai, Cao and Lapenta2023) that aim to retain the geometric properties of the continuous flow at the discrete level. The design of such structure-preserving algorithms for the constrained QN Vlasov model lies beyond the scope of the present study but will be pursued in subsequent work.

4.3. Simulation set-up

4.3.1. Electron–proton plasma ( $\mu \ll 1$ )

We initialise the simulation with two counter-streaming plasma beams having equal ion and electron densities, perturbed by a cosine modulation around the characteristic density $n_0=1$ and different drift velocities $V_e$ and $V_i$ :

(4.21) \begin{eqnarray} f_e = \frac {1}{2\sqrt {2\pi }} \left ({\rm e}^{-(v-V_e)^2/2}+{\rm e}^{-(v+V_e)^2/2}\right )\left [1+\epsilon \cos\left (\frac {2k\pi x}{L}+\pi \right )\right ] + \delta \frac {v{\rm e}^{-v^2/2}}{\sqrt {2\pi }},\qquad\qquad \end{eqnarray}

(4.22) \begin{eqnarray} f_i = \frac {1}{2\tau \sqrt {2\pi }} \left ({\rm e}^{-(v-V_i)^2/2\tau ^2}+{\rm e}^{-(v+V_i)^2/2\tau ^2}\right )\left [1+\epsilon \cos\left (\frac {2k\pi x}{L}+\pi \right )\right ]\!. \end{eqnarray}

Here, $\tau = v_{th,i}/v_{th,e}$ is the ratio of the ion to the electron thermal velocities; $\epsilon$ is a small perturbation parameter; $\delta$ controls the constant current density; $L$ is the size of the periodic box measured in electron Debye lengths $\lambda _{D,e}$ ; and $k$ is an integer indicating the number of density peaks within the box $L$ . Here we select $L=150$ , $\tau =1$ , $\epsilon =0.05$ , $\delta =0$ and $k=3$ . The drift velocities are $V_e=2$ and $V_i=0.5$ . To simulate the evolution of the distribution functions $f_i$ and $f_e$ we consider a $120\times 120$ grid discretising the simulation box, so that $\Delta x \sim \lambda _{D,e}$ while in velocity space the computational domain spans from $-4\pi$ to $4\pi$ , measured in electron thermal velocity units. The electron to ion mass ratio considered here is $\mu =1/1836$ , which corresponds to an electron–proton plasma. We also considered $\mu =1$ for an electron–positron plasma.

We simulate the time evolution of $f_i$ and $f_e$ up to 20 plasma periods $\omega _{p,e}^{-1}$ , while the time step is $\Delta t =10^{-4}$ . In figure 1 the initial conditions of the distribution functions $f_e$ and $f_i$ are plotted in phase space, while in figure 2 we provide snapshots at $t=5,10,15$ and $20$ $\omega _{p,e}^{-1}$ . A distinct evolution of $f_e$ is observed between the two scenarios, particularly during the nonlinear phase. In the QN case, phase-space vortices form earlier than in the VP simulation and have different shapes.

Figure 1. Initial conditions of the distribution functions $f_e$ and $f_i$ in phase space. The electron distribution forms two separated beams since $V_e$ is large enough in (4.21), whereas the two ion beams are so close to each other that they effectively form a single beam with zero macroscopic velocity. This broad, centralised ion beam was found to favour numerical stability over longer simulation times.

Figure 2. Snapshots of the electron and ion distribution functions in phase space $(x,v)$ at different times. The left-hand column corresponds to the VP case and the right-hand column to the QN case. A distinct evolution of $f_e$ is observed between the two scenarios, especially during the nonlinear phase.

Figure 3. Left: Evolution of the average charge density modulus $\langle |\rho |\rangle$ in the standard VP scenario (solid red line) versus the Dirac-constrained QN scenario (dashed blue line). The QN charge density, although non-zero, remains consistently three orders of magnitude smaller than in the VP case, even in the vortex saturation stage (right).

Figure 4. Left: Evolution of $\langle |\partial _x J|\rangle$ in the standard VP scenario (solid red line) versus the Dirac-constrained QN scenario (dashed blue line). Although non-zero, this quantity remains consistently at least three orders of magnitude smaller than in the VP case, even in the vortex saturation stage (right).

In figure 3, we show the temporal evolution of the space-averaged modulus of the charge density, $\langle |\rho | \rangle$ , as well as the ratio $\langle |\rho _{_{\mathrm{QN}}}| \rangle / \langle |\rho {_{_{\mathrm{VP}}}}| \rangle$ , while in figure 4 we show the corresponding diagrams for the second locally conserved quantity $\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{J}=\partial _x J$ . Both $\rho$ and $\partial _x J$ are three orders of magnitude smaller than their VP counterparts throughout the simulation. The initially vanishing charge density and the initial current incompressibility are not preserved to machine precision in the QN case, as the algorithm is not specifically designed to conserve energy and Casimir invariants with high accuracy. This limitation is also reflected in the relative energy and particle number errors in the VP simulation, shown in figure 5. These errors diminish with decreasing discretisation lengths, indicating that they primarily stem from discretisation and interpolation errors. Vlasov–Ampère simulations were also performed for the non-quasineutral system yielding similar results.

Figure 5. Relative energy error (left) and relative particle error (right) in the VP simulation. The particle number and energy are not conserved with high precision due to the non-conservative nature of the simulation algorithm, interpolation errors and the applied filtering method.

4.3.2. Electron–positron plasma ( $\mu =1$ )

For the case of an electron–positron plasma, where $\mu = 1$ , we assume that the distribution functions $f_e$ and $f_i$ are given by (4.21) and (4.22), respectively, introducing a shift by a phase $\pi$ in the $x$ direction in the positron distribution, $f_i$ , resulting in an initially non-vanishing charge density and a corresponding electric field. All other simulation parameters remain unchanged, except for the time step, which is set to $\Delta t=10^{-3}$ in this case. Thus, this set-up is intended to test the conservation of the Casimir $\varPhi _1$ , not to mimic a physically QN regime. It ultimately demonstrates that the constrained algorithm correctly handles non-zero initial charge distributions by enforcing their preservation, in contrast to the standard VP dynamics where the charge density evolves freely.

Figure 6. Snapshots of the electron and positron distribution functions in phase space $(x,v)$ at different times. The left column corresponds to the VP case and the right column to the QN case. We observe distinct evolution for both $f_e$ and $f_i$ .

Figure 7. Left: Evolution of the average charge density modulus $\langle |\rho |\rangle$ in the standard VP scenario (solid red line) versus the Dirac-constrained QN scenario (dashed blue line). We see that $\rho _{_{QN}}$ stays constant while $\rho _{_{VP}}$ shows significant change during the evolution. Right: The corresponding plots for the quantity $\langle |\partial _x J|\rangle$ showing that in the QN case the current incompressibility constraint is satisfied with acceptable accuracy.

In figure 6, we show contour plots of the distribution functions in phase space $(x, v)$ which reveal significant differences in the evolution of the distribution functions between the VP and the QN case. In figure 7, we confirm that the QN system approximately preserves the initial charge density, unlike the VP system, in which $\rho$ evolves and undergoes significant changes over time. The quantity $\partial _x J$ also exhibits significant temporal variations in the VP system, whereas in the QN case it remains approximately four orders of magnitude smaller. Thus, the current incompressibility constraint is satisfied with good precision.

The results presented in figures 3, 4 and 7 demonstrate that the numerical scheme preserves the key Casimir invariants imposed as Dirac constraints, namely the fixed charge density and the current incompressibility. This numerical preservation serves as a direct verification that the constrained dynamics has been implemented correctly. An additional verification step could, in principle, involve comparison with an analytical dispersion relation in the linear phase. However, deriving such a closed-form dispersion relation for the Dirac-constrained QN Vlasov system is significantly more involved than in the standard VP case. For this reason, a thorough dispersion-relation study, including comparison with theory, is deferred to work specifically focused on linear validation.

4.3.3. Estimating the significance of the Dirac forces

To estimate the significance of the Dirac forces we consider the case of electron–proton plasma ( $\mu \ll 1$ ) with immobile protons. To assess the relative strength of the Dirac forces acting on the electrons compared with the ‘fluid forces’ due to pressure and convection, we take the zeroth- and the first-order velocity moments of the electron Vlasov equation in the 1D-1V case to obtain

(4.23) \begin{eqnarray} \partial _t n_e &=& -\partial _x(n_eu_e)-\xi \partial _xn_e -n_e\partial _x\xi , \end{eqnarray}

(4.24) \begin{eqnarray} \partial _t(n_e u_e) &\,=\,& -\partial _x P_e-\xi \partial _x(n_eu_e)-2 n_eu_e\partial _x\xi +n_e(\zeta -\eta ), \end{eqnarray}

where $u_e = n_e^{-1} \int {d}^3v\, f_ev$ . Using (4.23) to reformulate (4.24) we find the following momentum equation:

(4.25) \begin{eqnarray} n_e \partial _t u_e= - n_e u_e \partial _x u_e -\partial _x \tilde P_e -n_e \xi \partial _xu_e -n_eu_e \partial _x\xi +n_e(\zeta -\eta ), \end{eqnarray}

where

(4.26) \begin{align} \tilde P_e = \int \text{d}v\, (v-u_e)^2f_e=P_e-n_eu_e^2. \end{align}

The first two terms on the right-hand side of (4.25) correspond to ‘fluid forces’ while the rest of the terms on the right-hand side are Dirac force densities. To quantify the relative strength of the Dirac forces compared with the fluid forces we investigate the evolution of the ratio

(4.27) \begin{eqnarray} \frac {\langle |n_e\xi \partial _xu_e +n_eu_e \partial _x\xi -n_e(\zeta -\eta )|\rangle }{\langle | n_eu_e \partial _x u_e +\partial _x \tilde P_e|\rangle } \end{eqnarray}

for various characteristic lengths $L$ , from $L=25\, \lambda _{De}$ to $L=250\,\lambda _{De}$ . For this comparison we consider a steady ion particle density $n_i=1-\epsilon \cos(2\pi k x/L)$ and the electrons are initially described by (4.21) with $V_e=0$ , in order to avoid the formation of phase-space vortices and the associated fine structures which are sources of instability and limit the simulation time window especially for small scales $L$ . The simulation ends at $t=10\omega _{pe}$ for which the quasineutrality condition is preserved with good accuracy. For later simulation times instability kicks in and quasineutrality is violated. In figure 8 we present the evolution of the ratio (4.27) for the various length scales $L$ and the time-averaged ratio. We observe that in all four cases $L=25$ , $L=50$ , $L=100$ and $L=250$ , the ratio of the Dirac forces over the fluid forces increases with time and becomes progressively smaller for larger length scales. The time-averaged ratio is over $10^{-1}$ for length scales $L\lt 50$ meaning that the Dirac forces responsible for imposing quasineutrality are significant and thus quasineutrality is not a good approximation while it becomes a better approximation for $L\gt 250$ where the averaged ratio is smaller than $10^{-2}$ . This comparison verifies a consistent QN scaling: Dirac forces become negligible at large system size, in agreement with the known asymptotic QN limit.

Figure 8. Left: The evolution of the ratio (4.27) of the Dirac forces over the fluid forces for a simulation with steady ion distribution. In all four cases $L=25$ , $L=50$ , $L=100$ and $L=250$ , this ratio increases with time and becomes progressively smaller for larger length scales. Right: The time-averaged ratio versus the length scale $L$ .