2.1.1. Scalar model

The scalar model, which was put forward in several earlier works (Marklund et al. Reference Marklund, Zamanian and Brodin2010; Zamanian et al. Reference Zamanian, Marklund and Brodin2010; Marklund & Morrison Reference Marklund and Morrison2011; Manfredi et al. Reference Manfredi, Hervieux and Hurst2019), is described by a scalar electron distribution function that depends on the extended phase space variable plus time,

(2.1)\begin{equation} f: (t,\boldsymbol{x},\boldsymbol{p},\boldsymbol{s})\in \mathbb{R}_+{\times} \mathbb{R}^3 \times \mathbb{R}^3 \times \mathbb{R}^3 \mapsto f(t,\boldsymbol{x},\boldsymbol{p},\boldsymbol{s}) \in \mathbb{R}, \end{equation}

together with the self-consistent electromagnetic fields $(\boldsymbol {E},\boldsymbol {B}): (t,\boldsymbol {x})\mapsto (\boldsymbol {E},\boldsymbol {B})(t,\boldsymbol {x})\in \mathbb {R}^3 \times \mathbb {R}^3$. These quantities obey the following Vlasov–Maxwell system (Crouseilles et al. Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021):

(2.2)\begin{gather} \left\{\!\begin{gathered} \frac{\partial f}{\partial t}+\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla} f+ \left[\boldsymbol{E}+\boldsymbol{p}\times\boldsymbol{B}+ {\beta}\mathfrak{h} \boldsymbol{\nabla}(\boldsymbol{s}\boldsymbol{\cdot}\boldsymbol{B})+ \gamma \mathfrak{h} \boldsymbol{\nabla} \left(\boldsymbol{B}\boldsymbol{\cdot} \frac{\partial}{\partial \boldsymbol{s}}\right)\right]\boldsymbol{\cdot} \frac{\partial f}{\partial\boldsymbol{p}}+\boldsymbol{s}\times\boldsymbol{B} \boldsymbol{\cdot} \frac{\partial f}{\partial\boldsymbol{s}}=0,\\ \frac{\partial\boldsymbol{E}}{\partial t}=\boldsymbol{\nabla}\times\boldsymbol{B}- \int_{\mathbb{R}^6} \boldsymbol{p}f \,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s} -{ \alpha} \mathfrak{h} \boldsymbol{\nabla} \times \int_{\mathbb{R}^6} \boldsymbol{s} {f} \,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s},\\ \frac{\partial\boldsymbol{B}}{\partial t} ={-} \boldsymbol{\nabla}\times\boldsymbol{E},\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E}={\rho-1},\quad \rho = \int_{\mathbb{R}^6} f \,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}=0, \end{gathered}\right. \end{gather}

with initial conditions

(2.2a–c)\begin{equation} f(t=0)=f^0,\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E}^0 = \int_{\mathbb{R}^6} f^0 \,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s}-1,\quad \boldsymbol{B}(t=0)=\boldsymbol{B}^0 \mbox{ such that } \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{B}^0=0, \end{equation}

where we have introduced the artificial parameters $\alpha, \beta, \gamma$ in order to facilitate the comparison with the vectorial model. Our earlier work (Crouseilles et al. Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021) used the above equations with $\alpha =1$, $\beta =1$ and $\gamma =0$ (i.e. no cross-derivative in $\boldsymbol {s}$ and $\boldsymbol {p}$).

2.1.2. Vector model

The vectorial model follows the standard representation of quantum mechanics in terms of density matrices, which are $2 \times 2$ matrices for spin-1/2 fermions. The corresponding Wigner function is also a $2 \times 2$ matrix $W(t, \boldsymbol {x}, \boldsymbol {p})\in {\mathcal {M}}_{2, 2}(\mathbb {C})$. The scalar distribution function is related to this matrix Wigner function by the formula (Manfredi et al. Reference Manfredi, Hervieux and Hurst2019)

(2.3)\begin{align} f(t, \boldsymbol{x}, \boldsymbol{p}, \boldsymbol{s}) & = \frac{1}{2{\rm \pi}}\sum_{j, k=1}^2 \frac{1}{2}[\delta_{j,k} + \boldsymbol{s}\boldsymbol{\cdot} {\boldsymbol \sigma}_{j,k}]W_{k, j} (t, \boldsymbol{x}, \boldsymbol{p}) \nonumber\\ & = \frac{1}{4{\rm \pi}}(1+s_3)W_{1, 1}(t, \boldsymbol{x}, \boldsymbol{p}) + \frac{1}{4{\rm \pi}}(s_1 -is_2)W_{2, 1}(t, \boldsymbol{x}, \boldsymbol{p}) \nonumber\\ & \quad + \frac{1}{4{\rm \pi}}(s_1+is_2)W_{1, 2}(t, \boldsymbol{x}, \boldsymbol{p})+\frac{1}{4{\rm \pi}}(1-s_3)W_{2, 2} (t, \boldsymbol{x}, \boldsymbol{p}), \end{align}

where ${\boldsymbol \sigma }=(\sigma _x, \sigma _y, \sigma _z)$ denote the Pauli matrices, $\boldsymbol {s}=(s_1, s_2, s_3)\in {\mathbb {S}^2}$ (${\mathbb {S}^2}$ being the sphere in ${\mathbb {R}^3}$), and $W_{\beta, \alpha }$ are the components of the Wigner function matrix (Zamanian et al. Reference Zamanian, Marklund and Brodin2010). Note that the scalar distribution is a linear function of the spin variable $\boldsymbol s$. Hence, one can write

(2.4)\begin{equation} \begin{pmatrix} W_{1, 1}(t, \boldsymbol{x}, \boldsymbol{p}) & W_{ 1, 2}(t, \boldsymbol{x}, \boldsymbol{p}) \\ W_{2, 1}(t, \boldsymbol{x}, \boldsymbol{p}) & W_{ 2, 2}(t, \boldsymbol{x}, \boldsymbol{p}) \end{pmatrix}=\int_{\mathbb{S}^2} f(t, \boldsymbol{x}, \boldsymbol{p},\boldsymbol{s})\frac{1}{2} \begin{pmatrix} 1+3s_3 & 3(s_1-{\rm i}s_2) \\ 3(s_1+{\rm i}s_2) & 1-3s_3 \end{pmatrix} \mathrm{d}\boldsymbol{s}. \end{equation}

Transforming to spherical coordinates on ${\mathbb {S}^2}$, $s_1=\sin \theta \cos \varphi, s_2=\sin \theta \sin \varphi, s_3=\cos \theta$ with $\theta \in [0, {\rm \pi}], \varphi \in [0, 2{\rm \pi} ]$, we have

(2.5a–d)\begin{equation} \int_{{\mathbb{S}}^2} s_i\,\mathrm{d}{\boldsymbol{s}} = 0,\quad \int_{{\mathbb{S}}^2} s_is_j\,\mathrm{d}{\boldsymbol{s}} =\delta_{i,j} \frac{4{\rm \pi}}{3},\quad \int_{{\mathbb{S}}^2} s_is_js_k\,\mathrm{d}{\boldsymbol{s}} =0,\quad \int_{{\mathbb{S}}^2} \mathrm{d}{\boldsymbol{s}} =4{\rm \pi}. \end{equation}

Then, it is possible to make the link between the scalar function $f(t, \boldsymbol {x}, \boldsymbol {p}, \boldsymbol {s})$ and its zeroth- and first-order ‘spin moments’ $(f_0, {\boldsymbol {f}})(t, \boldsymbol {x}, \boldsymbol {p})$, defined by

(2.6a,b)\begin{equation} f_0 = \int_{\mathbb{S}^2} f \,\mathrm{d}{\boldsymbol{s}},\quad {\boldsymbol{f}} = 3 \int_{\mathbb{S}^2} \boldsymbol{s} f\,\mathrm{d}{\boldsymbol{s}}. \end{equation}

Combining (2.3), (2.4) and (2.6a,b), we obtain the following relation between the scalar and vector representations:

(2.7)\begin{equation} f(t, \boldsymbol{x}, \boldsymbol{p},{\boldsymbol{s}})= \frac{1}{4{\rm \pi}} (f_0(t, \boldsymbol{x}, \boldsymbol{p}) + {\boldsymbol{s}}\boldsymbol{\cdot} {\boldsymbol{f}}(t, \boldsymbol{x}, \boldsymbol{p})). \end{equation}

Again, we note the linear relationship in the spin variable. We also note that $(f_0, {\boldsymbol {f}})$ are nothing but the components of the Wigner function $W_{\alpha, \beta }$ written in the Pauli basis, i.e. $f_0 = {\rm Tr} (W)$, $\boldsymbol {f} = {\rm Tr} ({\boldsymbol \sigma } W)$.

Then, the evolution equations for the vectorial model can be seen as the evolution equations for the ‘spin moments’ of the scalar distribution $f(t, \boldsymbol {x}, \boldsymbol {p},{\boldsymbol {s}})$, which are obtained by integrating (2.2) in spin space and assuming the relationship (2.7). This yields

(2.8)\begin{equation} \left\{\begin{gathered} \frac{\partial f_0}{\partial t}+\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla} f_0+ \left(\boldsymbol{E}+\boldsymbol{p} \times\boldsymbol{B}\right) \boldsymbol{\cdot} \frac{\partial f_0}{\partial\boldsymbol{p}}+ \frac{({\beta}+2\gamma)\mathfrak{h}}{3} \sum_{j=1}^3\boldsymbol{\nabla}{B_j}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}} {f_j}=0,\\ \frac{\partial f_j}{\partial t}+\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla} f_j+ \left(\boldsymbol{E}+\boldsymbol{p} \times\boldsymbol{B}\right) \boldsymbol{\cdot} \frac{\partial f_j}{\partial\boldsymbol{p}}+{{ \beta}}\mathfrak{h} \boldsymbol{\nabla}{B_j}\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{p}} {f_0}+({\boldsymbol{B}}\times {\boldsymbol{f}})_j=0, \quad j = 1, 2, 3,\\ \frac{\partial\boldsymbol{E}}{\partial t}=\boldsymbol{\nabla}\times\boldsymbol{B}-\left(\int_{\mathbb{R}^3} \boldsymbol{p}f_0 \mathrm{d}\boldsymbol{p} +\frac{{\alpha}\mathfrak{h}}{3} \boldsymbol{\nabla} \times \int_{\mathbb{R}^3} {\boldsymbol{f}} \mathrm{d}\boldsymbol{p}\right),\\ \frac{\partial\boldsymbol{B}}{\partial t} ={-} \boldsymbol{\nabla}\times\boldsymbol{E},\\ \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E}={\rho-1},\quad \rho = \int_{\mathbb{R}^3} f_0 \,\mathrm{d}\boldsymbol{p},\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}=0, \end{gathered}\right. \end{equation}

with initial conditions

(2.9a–d)\begin{equation} \left.\begin{gathered} (f_0(t=0),\quad \boldsymbol{f}(t=0))=(f_0^0,\boldsymbol{f}^0),\quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{E}^0=\int_{\mathbb{R}^3} f_0^0 \,\mathrm{d}\boldsymbol{p}-1,\\ \boldsymbol{B}(t=0)=\boldsymbol{B}^0 \mbox{ such that } \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B}^0=0. \end{gathered}\right\} \end{equation}

However, as we shall see, the structure (2.7) is not satisfied for all times by the solution of the scalar model (2.2) for any arbitrary choice of $\alpha, \beta, \gamma$. When this structure is not preserved by the time evolution, then the equivalence between the scalar model and its spin moments (which constitute the vector model) is not satisfied. Only for a particular choice of the parameters can one establish the equivalence, as stated in the following Proposition.

Proof. Assuming that the initial condition satisfies $f(t=0) = ({1}/{4{\rm \pi} }) (f_0(t=0) + \boldsymbol {s}\boldsymbol {\cdot } \boldsymbol {f}(t=0))$, we check that $f(t) =({1}/{4{\rm \pi} })(f_0(t)+{\boldsymbol {s}}\boldsymbol {\cdot }{\boldsymbol {f}}(t))=({1}/{4{\rm \pi} })(f_0(t)+\sum _i s_if_i(t))$ is true for all time $t>0$. To do so, we insert the decomposition $f(t) =({1}/{4{\rm \pi} })(f_0(t)+{\boldsymbol {s}}\boldsymbol {\cdot }{\boldsymbol {f}}(t))$ into the scalar equation (2.2). After some calculations detailed in Appendix A, we obtain

(2.10)\begin{align} 0 & = \frac{1}{4{\rm \pi}}\frac{\partial (f_0+{\boldsymbol{s}}\boldsymbol{\cdot}{\boldsymbol{f}})}{\partial t}+ \frac{1}{4{\rm \pi}}\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla} (f_0+{\boldsymbol{s}}\boldsymbol{\cdot}{\boldsymbol{f}})+ \frac{1}{4{\rm \pi}}(\boldsymbol{s}\times\boldsymbol{B}) \boldsymbol{\cdot} \frac{\partial (f_0+{\boldsymbol{s}}\boldsymbol{\cdot}{\boldsymbol{f}})}{\partial\boldsymbol{s}} \nonumber\\ & \quad +\frac{1}{4{\rm \pi}}\left[\left(\boldsymbol{E}+\boldsymbol{p}\times\boldsymbol{B}\right)+ \beta\mathfrak{h}\boldsymbol{\nabla}(\boldsymbol{s}\boldsymbol{\cdot}\boldsymbol{B})+ \beta\mathfrak{h} \boldsymbol{\nabla} \left(\boldsymbol{B}\boldsymbol{\cdot} \frac{\partial}{\partial \boldsymbol{s}}\right)\right]\boldsymbol{\cdot} \frac{\partial (f_0+{\boldsymbol{s}}\boldsymbol{\cdot}{\boldsymbol{f}})}{\partial\boldsymbol{p}}\nonumber\\ & =\frac{1}{4{\rm \pi}}\left(\frac{\partial f_0}{\partial t}+\boldsymbol{p} \boldsymbol{\cdot}\boldsymbol{\nabla} f_0+\left(\boldsymbol{E}+\boldsymbol{p}\times\boldsymbol{B}\right) \boldsymbol{\cdot}\frac{\partial f_0}{\partial\boldsymbol{p}}+\beta \mathfrak{h}\sum_{i} \boldsymbol{\nabla} {B}_i \boldsymbol{\cdot}\frac{\partial f_i}{\partial\boldsymbol{p}}\right)\nonumber\\ & \quad +\frac{1}{4{\rm \pi}}\sum_{i} s_i \left(\frac{\partial f_i}{\partial t}+ \boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{\nabla} f_i+ \left(\boldsymbol{E}+\boldsymbol{p}\times\boldsymbol{B}\right) \boldsymbol{\cdot}\frac{\partial f_i}{\partial\boldsymbol{p}}+ (\boldsymbol{B}\times\boldsymbol{f})_i+\beta \mathfrak{h} \boldsymbol{\nabla} {B}_i \cdot\frac{\partial f_0}{\partial\boldsymbol{p}} \right). \end{align}

Since $(1, s_1, s_2, s_3)$ is a basis, we can identify the coefficients to be zero which leads to the four equations of the vector model satisfied by $(f_0, f_1, f_2, f_3)$. Thus, since $f$ and $({1}/{4{\rm \pi} }) (f_0+{\boldsymbol {s}}\boldsymbol {\cdot }{\boldsymbol {f}})$ share the same initial condition, we get that $f=({1}/{4{\rm \pi} })(f_0+{\boldsymbol {s}}\boldsymbol {\cdot }{\boldsymbol {f}})$ is a solution of the scalar model if and only if $(f_0,f_1,f_2,f_3)$ is the solution of vector model (2.8).

This Proposition provides a link between the vector and the scalar models, which implies that we can solve either a six-dimensional phase space vector model or an extended nine-dimensional scalar model. Let us remark that the additional term (preceded by the coefficient $\gamma$) is not easy to approximate using PIC methods, since it is a second-order cross-derivative term, whereas its counterpart in the vector model is a simple transport term.

2.1.3. Poisson structure

It has been shown in (Marklund & Morrison Reference Marklund and Morrison2011) (see also Qin et al. Reference Qin, Liu, Xiao, Zhang, He, Wang, Sun, Burby, Ellison and Zhou2015; Xiao et al. Reference Xiao, Qin, Liu, He, Zhang and Sun2015) that the scalar model (2.2) with $\gamma =0,\ \alpha =\beta$ enjoys a non-canonical Poisson structure with the following Poisson bracket:

(2.11)\begin{align} \{\mathcal{F},\mathcal{G}\} & = \int f\left[\frac{\delta\mathcal{F}}{\delta f},\frac{\delta\mathcal{G}} {\delta f}\right]_{\boldsymbol{xp}}\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p}\,{\rm d}\boldsymbol{s} +\int\left(\frac{\delta\mathcal{F}}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \frac{\partial f}{\partial\boldsymbol{p}}\frac{\delta\mathcal{G}}{\delta f}- \frac{\delta\mathcal{G}}{\delta\boldsymbol{E}}\boldsymbol{\cdot}\frac{\partial f}{\partial\boldsymbol{p}} \frac{\delta\mathcal{F}}{\delta f}\right)\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s} \nonumber\\ & \quad +\int\left(\frac{\delta\mathcal{F}}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \left(\triangledown\times\frac{\delta\mathcal{G}}{\delta\boldsymbol{B}}\right)- \frac{\delta\mathcal{G}}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \left(\triangledown\times\frac{\delta\mathcal{F}}{\delta\boldsymbol{B}}\right)\right)\mathrm{d}\boldsymbol{x} \nonumber\\ & \quad +\int f\boldsymbol{B}\boldsymbol{\cdot}\left(\frac{\partial}{\partial\boldsymbol{p}} \frac{\delta\mathcal{F}}{\delta f}\times \frac{\partial}{\partial\boldsymbol{p}}\frac{\delta\mathcal{G}}{\delta f}\right)\mathrm{d}\boldsymbol{x}\, \mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s} \nonumber\\ & \quad + {\frac{1}{\beta\mathfrak{h}}} \int f\boldsymbol{s} \boldsymbol{\cdot} \left( \frac{\partial }{\partial \boldsymbol{s}} \frac{\delta \mathcal{F}}{\delta f} \times \frac{\partial }{\partial \boldsymbol{s}} \frac{\delta \mathcal{G}}{\delta f}\right)\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p}\, \mathrm{d}\boldsymbol{s}, \end{align}

and the Hamiltonian functional

(2.12)\begin{equation} \mathcal{H}(f,\boldsymbol{E},\boldsymbol{B})=\frac{1}{2}\int \boldsymbol{|p|}^{2}f \,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s}- {\beta\mathfrak{h}}\int \boldsymbol{s}\boldsymbol{\cdot} \boldsymbol{B} f \,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p}\,\mathrm{d}\boldsymbol{s}+ \frac{1}{2}\int \left(|\boldsymbol{E}|^{2}+|\boldsymbol{B}|^{2}\right)\mathrm{d}\boldsymbol{x}, \end{equation}

so that the scalar model can be reformulated as $\partial _t {\mathcal {Z}} = \{{\mathcal {Z}}, {\mathcal {H}}\}$, with ${\mathcal {Z}} = (f, \boldsymbol {E}, \boldsymbol {B})$.

We can also construct a geometric Poisson structure for the vector model introduced in Hurst et al. (Reference Hurst, Morandi, Manfredi and Hervieux2014, Reference Hurst, Morandi, Manfredi and Hervieux2017). For that, we need $\beta +2\gamma =\alpha$ in (2.8) and the Poisson bracket is

(2.13)\begin{align} \{\mathcal{F},\mathcal{G}\} & = \int {f_0} \left[\frac{\delta \mathcal{F}}{\delta f_0}, \frac{\delta \mathcal{G}}{\delta f_0}\right]\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} + \sum_{i=1}^3\int {f_i} \left[\frac{\delta \mathcal{F}}{\delta f_i},\frac{\delta \mathcal{G}} {\delta f_0}\right]\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} \nonumber\\ & \quad +\sum_{i=1}^3\int {f_i} \left[\frac{\delta \mathcal{F}}{\delta f_0},\frac{\delta \mathcal{G}} {\delta f_i}\right]\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} + \sum_{i=1}^3 {\frac{3\beta}{\alpha}} \int {f_0} \left[\frac{\delta \mathcal{F}}{\delta f_i},\frac{\delta \mathcal{G}} {\delta f_i}\right]\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} \nonumber\\ & \quad + \sum_{i=0}^3 \int \left(\frac{\delta \mathcal{F}}{\delta \boldsymbol{E}}\boldsymbol{\cdot} \frac{\partial f_i}{\partial \boldsymbol{p}}\frac{\delta \mathcal{G}}{\delta f_i}- \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}}\boldsymbol{\cdot} \frac{\partial f_i}{\partial \boldsymbol{p}} \frac{\delta \mathcal{F}}{\delta f_i}\right)\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} \nonumber\\ & \quad +\int f_0\boldsymbol{B}\boldsymbol{\cdot} \left(\frac{\partial}{\partial\boldsymbol{p}} \frac{\delta\mathcal{F}}{\delta f_0}\times \frac{\partial}{\partial\boldsymbol{p}}\frac{\delta\mathcal{G}}{\delta f_0}\right)\mathrm{d}\boldsymbol{x}\, \mathrm{d}\boldsymbol{p} \nonumber\\ & \quad +\sum_{i=1}^3\int f_i\boldsymbol{B}\boldsymbol{\cdot}\left(\frac{\partial}{\partial\boldsymbol{p}} \frac{\delta\mathcal{F}}{\delta f_i}\times \frac{\partial}{\partial\boldsymbol{p}}\frac{\delta\mathcal{G}}{\delta f_0} -\frac{\partial}{\partial\boldsymbol{p}}\frac{\delta\mathcal{F}}{\delta f_0}\times \frac{\partial}{\partial\boldsymbol{p}}\frac{\delta\mathcal{G}}{\delta f_i}\right)\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} \nonumber\\ & \quad +{\frac{3}{\alpha}}\frac{1}{\mathfrak{h}}\int \boldsymbol{f} \boldsymbol{\cdot} \left(\frac{\delta \mathcal{F}}{\delta \boldsymbol{f}} \times \frac{\delta \mathcal{G}}{\delta \boldsymbol{f}}\right)\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p}\nonumber\\ & \quad + \int\left(\frac{\delta\mathcal{F}}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \left(\triangledown\times\frac{\delta\mathcal{G}}{\delta\boldsymbol{B}}\right) -\frac{\delta\mathcal{G}}{\delta\boldsymbol{E}}\boldsymbol{\cdot} \left(\triangledown\times\frac{\delta\mathcal{F}}{\delta\boldsymbol{B}}\right)\right)\mathrm{d}\boldsymbol{x}, \end{align}

while the Hamiltonian functional is

(2.14)\begin{equation} \mathcal{H}(f_0, \boldsymbol{f}, \boldsymbol{E}, \boldsymbol{B}) = \frac{1}{2}\int |\boldsymbol{p}|^2 f_0 \,\mathrm{d}\boldsymbol{x}\,\mathrm{d}\boldsymbol{p} + \frac{1}{2}\int \left(|\boldsymbol{E}|^2 + |\boldsymbol{B}|^2\right) \mathrm{d}\boldsymbol{x} - {\frac{\alpha}{3}}\mathfrak{h} \int \boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{f} \,\mathrm{d} \boldsymbol{x}\,\mathrm{d}\boldsymbol{p}, \end{equation}

so that the system can be written as $\partial _t \mathcal {Z} = \{ \mathcal {Z}, \mathcal {H} \}$. It is easy to check that this bracket is bilinear, skew-symmetric and verifies Leibniz's rule, but it is not clear whether the Jacobi identity is satisfied. Hence, this bracket is not strictly speaking a Poisson bracket such as the one occurring in the scalar model; nevertheless, we will still refer to it as a Poisson bracket for the sake of simplicity. We note that the case $\beta =1, \gamma =1$ and $\alpha =3$ yields the vector model introduced in (Hurst et al. Reference Hurst, Morandi, Manfredi and Hervieux2014, Reference Hurst, Morandi, Manfredi and Hervieux2017; Manfredi et al. Reference Manfredi, Hervieux and Crouseilles2022) and ensures, by Proposition 1, the equivalence between the two models.

2.1.4. Summary

In summary, we have found that, in order to ensure the equivalence between the scalar model (2.2) and the vector model (2.8), one needs to take $\beta = \gamma$ and arbitrary $\alpha$. Furthermore, the vector model enjoys a geometric (antisymmetric-bracket) structure when $\beta + 2\gamma = \alpha$. The choice $\beta =1, \gamma =1$ and $\alpha =3$ satisfies both conditions, and recovers the ‘standard’ vector model commonly used in the literature (Hurst et al. Reference Hurst, Morandi, Manfredi and Hervieux2014, Reference Hurst, Morandi, Manfredi and Hervieux2017; Manfredi et al. Reference Manfredi, Hervieux and Crouseilles2022). However, as the $\gamma$-term is difficult to implement in a PIC method because of the double derivative, it has been common to take $\gamma =0$ when simulating the scalar model (Crouseilles et al. Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021). In that case, it is still possible to recover the Poisson structure when $\alpha =\beta$, but the equivalence with the standard vectorial model no longer holds.

2.2.1. The 1-D scalar model

Following Ghizzo et al. (Reference Ghizzo, Bertrand, Shoucri, Johnston, Fualkow and Feix1990), Crouseilles et al. (Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021) and Manfredi et al. (Reference Manfredi, Hervieux and Crouseilles2022), one can derive a reduced 1-D spin-Vlasov–Maxwell model by considering the case of a plasma interacting with an electromagnetic wave propagating in the longitudinal $x$ direction and assuming that all fields depend spatially on $x$ only. Choosing the Coulomb gauge ${\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol A =0}$, the vector potential $\boldsymbol {A}$ lies in the perpendicular (transverse) plane, i.e. $\boldsymbol {A} = (0, A_y, A_z) = (0, \boldsymbol {A}_\perp)$. Using $\boldsymbol {E} = -\boldsymbol {\nabla } \phi - \partial _t \boldsymbol {A}$, we obtain, using the notation $\boldsymbol {E} = (E_x, E_y, E_z) =(E_x, \boldsymbol {E}_\perp )$: $\boldsymbol {E}_\perp = - \partial _t \boldsymbol {A}_\perp \mbox { and } E_x=-\partial _x \phi$, and $\boldsymbol {B} =\boldsymbol {\nabla }\times \boldsymbol {A} = (0,- \partial _x A_z, \partial _x A_y)$. We then consider an electron distribution that is ‘cold’ in the transverse plane, i.e. $\delta (\boldsymbol {p}_\perp - \boldsymbol {A}_\perp ) f(t, x, p_x, \boldsymbol {s})$, where $\boldsymbol {p}=(p_x,p_y,p_z)=(p_x,\boldsymbol {p}_\perp )$ is the linear momentum and $\boldsymbol {p}- \boldsymbol {A}$ is the canonical momentum. After integration with respect to $\boldsymbol {p}_\perp$, the relevant extended phase space is reduced to five dimensions, instead of nine dimensions for the general case. in the following, the longitudinal variable $p_x$ will be simply denoted by $p$.

The scalar spin-Vlasov–Maxwell system (2.2) becomes

(2.15)\begin{equation} \left\{\begin{gathered} \frac{\partial f}{\partial t} + p \frac{\partial f}{\partial x} + \left[E_x - \boldsymbol{A}_\perp \boldsymbol{\cdot} \frac{\partial \boldsymbol{A}_\perp}{\partial x} - \mathfrak{h} \frac{\partial^2 A_z}{\partial x^2}\left(\beta s_2+\gamma \frac{\partial }{\partial s_2}\right) + \mathfrak{h} \frac{\partial^2 A_y}{\partial x^2} \left(\beta s_3+\gamma \frac{\partial }{\partial s_3}\right)\right] \frac{\partial f}{\partial p} \\ \quad + \left[s_3 \frac{\partial A_z}{\partial x} + s_2 \frac{\partial A_y}{\partial x}, -s_1 \frac{\partial A_y}{\partial x}, -s_1 \frac{\partial A_z}{\partial x}\right] \cdot \frac{\partial f}{\partial \boldsymbol{s}} = 0,\\ \frac{\partial E_x}{\partial t} ={-}\int_{\mathbb{R}^4} p f \,\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s},\\ \frac{\partial E_y}{\partial t} ={-} \frac{\partial^2 A_y}{\partial x^2} + A_y \int_{\mathbb{R}^4} f \,\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s} + \alpha \mathfrak{h}\int_{\mathbb{R}^4} s_3 \frac{\partial f}{\partial x}\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s},\\ \frac{\partial E_z}{\partial t} ={-} \frac{\partial^2 A_z}{\partial x^2} + A_z \int_{\mathbb{R}^4} f \,\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s} - \alpha \mathfrak{h}\int_{\mathbb{R}^4} s_2 \frac{\partial f}{\partial x}\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s},\\ \frac{\partial \boldsymbol{A}_\perp}{\partial t} ={-} \boldsymbol{E}_\perp,\\ \frac{\partial E_x}{\partial x} = \int_{\mathbb{R}^4} f \,\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s} - 1, \end{gathered}\right. \end{equation}

with initial condition

(2.16a–d)\begin{equation} \left.\begin{gathered} f(t=0)=f^0,\quad \frac{\partial E_x^0}{\partial x} = \int_{\mathbb{R}^4} f^0 \,\mathrm{d}{p}\, \mathrm{d}\boldsymbol{s} - 1,\\ (E_y, E_z)(t=0)=(E_y^0, E_z^0),\quad \boldsymbol{A}_\perp(t=0)=\boldsymbol{A}_\perp^0. \end{gathered}\right\} \end{equation}

Taking $\gamma =0$ and $\alpha =\beta$, the scalar spin-Vlasov–Maxwell system (2.15) can be expressed with one Poisson bracket as follows. For any two functionals $\mathcal {F}$, and $\mathcal {G}$ depending on $f, \boldsymbol {E}$ and $\boldsymbol {A}_\perp$, we have

(2.17)\begin{align} \{\mathcal{F}, \mathcal{G}\} & = \int f\left[\frac{\delta\mathcal{F}}{\delta f},\frac{\delta\mathcal{G}} {\delta f}\right]_{{xp}}\mathrm{d}{x}\,\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s}+\int \left(\frac{\delta\mathcal{F}}{\delta {E_x}}\frac{\partial f}{\partial {p}} \frac{\delta\mathcal{G}}{\delta f}-\frac{\delta\mathcal{G}}{\delta {E_x}} \frac{\partial f}{\partial {p}}\frac{\delta\mathcal{F}}{\delta f}\right)\mathrm{d}{x}\,\mathrm{d}{p}\,\mathrm{d}\boldsymbol{s} \nonumber\\ & \quad + \int \left(\frac{\delta \mathcal{G}}{\delta \boldsymbol{A}_\perp} \boldsymbol{\cdot} \frac{\delta \mathcal{F}}{\delta \boldsymbol{E}_\perp} - \frac{\delta \mathcal{F}}{\delta \boldsymbol{A}_\perp} \boldsymbol{\cdot} \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}_\perp}\right) \mathrm{d}x \nonumber\\ & \quad + \frac{1}{\beta\mathfrak{h}}\int f \boldsymbol{s}\boldsymbol{\cdot} \left(\frac{\partial}{\partial \boldsymbol{s}}\frac{\delta \mathcal{F}}{\delta {f}} \times \frac{\partial}{\partial \boldsymbol{s}}\frac{\delta \mathcal{G}}{\delta {f}}\right) \mathrm{d}x\,\mathrm{d}p\,\mathrm{d}\boldsymbol{s}. \end{align}

With the Hamiltonian functional defined by

(2.18)\begin{align} \mathcal{H}(f, \boldsymbol{E}, \boldsymbol{A}_\perp) & = \frac{1}{2}\int p^2 f \,\mathrm{d}x\,\mathrm{d}p\,\mathrm{d}\boldsymbol{s} + \frac{1}{2}\int |\boldsymbol{A}_\perp|^2 f \,\mathrm{d}x\,\mathrm{d}p\,\mathrm{d}\boldsymbol{s} \nonumber\\ & \quad + \frac{1}{2}\int {\left( |\boldsymbol{E}|^2 + \left|\frac{\partial \boldsymbol{A}_\perp}{\partial x}\right|^2 \right)} \mathrm{d}x + \beta \mathfrak{h}\int \left( s_2 \frac{\partial A_z}{\partial x} - s_3 \frac{\partial A_y}{\partial x} \right) f \,\mathrm{d}x\,\mathrm{d}p\,\mathrm{d}\boldsymbol{s}, \end{align}

the scalar spin Vlasov–Maxwell system of equations (2.15) can thus be reformulated as $\partial _t\mathcal {Z} = \{\mathcal {Z}, \mathcal {H} \}$, where $\mathcal {Z}=(f, E_x, E_y,E_z,A_y,A_z)$ denotes the vector of the dynamical variables.

2.2.2. The 1-D vector model

Here, the goal is to derive an equivalent model of (2.15) satisfied by $(f_0,\boldsymbol {f})(t,x,p)$. To do so, we use Proposition 1 which states that the scalar model is equivalent to the vector model whose unknown are given by some moments in the $\boldsymbol {s}$ variable. We now state the following Proposition 2, which is adapted to the 1-D laser–plasma model we consider here. Proposition 2 is a direct consequence of Proposition 1, hence the derivation of the vector model is postponed to Appendix B.

Proposition 2 Let $f(t)$ the solution to the scalar model (2.15) (with $\gamma =\beta$ and $\alpha \in \mathbb {R}$) with initial conditions (2.16a–d) with $f^0=({1}/{4{\rm \pi} }) (f_0^0 + \boldsymbol {s}\boldsymbol {\cdot } \boldsymbol {f}^0)$. Then, $(f_0(t), \boldsymbol {f}(t))$ is the solution of the following vector model where $\boldsymbol {B} = (0, -\partial _x A_z, \partial _x A_y)$:

\begin{gather} \left\{\begin{gathered} \frac{\partial f_0}{\partial t} + p \frac{\partial f_0}{\partial x} + \left(E_x - \boldsymbol{A}_\perp \boldsymbol{\cdot} \frac{\partial \boldsymbol{A}_\perp}{\partial x} \right) \frac{\partial f_0}{\partial p} - \frac{(\beta+2\gamma )}{3}\mathfrak{h} \frac{\partial^2 A_z}{\partial x^2}\frac{\partial f_2}{\partial p} + \beta\mathfrak{h} \frac{\partial^2 A_y}{\partial x^2} \frac{\partial f_3}{\partial p} = 0,\\ \frac{\partial \boldsymbol f}{\partial t} + p\frac{\partial \boldsymbol f}{\partial x} + \left(E_x - \boldsymbol{A}_\perp \boldsymbol{\cdot} \frac{\partial \boldsymbol{A}_\perp}{\partial x} \right) \frac{\partial \boldsymbol f}{\partial p}+ \boldsymbol{B} \times \boldsymbol{f} + \beta \mathfrak{h} \frac{\partial \boldsymbol{B}}{\partial x} \frac{\partial f_0}{\partial p}=0, \\ \frac{\partial E_x}{\partial t} ={-}\int_{\mathbb{R}} p f_0 \,\mathrm{d}p,\\ \frac{\partial E_y}{\partial t} ={-} \frac{\partial^2 A_y}{\partial x^2} + A_y \int_{\mathbb{R}} f_0 \,\mathrm{d}p + \frac{{ \alpha}\mathfrak{h}}{3} \int_{\mathbb{R}} \frac{\partial f_3}{\partial x}\mathrm{d}{p},\\ \frac{\partial E_z}{\partial t} ={-} \frac{\partial^2 A_z}{\partial x^2} + A_z \int_{\mathbb{R}} f_0 \,\mathrm{d}p -\frac{{ \alpha}\mathfrak{h}}{3} \int_{\mathbb{R}} \frac{\partial f_2}{\partial x}\mathrm{d}{p},\\ \frac{\partial \boldsymbol{A}_\perp}{\partial t} ={-} \boldsymbol{E}_\perp,\\ \frac{\partial E_x}{\partial x} = \int_{\mathbb{R}} f_0 \,\mathrm{d}{p} - 1, \end{gathered}\right. \end{gather}

with the initial conditions $(f_0(t=0), \boldsymbol {f}(t=0)) = (f_0^0, \boldsymbol {f}^0), {\partial E_x}/{\partial x}$ $= \int _{\mathbb {R}} f^0_0 \,\mathrm {d}{p} - 1, (E_y, E_z)(t=0)=(E_y^0, E_z^0), \boldsymbol {A}_\perp (t=0) = \boldsymbol {A}_\perp ^0$. Conversely, if $(f_0(t), \boldsymbol {f}(t))$ is the solution of (2.19), thus $f(t)=({1}/{4{\rm \pi} })(f_0(t) + \boldsymbol {s}\boldsymbol {\cdot } \boldsymbol {f}(t))$ is the solution of the scalar model (2.15) (with $\gamma =\beta, \alpha \in \mathbb {R}$).

2.2.3. Poisson structure

The model (2.19) with $\alpha =\beta +2\gamma$ enjoys a Poisson structure with the following (antisymmetric) bracket:

(2.20)\begin{align} \{\mathcal{F},\mathcal{G}\} & = \int {f_0} \left[\frac{\delta \mathcal{F}}{\delta f_0}, \frac{\delta \mathcal{G}}{\delta f_0}\right]_{{xp}}\mathrm{d}{x}\,\mathrm{d}{p} + \sum_{i=1}^3\int {f_i} \left[\frac{\delta \mathcal{F}}{\delta f_0},\frac{\delta \mathcal{G}} {\delta f_i}\right]_{{xp}}\mathrm{d}{x}\,\mathrm{d}{p} \nonumber\\ & \quad +\sum_{i=1}^3\int {f_i} \left[\frac{\delta \mathcal{F}}{\delta f_i},\frac{\delta \mathcal{G}} {\delta f_0}\right]_{{xp}}\mathrm{d}{x}\,\mathrm{d}{p} \nonumber\\ & \quad + \frac{3\beta}{\alpha} \sum_{i=1}^3 \int {f_0} \left[\frac{\delta \mathcal{F}}{\delta f_i}, \frac{\delta \mathcal{G}} {\delta f_i}\right]_{{xp}}\mathrm{d}{x}\,\mathrm{d}{p} \nonumber\\ & \quad + \sum_{i=0}^3 \int \left(\frac{\delta \mathcal{F}}{\delta {E_x}}\frac{\partial f_i}{\partial {p}} \frac{\delta \mathcal{G}}{\delta f_i}-\frac{\delta \mathcal{G}}{\delta {E_x}} \frac{\partial f_i}{\partial {p}}\frac{\delta \mathcal{F}}{\delta f_i}\right)\mathrm{d}{x}\,\mathrm{d}{p} \nonumber\\ & \quad +\frac{3}{\alpha \mathfrak{h}}\int \boldsymbol{f} \boldsymbol{\cdot} \left(\frac{\delta \mathcal{F}}{\delta \boldsymbol{f}} \times \frac{\delta \mathcal{G}}{\delta \boldsymbol{f}}\right)\mathrm{d}{x}\,\mathrm{d}{p}+ \int \left( \frac{\delta \mathcal{G}}{\delta \boldsymbol{A}_\perp} \boldsymbol{\cdot} \frac{\delta \mathcal{F}}{\delta \boldsymbol{E}_\perp} - \frac{\delta \mathcal{F}}{\delta \boldsymbol{A}_\perp} \boldsymbol{\cdot} \frac{\delta \mathcal{G}}{\delta \boldsymbol{E}_\perp}\right) \mathrm{d}x, \end{align}

where $\boldsymbol {f} = (f_1, f_2,f_3)$. The Hamiltonian functional is given by

(2.21)\begin{align} \mathcal{H}(f_0, \boldsymbol{f}, \boldsymbol{E}, \boldsymbol{A}_\perp) & = \frac{1}{2}\int p^2 f_0 \,\mathrm{d}x\,\mathrm{d}p+ \frac{1}{2} \int |\boldsymbol{A}_\perp|^2 f_0 \,\mathrm{d}x\,\mathrm{d}p \nonumber\\ & \quad + \frac{1}{2}\int \left(|\boldsymbol{E}|^2 + \left|\frac{\partial \boldsymbol{A}_\perp}{\partial x}\right|^2\right) \mathrm{d}x + \frac{\alpha \mathfrak{h}}{3}\int \left( f_2 \frac{\partial A_z}{\partial x} - f_3 \frac{\partial A_y}{\partial x} \right) \mathrm{d}x\,\mathrm{d}p, \end{align}

so that the system can be written in a compact form as

(2.22)\begin{equation} \frac{\partial \mathcal{Z}}{\partial t} = \{ \mathcal{Z}, \mathcal{H} \}, \end{equation}

where $\mathcal {Z}=(f_0, \boldsymbol {f}, E_x, E_y,E_z, A_y,A_z)$ denotes the vector of the dynamical variables.

3.2.1. Full spin-dependent models at short times

In this part, we focus on the solutions for short times (approximately $100 \omega _p^{-1}$), which enables us to employ refined numerical parameters in order to compare as accurately as possible the two models.

From Proposition 1, we proved the equivalence between the scalar and vector models. However, this equivalence is only true for $\gamma \neq 0$, so that the term $\gamma \mathfrak {h}\boldsymbol {\nabla } (\boldsymbol {B}\boldsymbol {\cdot }{\partial }/{\partial \boldsymbol {s}})\boldsymbol {\cdot } {\partial f}/{\partial \boldsymbol {p}}$ is present in the scalar Vlasov equation (2.15). This term cannot be easily described using standard PIC methods based on trajectories, hence, it is neglected most of the time. However, when $\mathfrak {h}=0$ this term disappears, so that the two models are again equivalent for any value of $\gamma$. Then, by progressively increasing $\mathfrak {h}$, we can show how the results of the scalar and vector models depart from each other. The numerical solutions are obtained using the PIC code described in our earlier work (Crouseilles et al. Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021) and the grid-based Eulerian code presented in the preceding section of the present paper.

To control the influence of the numerical error, we use a refined mesh for both methods: $N_x=258$ and $N_p=10^3$ for the grid-based method, and $N_x=258$ with number of particles $N_{\textrm {part}}=10^5$ for the PIC method. The initial conditions are given by (3.10) with $\eta =0.5, \alpha =0.02, v_{\textrm {th}}=0.17$ and $E_0=E_{\textrm {ref}}$. Finally, the time step is $\Delta t =0.004$ for both cases (both methods use a Strang splitting) and we consider the following values for the scaled Planck constant $\mathfrak {h}=0,\ 0.001,\ 0.1,\ 0.25,\ 0.5,\ 0.75$ and $1$. We note that, strictly speaking, for such large values of $\mathfrak {h}$ a relativistic treatment of the motion should be adopted. However, here the purpose is simply to study the agreement between the scalar and vector models, so that all considerations about relativistic corrections are ignored.

In figure 3, we study the difference (‘error’) between the scalar and vector models by considering the perpendicular electric energy $\mathcal {H}_{E,\perp }=\frac {1}{2} \int (E_y^2+E_z^2)\,\mathrm {d}x$ and the magnetic energy $\mathcal {H}_B=\frac {1}{2} \int (B_y^2+B_z^2)\,\mathrm {d}x$. We plot, on a log–log scale, the error (in the $L^\infty$ norm) of these energies obtained from the scalar and vector models, for different values of $\mathfrak {h}$. For both energies, the difference decreases with decreasing $\mathfrak {h}$ in an approximately linear fashion. This is in agreement with the fact that, for $\gamma =0$, the difference between two models is expected to be $O(\mathfrak {h})$. For very small values of ${\mathfrak {h}}$, the difference between the two models saturates to the numerical error due to the different numerical methods (PIC and Eulerian). The horizontal dashed line in figure 3 corresponds to such numerical error observed for $\mathfrak {h}=0$.

Figure 3. Difference ($L^\infty$ norm) between the scalar (superscript $S$) and vector (superscript $V$) models for the perpendicular electric energy ($\mathcal {H}_{E,\perp }$, a) and magnetic energy ($\mathcal {H}_B$, b), for different values of the scaled Planck constant $\mathfrak {h}$. The horizontal dashed lines correspond to the numerical error observed for $\mathfrak {h}=0$. The green straight line has a slope equal to unity.

3.2.2. Spin-dependent models without wave self-consistency at long times

Here, we study the propagation of a circularly polarized wave into a plasma that does not retroact on the wave in the transverse direction, although the self-consistency along the propagation direction is maintained through the Poisson equation (Crouseilles et al. Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021). To simulate this, we consider a reduced model. Using the above notations, we remove the term $\frac {1}{2}\int _{\varOmega }|\boldsymbol {A}_\perp |^2 f_0\,\mathrm {d}x\,\mathrm {d}p$ from the Hamiltonian $\mathcal {H}$ in (2.21), then through the Poisson bracket representation (2.22) we can derive a reduce system similar to (2.19), but without the terms $\boldsymbol {A}_\perp \boldsymbol {\cdot } ({\partial \boldsymbol {A}_\perp }/{\partial x}) ({\partial f_0}/{\partial p})$ and $\boldsymbol {A}_\perp \boldsymbol {\cdot } ({\partial \boldsymbol {A}_\perp }/{\partial x}) ({\partial \boldsymbol {f}}/{\partial p})$ in the Vlasov equations, and the terms $A_y \int _{\mathbb {R}} f_0 \,\mathrm {d}{p}$ and $A_z \int _{\mathbb {R}} f_0 \,\mathrm {d}{p}$ in the sources of the Maxwell equations. In this case, the electromagnetic wave is not coupled to the plasma and its evolution can be determined exactly by solving the corresponding Maxwell equations for $E_y$ and $E_z$. In contrast, the longitudinal nonlinearity is kept in the model, hence, $E_x$ is a solution of Poisson's equation.

We consider the initial condition (3.10) for the distribution functions and we use ${\eta =0.5,\ \alpha =0.02, v_{\textrm {th}}=0.17}$ and $\mathfrak {h}=0.00022$. The electromagnetic fields are initialized as in the spinless case (see § 3.1), but here we have different matching conditions for the circularly polarized wave, since the wave propagates in vacuum in this case. We still have the relation (3.3a,b) together with the following ones:

(3.11a,b)\begin{equation} \omega_{0,s}=k_{0,s},\quad \omega_e^2=1+3k_e^2v_{{\rm th}}^2. \end{equation}

We still consider $k_0=2k_e$, so that we obtain from the matching relations (3.3a,b) and (3.11a,b),

(3.12a,b)\begin{equation} \omega_0=2.0928,\quad k_s=k_e=\omega_s=\omega_e=1.0464. \end{equation}

The amplitude of the incident wave is $E_0=E_{\textrm {ref}}=0.325$.

Neglecting the spatial dependence and considering the following time dependent magnetic field $\boldsymbol {B}=(0,B_0\sin (\omega _0 t),-B_0\cos (\omega _0 t))$ with $B_0=E_0k_0/\omega _0$, one can obtain an approximate closed equation for the dynamics of the macroscopic spin ${S}_z(t)$ (with ${S}_z(t)=\frac {1}{3}\int f_3 \,\mathrm {d}x\,\mathrm {d}p$ for the vector case and ${S}_z(t)=\int s_3 f \,\mathrm {d}x\,\mathrm {d}p\, \mathrm {d}\boldsymbol {s}$ for the scalar case). In the regime $B_0/\omega _0 \ll 1$ (i.e. when the Larmor precession frequency $eB_0/m$ is much smaller than the laser frequency), it can be shown that the spin $S_z(t)$ oscillates with a frequency $\omega _\textrm {spin}=B_0^2/(2\omega _0)$ (see Crouseilles et al. (Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021) for more details).

On the numerical side, we use the same parameters as in § 3.1 (i.e. $N_x=N_p=129$ and $\Delta t=0.04$) for the Eulerian method and $N_x=128,\ N_{\textrm {part}}=2\times 10^4,\ \Delta t=0.04$ for the PIC method. Different values of the incident wave amplitude $E_0$ are studied to compare the numerical results with the analytic frequency $\omega _{\textrm {spin}}=B_0^2/(2\omega _0)$. Indeed, as the transverse retroaction terms were removed from the Vlasov equation, we do not expect an instability here, so that the amplitude of the incident wave should remain approximately constant, which justifies the analytical calculations.

In figure 4, the time evolution of $S_z$ obtained by the vector model (Eulerian method) and scalar model (PIC method) are displayed for three values of the amplitude of the incident electromagnetic wave $E_0$. The corresponding spectrum, normalized to its peak value, is also shown in figure 4(b,d, f). We observe that $S_z$ displays an oscillatory damped behaviour for all cases, with a damping rate that is significantly higher when $E_0$ increases. The damping rate is slightly larger for the vector model at long times, but the agreement between the two approaches is very good for short times (up to ${\approx }1500 \omega _p^{-1}$). In the spectrum, the expected analytical frequencies $\omega _{\textrm {spin}}= 0.0039,\ 0.0158$ and $0.063$ are also recovered with good precision. In particular, the quadratic scaling between $\omega _{\textrm {spin}}$ and $B_0$ is correctly reproduced.

Figure 4. Spin-dependent reduced model without wave self-consistency. Time evolution of $S_z$ (a,c,e) and frequency spectrum of $S_z$ (b,d, f), for the scalar (red lines) and vector (blue lines) models. We recall that $S_z=\int s_3 f \,\mathrm {d}x\,\mathrm {d}p\,\mathrm {d}\boldsymbol {s}$ for the vector model (Eulerian code) and $S_z=\frac {1}{3}\int f_3 \,\mathrm {d}x\,\mathrm {d}p$ for the scalar model (PIC code). The amplitude of the incident wave is, from top to bottom, $E_0=0.5E_{\textrm {ref}}$, $E_0=E_{\textrm {ref}}$ and $E_0=2E_{\textrm {ref}}$. The analytical predictions for the peak frequency are $\omega _{\textrm {spin}}= 0.0039,\ 0.0158$ and $0.063$, from top to bottom.

3.2.3. Full spin-dependent models at long times

Here, we present numerical results corresponding to $\alpha =1,\ \beta =1, \gamma =0$ for the scalar and vector models for large times (the short time behaviour was investigated in § 3.2.1). In this case, the models are not equivalent as soon as $\mathfrak {h}\neq 0$; however, the Poisson bracket structure is ensured for both models, so that the Hamiltonian splitting can be used (Strang splitting) to get a good total energy conservation for both codes. The initial condition for the distribution functions is that of (3.10) and the electromagnetic waves are circularly polarized, with the same parameters as in the spinless case of § 3.1.

Eulerian code results (vector model). First, we illustrate the general results obtained with the Eulerian (vector approach), before proceeding to the systematic comparison with the PIC code (scalar approach). The results, which were obtained with $\mathfrak {h}=0.00022$, are shown in figures 5 and 6.

Figure 5. Full spin-dependent vector model (Eulerian method with $N_x=259, N_p=512,\ p\in [-2.5,2.5],\ \Delta t=0.02$). (a) Longitudinal electric field norm $\|E_x \|$ in semilog scale (the inset shows a zoom on short times to highlight the linear instability phase with slope ${\approx }0.03$). (b) Magnetic energy $\frac {1}{2}\int |\boldsymbol {B}|^2\,\mathrm {d}x$. (c) Spin component ${S}_z(t)$.

Figure 6. Full spin-dependent vector model (Eulerian method with $N_x=129, N_p=256,\ p\in [-2.5,2.5],\ \Delta t=0.04$). From (a) to (d): contour plots of the four components of the electron distribution function $f_0$, $f_1$, $f_2$, $f_3$ in the phase space $(x,p)$, at time $t=320 \omega _p^{-1}$.

The instability rate measured from the growth of the longitudinal electric field amplitude $\|E_x(t)\|$ is close to the one observed in the spinless simulations ($\gamma _{\textrm {inst}} \approx 0.03$). Moreover, the magnetic field energy strongly decreases during the linear phase (where the longitudinal electric energy increases exponentially) to reach a plateau, before another strong drop around $t\approx 2000 \omega _p^{-1}$. Finally, the spin component ${S}_z(t)= \frac {1}{3}\int f_3 \,\mathrm {d}x\,\mathrm {d}p$ has a damped oscillatory behaviour which is qualitatively similar to the results obtained with the PIC approach (Crouseilles et al. Reference Crouseilles, Hervieux, Li, Manfredi and Sun2021).

We also show, in figure 6, the phase space portraits of the distribution functions $f_0, f_1, f_2, f_3$ at time $t=320\omega _p^{-1}$. Two phase-space vortices can be clearly seen in these figures thanks to the high-order numerical methods used for these Eulerian simulations. As expected, the size of each vortex is close to the wavelength of the unstable plasma wave $2{\rm \pi} /k_e \approx 5.15$ (see § 3.1). We note that the statistical noise inherent to PIC codes would have precluded such high-resolution diagnostics.

Comparison with PIC code (scalar model) for $\mathfrak {h}=0$. In figure 7, we compare the vector and scalar models with $E_0=2E_{\textrm {ref}}$ and $\mathfrak {h}=0$. When $\mathfrak {h}=0$ and $\alpha =\beta =1,\ \gamma =0$, the scalar and vector models are mathematically equivalent, but the results may differ because of the different numerical methods that are used. Therefore, this is a good comparative test of the two numerical methods.

Figure 7. Comparison of the full spin-dependent vector and scalar models with $\mathfrak {h}=0.0$ and $E_0=2 E_{\textrm {ref}}$ (Eulerian method with $p\in [-5,5]$ $N_x=259,\ N_p= 512,\ \Delta t=0.02$ and PIC method with $N_x=259,\ N_{\textrm {part}}= 10^5,\ \Delta t=0.02$). (a,c,e,g) Vector model with Eulerian method; (b,d, f,h) scalar model with PIC method. From top to bottom are shown (a,b) the longitudinal electric field, (c,d) the magnetic energy, (e, f) the spin component $S_z$, and (g,h) its frequency spectrum.

First, we observe that the instability rate on $\|E_x(t)\|$ is similar for the two approaches, and very close to the one observed in the spinless case ($\gamma _\textrm {inst}\approx 0.06$). We also checked that the growth rate is linear with respect to $E_0$ by running the case $E_0=E_{\textrm {ref}}$ (twice as small as in figure 7) and that the total energy is very well conserved (with an error ${\approx }10^{-5}$ even for long times, for both codes). Second, we can see that the spin frequencies obtained from the PIC and Eulerian simulations are very similar (${\approx }0.08\omega _p$).

Generally speaking, the results of the two codes are in good agreement, particularly during the initial linear phase. The remaining observed differences should be attributed to the two numerical methods (PIC and Eulerian), which are intrinsically different in their approaches.

Effect of the scaled Planck constant. In figures 8, 9 and 10, the influence of $\mathfrak {h}$ ($\mathfrak {h}=0.01, 0.1$ and $0.5$) is studied for the vector and scalar models for $E_0=2E_{\textrm {ref}}$. The same initial condition as before is used and $\alpha =\beta =1,\ \gamma =0$. Let us recall that here, since $\mathfrak {h}\neq 0$, the two models are not mathematically equivalent. The linear phase becomes quite different as $\mathfrak {h}$ increases: in particular, the electric energies do not saturate at the same level. However, the evolution of ${S}_z$ is similar for both approaches: the oscillations are damped more quickly as $\mathfrak {h}$ increases. For both codes, the total energy is well preserved (the relative error is approximately $10^{-4}$, not shown). Finally, the magnetic energies, even though they both decrease in the two approaches, have rather different behaviours, with a bigger drop observed for the vector model. Of course, the differences can be explained by the fact that in this case the underlying models are not equivalent. We emphasize that the Eulerian code has been tested on several meshes to ensure the numerical convergence, whereas the PIC code would require, in order to achieve convergence, a large number of particles $N_p$ that would lead to very time-consuming simulations.

Figure 8. Comparison of the full spin-dependent vector and scalar models for $\mathfrak {h}=0.01$ and $E_0=2E_{\textrm {ref}}$ (Eulerian method with $p\in [-5,5]$ $N_x=259,\ N_p= 512,\ \Delta t=0.02$ and PIC method with $N_x=259,\ N_{\textrm {part}}= 10^5,\ \Delta t=0.02$). (a,c,e) Vector model with Eulerian method; (b,d, f) scalar model with PIC method. From top to bottom, we show (a,b) the time evolution of the longitudinal electric field, (c,d) the magnetic energy, and (e, f) the spin component $S_z$.

Figure 9. Same as figure 8 for $\mathfrak {h}=0.1$.

Figure 10. Same as figure 8 for $\mathfrak {h}=0.5$.

Finally, we note, by comparing figures 8(e, f), 9(e,f) and 10(e, f), that the spin oscillations are more strongly damped for larger values of $\mathfrak {h}$.

Temperature effects. In figure 11, we study the influence of a higher electron temperature by taking $v_{\textrm {th}}=0.51$ in the initial condition (3.10) for the Eulerian and PIC codes (instead of $v_{\textrm {th}}=0.17$ in figure 5). The field amplitude is fixed to $E_0=E_{\textrm {ref}}$ and $\mathfrak {h}=0.00022$. For the present case $v_{\textrm {th}}=0.51$, the matching relations (3.3a,b) and (3.4a,b) yield the wavenumbers $k_e=1.46$ and $k_0=2k_e$. In figure 11, we show the time evolution of the electric energy and the spin component ${S}_z$ for the Eulerian and PIC methods. These results have to be compared with the case $v_{\textrm {th}}=0.17$ presented in figure 5.

Figure 11. Influence of the electron temperature using $v_{\textrm {th}}=0.51$. (a,b) Time evolution of the electric energy (semilog scale). (c,d) Time evolution of the spin $S_z$. (a,c) Full spin-dependent vector model (Eulerian method with $N_x=N_p=129,\ \Delta t=0.04$). (b,d) Full spin-dependent scalar model (PIC method with $N_x=129,\ N_{\textrm {part}}=2\times 10^4,\ \Delta t=0.04$).

Thanks to its superior resolution, the Eulerian code allows the initial longitudinal electric field to be very small (figure 11a), in accordance with the initial condition (3.10). Subsequently, this initial perturbation grows unstable and saturates nonlinerarly at a certain (still rather small) value. The spikes observed on the figure correspond to the well-known recurrences occurring at times multiple of $2{\rm \pi} /(k_e \Delta p)$, where $\Delta p$ is the grid step in momentum space. These recurrences are due to the discretization of momentum space. In contrast, the PIC code (figure 11b) cannot resolve the initial small perturbation because of its inherent noise, so the instability is lost. Indeed, the noise level of the PIC code is even higher than the saturation level observed for the Eulerian code.

In spite of these limitations, the spin component $S_z$ decreases in a similar fashion for both methods, and becomes almost zero after $\omega _p t=1000$.

By comparing the results of figure 11 ($v_{\textrm {th}}=0.17$) and figure 5 ($v_{\textrm {th}}=0.51$), we also note that the spin component $S_z$ is more strongly damped for the case at higher temperature, as should be expected.

Full vector model with $\gamma \neq 0$. In this last part, we simulate the full spin-Vlasov–Maxwell vector model introduced in Hurst et al. (Reference Hurst, Morandi, Manfredi and Hervieux2014, Reference Hurst, Morandi, Manfredi and Hervieux2017) and Manfredi et al. (Reference Manfredi, Hervieux and Crouseilles2022), i.e. for $\gamma =1,\alpha =3,\beta =1$, for which we cannot perform simulations for the scalar model, since the cross-derivative term is not implemented in our PIC code. In particular, we compare the results of the Eulerian code for this set of parameters (with $\gamma \neq 0$) with those obtained with the same code but $\gamma = 0$ and $\alpha =\beta =1$ (see figure 9).

In figure 12, we show the comparison for a typical case with $E_0=2E_{\textrm {ref}}$ and $\mathfrak {h}=0.1$. The results are in reasonable agreement: in particular, the longitudinal electric fields and the magnetic energies saturate nonlinearly at similar levels. The $z$ component of the spin appears to decay faster and with fewer oscillations in the $\gamma \neq 0$ case. For both cases, the total energy is very well preserved (around $10^{-4}$, not shown).

Figure 12. Eulerian method for the vector model with $\alpha =\beta =1, \gamma =0$ (a,c,e) and ${\alpha =3,\beta =\gamma =1}$ (b,d, f) with $\mathfrak {h}=0.1$. $N_x=259, N_p=10^3,\ \Delta t=0.02$, $E_0=2E_{\textrm {ref}}$. From top to bottom, we show (a,b) the norm of the longitudinal electric field, (c,d) the magnetic energy, and (e, f) the $z$ component of the spin.