6D phase space collective modes in Vlasov-Maxwell system

According to MEs, scalar potential $\phi$ and vector potential $\vec{A}$ satisfy wave equations with different 'source' terms: $n\equiv \displaystyle \int f{d}^{3}\upsilon$ for $\phi$ and $\vec{j}\equiv \displaystyle \int \upsilon f{d}^{3}\upsilon$ for $\vec{A},$ where $f$ is the PDF. Lorentz gauge condition ${\partial }_{t}\phi +{\rm{\nabla }}\cdot \vec{A}=0$ relates $\phi$ and $\vec{A}$ and agrees with continuity equation ${\partial }_{t}n+{\rm{\nabla }}\cdot \vec{j}=0$ reflecting conservation of total particles number. According to standard partial differential equation (PDE) theory, $\phi$ and $\vec{A}$ solutions can be expressed in a general form [51]:

$$\begin{eqnarray}&&F\left(\vec{r},t\right)=\iiint \frac{{FS}\left({\vec{r}}^{{\prime} },{t}^{{\prime} }\right)}{{\rm{| | }}\vec{r}-{\vec{r}}^{{\prime} }{\rm{| }}-{ct}+c{t}^{{\prime} }{\rm{| }}}{d}^{3}{\vec{r}}^{{\prime} }d{t}^{{\prime} }+{\partial }_{t}\left[t\cdot \iint F{{\rm{| }}}_{t=0}{d}^{2}S\right]+t\cdot \iint {\partial }_{t}F{{\rm{| }}}_{t=0}{d}^{2}S,\end{eqnarray} \tag{ 1 }$$

where $F$ stands for $\phi$ or $\vec{A},$ ${FS}$ stands for its source term charge density $n$ or current density $\vec{j},$ and ${dS}$ is a sphere surface which is centered at $r$ and of a radius ${ct}.$ This general form indicates that $F$ depends on its two initial conditions: $F{{\rm{| }}}_{t=0}$ and ${\partial }_{t}F{{\rm{| }}}_{t=0}.$ Due to Lorentz gauge condition, ${\partial }_{t}\phi {{\rm{| }}}_{t=0}$ can be simultaneous with $\vec{A}{{\rm{| }}}_{t=0},$ and further with $f{{\rm{| }}}_{t=0}.$ But ${\partial }_{t}\vec{A}{{\rm{| }}}_{t=0},$ or initial transverse electric field profile ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0},$ has a different simultaneous partner. Now that $\vec{A}$ satisfies a wave equation with a source $\vec{j},$ ${\partial }_{t}\vec{A}$ will satisfy a wave equation with a source ${\partial }_{t}\vec{j}.$ Thus, ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0}\equiv -{\partial }_{t}\vec{A}{{\rm{| }}}_{t=0}$ is simultaneous with ${\partial }_{t}f{{\rm{| }}}_{t=0}.$ Moreover, when ${FS}$ in equation (1) obeys an equation dependent on ${\partial }_{t}F{{\rm{| }}}_{t=0},$ such as Vlasov equation (VE) containing ${E}_{{Tr}}\equiv {\partial }_{t}A,$ by setting the result of the operator ${\partial }_{t}$ on the formula equation (1) of $F\equiv A$ at $t=0,$ we can easily obtain a complicated self-consistent equation of ${\partial }_{t}A{{\rm{| }}}_{t=0}$ (because ${\partial }_{t}{FS}$ depends on ${\partial }_{t}A{{\rm{| }}}_{t=0}$). Namely, it is not an easy task to solve an initial condition ${\partial }_{t}A{{\rm{| }}}_{t=0}$ consistent with others such as $f{{\rm{| }}}_{t=0},$ $A{{\rm{| }}}_{t=0}$ etc

$\phi$ and $\vec{A}$ appearing in the VE of $f,$ which formally is a PDE defined over a $6D$ phase space plus a time coordinate, can be completely expressed, through above formula, in terms of $f.$ This makes the VE finally to become a differential-integral equation of $f,$ or an integral equation of local time growth-rate ${\partial }_{t}f$ because ${\vec{E}}_{{Tr}}$ is simultaneous with ${\partial }_{t}f:$

$$\begin{eqnarray}&&{\partial }_{t}f+e{\vec{E}}_{{Tr}}\left[\int \upsilon {\partial }_{t}f{d}^{3}\upsilon \right]\cdot {\partial }_{p}f+O\left[f\right]=0,\end{eqnarray} \tag{ 2 }$$

where $e$ is electronic charge, and ${\vec{E}}_{{Tr}}\left[\int \upsilon {\partial }_{t}f{d}^{3}\upsilon \right]$ represents ${\vec{E}}_{{Tr}}$ being a functional of $\int \upsilon {\partial }_{t}f{d}^{3}\upsilon .$ Moreover, all other terms in the VE are represented by a functional of $f,$ which is denoted as $O\left[f\right]$ and $O$ is the initial letter of ${Others}.$ Its detailed form is easily written out but somewhat long for being presented here. Here, magnetic field $\vec{B}$ and longitudinal electric field ${\vec{E}}_{{Lo}}$ are simultaneous with $f$ and hence contained in ${Others}$-functional.

Once specifying the value of ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0}\equiv -{\partial }_{t}\vec{A}{{\rm{| }}}_{t=0},$ we have implicitly specified the value of ${\partial }_{t}f{{\rm{| }}}_{t=0}$ to some extent. Thus, when conducting the V-M simulation in Eulerian approach, substituting this ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0},$ (i.e., the specified ${\partial }_{t}f{{\rm{| }}}_{t=0},$) into the VE which is treated as an 'algebra' equation of ${\partial }_{t}f{{\rm{| }}}_{t=0},$ it is necessary to examine whether the implicitly-specified ${\partial }_{t}f{{\rm{| }}}_{t=0}$ agrees with the calculated ${\partial }_{t}f{{\rm{| }}}_{t=0}.$ For brevity, we denote them as ${\partial }_{t}f{{\rm{| }}}_{t=0}^{{MEs}}$ and ${\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}}$ (where upper index 'MEs' means it being determined by MEs and 'VE' by VE), respectively. If they are different or ${\rm{| }}{\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}}-{\partial }_{t}f{{\rm{| }}}_{t=0}^{{MEs}}{\rm{| }}$ is not small enough, using the ${\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}}$ to calculate $f{{\rm{| }}}_{t={\rm{\Delta }}t}$ is futile. (In other words, because we can construct, according to the simultaneity relation above-mentioned, a ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0}^{{VE}}$ from ${\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}},$ we should examine whether ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0}^{{VE}}$ agreeing with the specified value of ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0},$ or ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0}^{{MEs}}.$) Once disagreement occur, we have to try a new trial of ${\vec{E}}_{{Tr}}{{\rm{| }}}_{t=0},$ as well as a new trial of ${\partial }_{t}f{{\rm{| }}}_{t=0}^{{MEs}}$-profile, and calculate its yielded ${\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}}$-profile, and then compare them. Repeating such a cycle until ${\rm{| }}{\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}}-{\partial }_{t}f{{\rm{| }}}_{t=0}^{{MEs}}{\rm{| }}$ being small enough (if can). For example, $\frac{{\rm{| }}{\partial }_{t}f{{\rm{| }}}_{t=0}^{{VE}}-{\partial }_{t}f{{\rm{| }}}_{t=0}^{{MEs}}{\rm{| }}}{{\rm{| }}{\partial }_{t}f{{\rm{| }}}_{t=0}^{{MEs}}{\rm{| }}}\lt 10 \%$ is a tolerable level of the relative error. Clearly, such iterative manner is very time-consuming in computer experiment but strictly speaking unable to be omitted.

Of course, if not specifying the ${\partial }_{t}f{{\rm{| }}}_{t=0}$ but taking the VE as an integral equation of ${\partial }_{t}f,$ discrete expression of the integral $\int \upsilon {\partial }_{t}f{d}^{3}\upsilon$ in the VE will lead to a high-dimension matrix, or those values of ${\partial }_{t}f$ at different $\left(t,r,\upsilon \right)$-coordinates will form, (because discretization of the integral with respect to $\upsilon$ needs $\upsilon$-points in sufficient amount), a ${\rm{\infty }}\times {\rm{\infty }}$ matrix equation nearly impossible to-be-solved. Moreover, singularities associated with $\frac{1}{{\rm{| | }}\vec{r}-{\vec{r}}^{{\prime} }{\rm{| }}-{ct}+c{t}^{{\prime} }{\rm{| }}}$ and $\frac{1}{{\rm{| | }}\vec{r}-{\vec{r}}^{{\prime} }{\rm{| }}-{ct}+c{t}^{{\prime} }{{\rm{| }}}^{2}}$ also make the calculation more complicated.

The V-M simulation in Lagrangian approach [43–50, 52–54] is to use displacements of $f$-elements in $6D$ phase space to describe time-evolution of $f.$ Here, the phrase 'displacement' represents a $1D$ 'coordinate' along a phase space contour defined by $f\left({r}_{s}-{\upsilon }_{s}{\rm{\Delta }}t,{\upsilon }_{s}-a\left({r}_{s},{\upsilon }_{s}\right){\rm{\Delta }}t\right)=f\left({r}_{s},{\upsilon }_{s}\right),$ where $a\left({r}_{s},{\upsilon }_{s}\right)$ stands for accelerating rate and hence depends on $\left({E}_{{Lo}},{E}_{{Tr}},B\right){{\rm{| }}}_{\left(r,t\right)=\left({r}_{s},0\right)}$ and relativistic factor ${\rm{\Gamma }}{{\rm{| }}}_{\upsilon ={\upsilon }_{s}}.$ These displacements also need to be examined whether their yielded ${E}_{{Tr}}$ (through MEs) agrees with ${E}_{{Tr}}$ in the accelerating rate $a\left({r}_{s},{\upsilon }_{s}\right)=a\left[\left({E}_{{Lo}},{E}_{{Tr}},B\right){{\rm{| }}}_{\left(r,t\right)=\left({r}_{s},0\right)},{\rm{\Gamma }}{{\rm{| }}}_{\upsilon ={\upsilon }_{s}}\right].$ Here $\left({E}_{{Lo}},B\right){{\rm{| }}}_{\left(r,t\right)=\left({r}_{s},0\right)}$ are functions of $\int f\left({r}_{s},{\upsilon }_{s}\right){d}^{3}{\upsilon }_{s}$ but ${E}_{{Tr}}{{\rm{| }}}_{\left(r,t\right)=\left({r}_{s},0\right)}$ is a function of $\int {\partial }_{t}f\left({r}_{s},{\upsilon }_{s}\right){d}^{3}{\upsilon }_{s}.$

Similar difficulty also exists for Newton-Maxwell system. In the PIC code [1–15], ${\vec{E}}_{{Tr}}$ on any macroparticle is of a simultaneity relation with 'accelerating rate' of other macroparticles within adjacent cells (because two interpolation functions, one for from particle to grid and the other for from grid to particle, are taken as only dependent on space coordinates). If using this ${\vec{E}}_{{Tr}}$ to calculate the 'accelerating rate' of a macroparticle, we will have a similar equation (where ${m}_{e}$ is static electronic mass and a macroparticle shares a same charge/mass ratio with an electron)

$$\begin{eqnarray}&&{m}_{e}{d}_{t}\frac{{d}_{t}{r}_{i}}{\sqrt{1-{\rm{| }}{d}_{t}{r}_{i}{{\rm{| }}}^{2}}}=e{\vec{E}}_{{Tr}}\left[{d}_{{tt}}{r}_{1\leqslant j\ne {\rm{i}}\leqslant N},{r}_{1\leqslant j\leqslant N}\right]+O\left[{r}_{1\leqslant j\leqslant N},{d}_{t}{r}_{1\leqslant j\leqslant N}\right]\end{eqnarray} \tag{ 3 }$$

and hence need to diagonalize a matrix whose dimension is the number of macroparticles within adjacent cells.

This consistency-examination of ${\vec{E}}_{{Tr}}$ is absent in various algorithms. For example, in the PIC code [1–16], although particles' information and fields' information are alternatively updated, such a consistency-examination should be made per two time-steps. Unfortunately, it is omitted. This omission might be an important factor causing various 'violations of physics laws' [17]. Moreover, it will also affect implicit EM PIC scheme which utilizes particle enslavement technique [55] because iteratively solving the root of enslaved residual $G\left({E}^{k}\right),$ where $k$ is iteration number, will encounter a difficulty arising from uncertain initial form of $G\left({E}^{0}\right)$ due to ${\vec{E}}_{{Tr}}.$

Because ${E}_{{tr}}$ is taken as playing less role in these issues [4–39], it is often not taken into account and hence EM kinetic simulations in these works are actually approximated [4–39]. For example, in some issues on plasmas in magnetic fusion device involving in kinetic simulation on ionic PDF, the self-consistent magnetic field generated by plasma currents, (which is mainly determined by electronic PDF), is neglected and hence ${E}_{{tr}}$ is indeed not taken into account [23, 28–30, 33, 34]. Some kinetic simulations are based on Vlasov-Poisson (VP) code to study electronic PDF [39–42] and ionic PDF [35, 38]. Some simulations are on 1D model and hence not take into account $B$ and ${E}_{{tr}}$ [38]. Some kinetic simulations on electronic PDF through Eulerian Vlasov code ignored ${E}_{{tr}}$ [18, 21], so did Semi-Lagrangian (SL) Vlasov code on discharge plasmas [31, 32] and PIC code for hot electron [22]. In some $3D$ EM PIC simulations on kinetic-scale plasma, even though fluctuating magnetic field is taken into account, ${E}_{{tr}}$ is still not strictly taken into account [15, 16].

Even though many authors have been aware of some inherent drawbacks of current popular simulation schemes, such as hiring truncation approximation [20], statistical noise [35, 38] etc, the most urgent task, i.e. strictly taking ${E}_{{Tr}}$ into account, have not win sufficient attention. This makes these simulation schemes under the risk of violating MEs [17] and hence discounts their effort in improving these simulation schemes [19, 23, 25]. For example, usually the kinetic simulation is taken as an initial problem evolution of an initial PDF $f{{\rm{| }}}_{t=0},$ but people often choose $f{{\rm{| }}}_{t=0}$ as having a non-zero vortex of the current density or ${\rm{\nabla }}\times \vec{j}{{\rm{| }}}_{t=0}\ne 0.$ For example, $f{{\rm{| }}}_{t=0}={\exp }\left(-\frac{{x}^{2}+{y}^{2}}{{R}^{2}}\right){\exp }\left(-\frac{1-0.01{\upsilon }_{z}/c}{\sqrt{1-{\rm{| }}\upsilon {{\rm{| }}}^{2}/{c}^{2}}}\right)$ means the initial PDF having a Gaussian space transverse shape and is a typical example of ${\rm{\nabla }}\times \vec{j}{{\rm{| }}}_{t=0}\ne 0,$ which will inevitably, according to MEs, trigger ${E}_{{Tr}}$ to be present. Namely, the choice of the initial PDF is inconsistent with hiring the ES approximation or using the VP code.

The V-M simulation in Eulerian approach is finally inevitably frustrated by the simultaneity relation between ${\vec{E}}_{{Tr}}$ and ${\partial }_{t}f.$ Ignoring this simultaneity relation can make the simulation going-on and speed-up but the reliability of yielded results is limited. Therefore, even though this Eulerian approach is applicable for electrostatic (ES) kinetic simulation [30–39, 41, 42], its suitability to EM kinetic simulation is still uncertain. This simultaneity relation also have profound effect on the V-M simulation in Lagrangian approach.

For illustrating the effect of the simultaneity relation described above on simulation, we display results based on strict treatment on the simultaneity relation and those based on approximate treatment in figures 1, 2. Here, the phrase 'strict treatment' refers to BOGRs being linked by a matrix, and the phrase 'approximate treatment' refers to simplifying this matrix as a diagonal one by approximating its non-diagonal elements as being of a different simultaneity relation (alike that between ${\vec{E}}_{{Lo}}$ and $f$) and hence merging them into the functional $O$ in equations (2), (3). Figure 1 is for the V-M simulation in Eulerian approach and figure 2 for the PIC simulation. In figure 1, $f$ starts from a familiar Gaussian initial profile $f{{\rm{| }}}_{t=0}={f}_{{\max }}^{0}{\exp }\left(-\frac{{x}^{2}+{y}^{2}}{{R}^{2}}\right){\exp }\left(-1/\sqrt{1-\frac{{\upsilon }_{x}^{2}+{\upsilon }_{y}^{2}}{{c}^{2}}}\right),$ and all discrete grid-points in the phase space are numbered. Likewise, in figure 2, all macroparticles are numbered. Moreover, the term 'driving strength' in figure 1(b) stands for $O(f)$ in equation (2) and ${F}_{x}$ in figure 2(b) for $x$-component of $O$ in equation (3). Parameters in figure 1 are chosen as $R=1{0}^{-6}m$ and ${f}_{{\max }}^{0}=1{0}^{25}{m}^{-3}{c}^{-3}$ (where $c$ is speed of light in vacuum).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Clearly, marked differences exist between results based on different treatments. As shown in figures 1(c), (d), the difference is so large that results in figure 1(c) is higher several orders of magnitude than those in figure 1(d). Moreover, although the curve in figure 1(c) and that in figure 1(d) are in same unit, their shapes are obviously different. Using finer mesh in phase space still faces same problem. Even though results in figures 2(c), (d) are at a same order of magnitude, their relative error is nearly $\gt 100 \% .$ Figures 1, 2 display examples in $2D2V$ format. Similar marked differences between results based on different treatments exist universally for EM problem in all formats (from the simplest $1D2V$ format to $3D3V$ one) because there is always a high-dimension matrix (due to ${E}_{{Tr}}$) no matter which of these formats is studied. For example, a typical $1D2V$ initial PDF $f{{\rm{| }}}_{t=0}={f}_{{\max }}^{0}{\exp }\left(-\frac{{x}^{2}}{{R}^{2}}\right){\exp }\left(-\left(1-v\left(t=0\right){\upsilon }_{y}\right)/\sqrt{1-\frac{{\upsilon }_{x}^{2}+{\upsilon }_{y}^{2}}{{c}^{2}}}\right),$ where $v$ is space-independent, can be represented by a mesh in the $1D2V$ phase space, and the number of points in this mesh are usually chosen to be up to $1000$ or higher. Because $v\left(t=0\right)\ne 0$ implies $f{{\rm{| }}}_{t=0}$ being an asymmetric function of ${\upsilon }_{y}$ and only a space variable $x$ exists, there are initially a current vortex ${\partial }_{x}{j}_{y}\ne 0$ and hence a $z$-direction component of $B$ and a fast growing $y$-direction component of ${E}_{{Tr}}.$ As previously mentioned, the presence of ${E}_{{Tr}}$ means a high-dimension matrix to be treated when calculating subsequent time evolution of $f.$ According to figures 3(c), (d), strict treatment and approximated one yield obviously different results of the local growth rate. This is similar to what are displayed by figures 1(c),(d). Figure 3(b) displays ${E}_{{Lo}}$-profile after two time steps $2{\rm{\Delta }}t,$ which exhibits marked space oscillation and hence implies marked density wave.

Zoom In Zoom Out Reset image size

The root reason for limiting scientific validity of older approximated methods/approaches is to adopt approximation on the coefficient matrix needing diagonalization in discrete expression of the VE, which consists of ${\partial }_{t}f$-values at ${N}_{r}\times {N}_{\upsilon }$ mesh points of $6D$ phase space:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \sum _{1\leqslant m\leqslant {N}_{r},1\leqslant n\leqslant {N}_{\upsilon }}coef\,g\left(i{\rm{\Delta }}r,j{\rm{\Delta }}\upsilon ;m{\rm{\Delta }}r,n{\rm{\Delta }}\upsilon \right){\partial }_{t}f\left(i{\rm{\Delta }}r,j{\rm{\Delta }}\upsilon ,t\right)\\ \,=\,FG\left[i{\rm{\Delta }}r,j{\rm{\Delta }}\upsilon ,f\left(i{\rm{\Delta }}r,j{\rm{\Delta }}\upsilon ,t\right)\right],\end{array}\end{eqnarray} \tag{ 4 }$$

where ${\rm{\Delta }}r$ and ${\rm{\Delta }}\upsilon$ are steps of the mesh/grid in the phase space, and detailed forms of the function ${coefg}$ and the functional ${FG}$ can be directly derived from MEs and the functional $O$ in equations (1), (2). After being aware of this root reason, we diagonalize this coefficient matrix, even though doing so is rather time-consuming, to obtain a set of equations, which are straightforward expressions of ${N}_{r}\times {N}_{\upsilon }$ functions ${\partial }_{t}f\left(r,\upsilon ,t\right)$ in terms of $f\left(r,\upsilon ,t\right):$

$$\begin{eqnarray}&&\begin{array}{l}{\partial }_{t}f\left(i{\rm{\Delta }}r,j{\rm{\Delta }}\upsilon ,t\right)=\displaystyle \sum _{1\leqslant m\leqslant {N}_{r},1\leqslant n\leqslant {N}_{\upsilon }}coefd\left(i{\rm{\Delta }}r,j{\rm{\Delta }}\upsilon ;m{\rm{\Delta }}r,n{\rm{\Delta }}\upsilon \right)\\ \,\times FD\left[m{\rm{\Delta }}r,n{\rm{\Delta }}\upsilon ,f\left(m{\rm{\Delta }}r,n{\rm{\Delta }}\upsilon ,t\right)\right],\end{array}\end{eqnarray} \tag{ 5 }$$

where detailed forms of the coefficient function ${coeffd}$ and the functional ${FD}$ are easily written out but somewhat long for being presented here. Although this straightforward expression is now available, another difficulty still hinders calculating the $f{{\rm{| }}}_{t={\rm{\Delta }}t}$-profile from the $f{{\rm{| }}}_{t=0}$-profile and the ${\partial }_{t}f{{\rm{| }}}_{t=0}$-profile (because the time-step ${\rm{\Delta }}t$ should ensure ${f}_{t={\rm{\Delta }}t}\geqslant 0\forall r,\upsilon .$) Because the discrete expression is through ${N}_{r}\times {N}_{\upsilon }$ mesh points in $6D$ phase space, the allowed value of ${\rm{\Delta }}t$ should be less than the minimum of ${N}_{r}\times {N}_{\upsilon }$ values of ${\rm{| }}\frac{0-f}{{\partial }_{t}f}{{\rm{| }}}_{r,\upsilon }{\rm{| }}.$ Usually, a comprehensive discrete expression of a $f$-profile should contain those boundary mesh points where $f=0$ satisfies. Because not all ${N}_{r}\times {N}_{\upsilon }$ values of ${\partial }_{t}f{{\rm{| }}}_{t=0}$ are $\geqslant 0,$ once at some boundary mesh points the values of ${\partial }_{t}f{{\rm{| }}}_{t=0}$ are negative, the allowed value of ${\rm{\Delta }}t$ will be confined to be $\to 0.$ This difficulty is analogous to the uniform convergence requirement mentioned before.

For the V-M simulation, two difficulties drive us to seek for new ways of describing the evolution of $f.$ As shown latter, in a geometric expression of $f,$ ${N}_{r}\times {N}_{\upsilon }$ expression formulas of ${\partial }_{t}f$ are converted into those of ${\partial }_{t}{H}_{0}$ in a same amount. Unlike $f$-values mandatory to be $\geqslant 0,$ ${H}_{0}$-value are allowed to be negative. This overcomes the difficulty arising from the requirement that ${\rm{\Delta }}t$ should warrant $f{{\rm{| }}}_{t={\rm{\Delta }}t}$ being non-negative anywhere.

A main difficulty in simulating the V-M system, both in Eulerian approach and Lagrangian one, is to find an efficient expression of $f$ free from being supplemented by unnecessary approximations. In the Eulerian approach, if $f$ is expressed by a power expansion formula

$$\begin{eqnarray}&&f\left(\upsilon ,r,t\right)=\displaystyle \sum _{i}{c}_{i}\left(r,t\right){\left[\upsilon /c-\vec{u}/c\right]}^{i},\end{eqnarray} \tag{ 6 }$$

the VE will act as a recurrence relation/formula defining a set of expansion coefficient functions $\{{c}_{i\geqslant 0}\left(r,t\right),i\geqslant 0\}.$ $f$ constructed from this $\{{c}_{{\rm{i}}}\}$-set obviously satisfy the VE. The recurrence relation/formula can lead to a general relation between two expansion coefficient functions with adjacent subindex

$$\begin{eqnarray}&&{c}_{i+1}\left(r,t\right)\sim \frac{1}{i+1}{c}_{i}^{\left(0,1\right)}\left(r,t\right){\rm{;}}\end{eqnarray} \tag{ 7 }$$

which implies ${c}_{i}\sim \frac{1}{i!}{c}_{0}^{\left(0,i\right)}.$ Here, ${c}_{0}^{\left(0,i\right)}\equiv {d}_{{t}^{i}}^{i}{c}_{0}$ and ${c}_{0}^{\left(i,j\right)}\equiv {d}_{{r}^{i}}^{i}{d}_{{t}^{j}}^{j}{c}_{0}$ are common expressions of high-order derivatives of ${c}_{0}\left(r,t\right).$ The convergence condition ${li}{m}_{i\to {\rm{\infty }}}{\rm{| }}\frac{{c}_{i+1}{\left[\upsilon -u\right]}^{i+1}}{{c}_{i}{\left[\upsilon -\vec{u}\right]}^{i}}{\rm{| }}\lt 1,$ whose detailed form is ${\rm{| }}\upsilon -\vec{u}{\rm{| }}/c\lt {li}{m}_{i\to {\rm{\infty }}}\frac{{c}_{i}}{{c}_{i+1}}\sim {li}{m}_{i\to {\rm{\infty }}}\left(i+1\right){li}{m}_{i\to {\rm{\infty }}}\frac{{c}_{0}^{\left(0,i\right)}}{{c}_{0}^{\left(0,i+1\right)}}\lt {li}{m}_{i\to {\rm{\infty }}}\left(i+1\right),$ can be warranted if ${c}_{0}$ can satisfy $c\lt {li}{m}_{i\to {\rm{\infty }}}{\rm{| }}\frac{{c}_{0}^{\left(0,i\right)}}{{c}_{0}^{\left(0,i+1\right)}}{\rm{| }},$ where $c$ is speed of light in vacuum. Namely, once ${c}_{0}$ is convergent with respect to $t,$ it can lead to a set $\{{c}_{i\gt 0}\}$ warranting the expansion series ${\sum }_{i\gt 0}{c}_{i}{\rm{| }}\upsilon -\vec{u}{{\rm{| }}}^{i}$ convergent over the region ${\rm{| }}\upsilon -\vec{u}{\rm{| }}/c\lt {\rm{\infty }}\supset {\rm{| }}\upsilon -\vec{u}{\rm{| }}/c\lt 2.$

The power expansion can be converted into another form

$$\begin{eqnarray}&&\begin{array}{l}f\left(\upsilon ,r,t\right)\equiv \displaystyle \sum _{i\geqslant 0}{c}_{i}\left(r,t\right)\cdot {\left[\upsilon /c-\vec{u}/c\right]}^{i}={b}_{0}\left(r,t\right)\cdot \displaystyle \int \cos \left(q* \left[\upsilon -\vec{u}\right]\right)dq\\ \,+\displaystyle \sum _{i}{b}_{i\geqslant 1}\left(r,t\right)\cdot {\left[\upsilon /c-\vec{u}/c\right]}^{i}={b}_{0}\left(r,t\right)\cdot \left[1-H\left(| \upsilon /c-\vec{u}/c| \right)\right]\\ \,+\displaystyle \sum _{i}{b}_{i\geqslant 1}\left(r,t\right)\cdot {\left[\upsilon /c-\vec{u}/c\right]}^{i},\end{array}\end{eqnarray} \tag{ 8 }$$

where Heviside function is defined as $H\left(x\gt 0\right)=1$ and $H\left(x\leqslant 0\right)=0$ and satisfies $2H\left({\rm{| }}x{\rm{| }}\right)=H\left(x\right)+H\left(-x\right),$ and the coefficient function $\{{b}_{i}\}$-set can be straightforward transformed from the $\{{c}_{i}\}$-set.

$$\begin{eqnarray}&&{b}_{2m}\equiv {c}_{2m}-\int \frac{{\left(-1\right)}^{m}{q}^{2m}}{2m!}{dq}{\rm{;}}\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{b}_{2m+1}\equiv {c}_{2m+1},\end{eqnarray} \tag{ 10 }$$

where $\frac{{\left(-1\right)}^{m}{q}^{2m}}{2m!}$ is from Taylor expansion of ${\cos }$ function. The non-negative-probability requirement demands ${b}_{1}\left(r,t\right)\equiv 0$ mandatory to be fulfilled [56–58]. Thus, complete microscopic description, which consists of the VE and the non-negative-probability requirement, implies following relation:

$$\begin{eqnarray}&&\begin{array}{l}f\in \{g:\hat{L}g=0,g{| }_{| \upsilon | /c=1}=0,g\in C\left({\rm{\Omega }}\right)\}{\cap }^{\,}\{g:g\geqslant 0\forall r,t,| \upsilon | /c\lt 1\}\\ \,\subset \{g:\hat{L}g=0,g{| }_{| \upsilon | /c=1}=0,g\in C\left({\rm{\Omega }}\right)\}\end{array}\end{eqnarray} \tag{ 11 }$$

where $C\left({\rm{\Omega }}\right)$ represents the set of all continuous functions defined over $7$-D space ${\rm{\Omega }}$ whose coordinates are $\upsilon ,r,t.$

Note that when the operator $\hat{L}$ acts on ${b}_{i}\cdot {\left(\upsilon -\vec{u}\right)}^{i},$ relevant terms will be in three classes: rising-order terms $\sim {\left(\upsilon -\vec{u}\right)}^{i+1},$ keeping-order terms $\sim {\left(\upsilon -\vec{u}\right)}^{i},$ and decreasing-order ones $\sim {\left(\upsilon -\vec{u}\right)}^{i-1}.$ Thus, there will be a mono-variable recurrence relation/formula chain

$$\begin{eqnarray}&&{b}_{i\geqslant 2}\left(r,t\right)=\frac{1}{i}{RF}\left(\vec{u},\vec{E},\vec{B},{b}_{i-1},{b}_{i-2}\right)=S{F}_{i}\left(\vec{u},\vec{E},\vec{B},{b}_{1},{b}_{0}\right)={F}_{i}\left(\vec{u}\left(\vec{E},\vec{B}\right),{b}_{0}\right)={H}_{i}\left(\vec{u}\left(\vec{E},\vec{B}\right)\right){\rm{;}}\end{eqnarray} \tag{ 12 }$$

where the functional ${RF}$ represents the recurrence formula reflected by the VE, and expressing $\vec{u}$ as $\vec{u}\left(\vec{E},\vec{B}\right)$ is due to the fact that ${b}_{1}\equiv 0$ can lead to an equation among $\vec{u},\vec{E},\vec{B}$ [56–58]

$$\begin{eqnarray}&&{\partial }_{t}\left[\frac{\vec{u}}{\sqrt{1-\vec{u}\cdot \vec{u}/{c}^{2}}}\right]=-\frac{e}{{m}_{e}}\left[\vec{E}+\vec{u}\times \vec{B}\right].\end{eqnarray} \tag{ 13 }$$

The relation between ${F}_{i}$ and ${H}_{i}$ is described as follows: substituting the power series expansion back into Gauss law (where these ${b}_{i}$ and $\vec{E}$ are space-time inhomogeneous functions and ${\varepsilon }_{0}$ is permittivity of free space)

$$\begin{eqnarray}&&{b}_{0}+\displaystyle \sum _{i\geqslant 2}\left[{b}_{i}\cdot {\int }_{{\rm{| }}\upsilon /c{\rm{| }}\lt 1}{\left[\upsilon /c-\vec{u}/c\right]}^{i}{d}^{3}\upsilon \right]={\varepsilon }_{0}{\rm{\nabla }}\cdot \vec{E}{\rm{;}}\end{eqnarray} \tag{ 14 }$$

and utilizing ${b}_{i}$'s expression equation (12), we can express ${b}_{0}$ in terms of $\vec{u},\vec{E},\vec{B}$ and hence convert ${F}_{i}$ into ${H}_{i}.$ Note that equation (14) is indeed an infinite-order PDE (with respect to $t$) of ${b}_{0}$ and hence often inevitably trigger approximations.

Semi-Lagrangian approaches [43–50, 52–54] face similar difficulty which refers to exact integration of contour of $f$ over $\upsilon$-space difficult to be conducted. As a result, even exact functional dependence of $f$ on $\left(\upsilon ,E,B\right)$ can be known from the contour form, it is difficult to obtain exact solution of $f$ in terms of $\left(\upsilon ,r,t\right).$

The VE reflects global conservation law $\int f{d}^{3}\upsilon {d}^{3}r={constant},$ and the non-negative-probability requirement $f\geqslant 0\forall \upsilon ,r,t$ reflects local positivity-preserving law/requirement [43–50] which is mandatory for applying probability theory to various application issues. Many authors try to take this local positivity-preserving law/requirement into account by adjusting key technical details used in Lagrangian approaches [43–50], such as using more interpolation functions and dividing a time-step into more sub-steps etc Mathematically, these effort choose to update $f$ in a more complicated manner in which the displacement of a $f$-element in $6D$ phase space is a nonlinear function of the time-step ${\rm{\Delta }}t.$ But these effort also cause computation cost increasing significantly.

Historically, in theoretical plasmas physics, investigations are merely around the VE and hence often drop in a logic contradiction: Starting from the VE,$\to$ deriving its equivalent equation ${EE}\{{VE}\},$ (where ${EE}$ is the acronym of 'equivalent equation' and the phrase 'equivalent' means no any mathematics approximation being applied and hence can share common solutions with the VE),$\to$ then truncating this equivalent equation to obtain an equation ${TEE}\{{VE}\},$ (where ${TEE}$ is the acronym of 'truncated equivalent equation' and the phrase 'truncating' means discarding some non-zero terms in the equation and hence cannot share common solutions with the equivalent equation),$\to$ then combining the amputated ${TEE}\{{VE}\}$ and the VE to obtain solutions of the VE. Obviously, such a procedure is contradict with the fact the a ${TEE}\{{VE}\}$ cannot share solutions with the equivalent equation ${EE}\{{VE}\},$ as well as the VE. Actually, people have been aware of that 'the 'truncation approximation' is a delicate operation, which might lead to a numerical failure of the solution for later times' [20]. Therefore, ignoring the non-negative-probability requirement delays a sound kinetic simulation scheme to be found.

There are different expressions of the PDF $f.$ For example, the set $\{\int {\left[\upsilon -\vec{u}\right]}^{i}f{d}^{3}\upsilon {\rm{;}}i\geqslant 0\},$ where $\vec{u}\left(r,t\right)\equiv \frac{\int \upsilon f{d}^{3}\upsilon }{\int f{d}^{3}\upsilon },$ is convenient to reflect the global conservation law, and the set $\{{f}^{\left(i\right)}{{\rm{| }}}_{\upsilon =\vec{u}}\equiv \frac{{d}^{i}f}{d{\upsilon }^{i}}{{\rm{| }}}_{\upsilon =\vec{u}}{\rm{;}}i\gt 0\}$ is convenient to reflect the local positivity-preserving law/requirement. Compared with sole usage of these expressions, combinative usage of them can avoid unnecessary approximations and is more efficient for establishing a universal and sound scheme of the EM kinetic simulation.

Among different expressions of $f,$ choosing the most suitable is important to obtain exact solutions of $f.$ An unsuitable expression will inevitably lead to unnecessary mathematical approximations supplemented. As far as $f$ is concerned, previously mentioned two Eulerian expressions: mesh expression and power series expressions are far from suitable ones. For example, to make the $\{{b}_{i}\}$-set certain needs to first make the ${b}_{0}$ certain from equation (14)), an infinite-order PDE (with respect to $t$) of ${b}_{0},$ and hence often triggers approximations. In contrast, respecting some mathematical features of $f$ and targeted introducing more suitable Eulerian expression is feasible. For example, because equation (14), or Gauss law, must be respected, the $\{{b}_{i}\}$-set expression can be further improved into a more suitable expression in which $f$ is expressed in term of $f$ itself [59]:

$$\begin{eqnarray}&&f\left(\upsilon ,r,t\right)={b}_{0}\left(r,t\right)\cdot \delta \left(\upsilon -\vec{u}\right)+{\rm{\Delta }}n\left(r,t\right)\cdot \left[{H}_{1}\left(\upsilon ,r,t\right)+{H}_{0}\left(\upsilon ,r,t\right)\right],\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Delta }}n\equiv n-{b}_{0}={\varepsilon }_{0}{\rm{\nabla }}\cdot \vec{E}-{b}_{0}$ is always $\geqslant 0,$ non-negative function $\left[{H}_{1}+{H}_{0}\right],$ which acts as a conditional probability density function (C-PDF), describes the allocation of particles contained in the ${\rm{\Delta }}n$-profile into different values of $\upsilon -\vec{u}\ne 0.$ It is straightforward, from the definition formula equation (15), to derive that the non-negative function $\left[{H}_{1}+{H}_{0}\right]$ satisfies $\int \left[{H}_{1}+{H}_{0}\right]{d}^{3}\upsilon =1,\int \left[\upsilon -\vec{u}\right]\cdot \left[{H}_{1}+{H}_{0}\right]{d}^{3}\upsilon =0,\left[{H}_{1}+{H}_{0}\right]\left(\vec{u},r,t\right)\equiv 0\equiv \left[{H}_{1}+{H}_{0}\right]\left({\rm{| }}\upsilon {\rm{| }}=c,r,t\right).$ Clearly, such a non-negative function can be constructed from functions in two different classes, one class is denoted as ${H}_{1}$ and the other ${H}_{0}$ [56]:

$$\begin{eqnarray}&&\int {H}_{1}{d}^{3}\upsilon =1,\int \left[\upsilon -\vec{u}\right]\cdot {H}_{1}{d}^{3}\upsilon =0,{H}_{1}\left(\vec{u},r,t\right)\equiv 0\equiv {H}_{1}\left({\rm{| }}\upsilon {\rm{| }}=c,r,t\right){\rm{;}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&\int {H}_{0}{d}^{3}\upsilon =0,\int \left[\upsilon -\vec{u}\right]\cdot {H}_{0}{d}^{3}\upsilon =0,{H}_{0}\left(\vec{u},r,t\right)\equiv 0\equiv {{\rm{H}}}_{0}\left({\rm{| }}\upsilon {\rm{| }}=c,r,t\right){\rm{;}}\end{eqnarray} \tag{ 17 }$$

and ${H}_{1}$ is $\geqslant 0\forall \upsilon ,r,t$ while ${H}_{0}$ is not. Moreover, it is easy to find, by applying ${\partial }_{t}$-operator and ${{\rm{\nabla }}}_{r}$-operator on the definition formulas equations (16), (17), that, $4$ functions: ${\partial }_{t}{H}_{1},$ ${{\rm{\nabla }}}_{r}{H}_{1},$ ${\partial }_{t}{H}_{0}$ and ${{\rm{\nabla }}}_{r}{H}_{0},$ (denoting them as $F$), are alike to ${H}_{0}$ or satisfy $\int F{d}^{3}\upsilon =0=\int \left[\upsilon -\vec{u}\right]* F{d}^{3}\upsilon .$ Note that any solution of equation (17) multiplying a $\upsilon$-independent function is still a solution of equation (17), but this important property does not hold for solutions of equation (16).

If an initial profile $f{{\rm{| }}}_{t=0}$ can correspond to a $\{{b}_{i\ne 0}{{\rm{| }}}_{t=0}\}$-set convergent with respect to $i,$ it will yield at least $5$ initial profiles: $n{{\rm{| }}}_{t=0},$ $\vec{u}{{\rm{| }}}_{t=0},$ ${b}_{0}{{\rm{| }}}_{t=0},$ ${H}_{1}{{\rm{| }}}_{t=0}$ and ${H}_{0}{{\rm{| }}}_{t=0}\equiv 0.$ Because $W\left(\vec{u},m\right)\equiv {\int }_{{\rm{| }}\upsilon {{\rm{| }}}^{2}\lt {c}^{2}}{\left[\upsilon -\vec{u}\right]}^{m}{d}^{3}\upsilon$ is a function of $\vec{u}$ and $m,$ ${\int }_{{\rm{| }}\upsilon {{\rm{| }}}^{2}\lt {c}^{2}}{H}_{1}{d}^{3}\upsilon$ is therefore a combination of these functions of $\vec{u}$ and $m:$

$$\begin{eqnarray}&&{\int }_{{\rm{| }}\upsilon {{\rm{| }}}^{2}\lt {c}^{2}}{H}_{1}{d}^{3}\upsilon =\displaystyle \sum _{m\geqslant 2}^{{\rm{\infty }}}{\tilde{b}}_{m\geqslant 2}{\int }_{{\rm{| }}\upsilon {{\rm{| }}}^{2}\lt {c}^{2}}{\left[\upsilon -\vec{u}\right]}^{m}{d}^{3}\upsilon =\displaystyle \sum _{m\geqslant 2}^{{\rm{\infty }}}\left[{\tilde{b}}_{m\geqslant 2}\cdot \vec{W}\left(\vec{u},m\right)\right],\end{eqnarray} \tag{ 18 }$$

where ${\tilde{b}}_{m\geqslant 2}\equiv \frac{{b}_{m\geqslant 2}{\rm{| }}\left(r,t\right)}{{\rm{\Delta }}n\left(r,t\right)},$ and its value equaling to $1$ implies an equation of combination coefficients ${\tilde{b}}_{m\geqslant 2}.$ Likewise, ${\int }_{{\rm{| }}\upsilon {{\rm{| }}}^{2}\lt {c}^{2}}\left[\upsilon -\vec{u}\right]* {H}_{1}{d}^{3}\upsilon =0,$ ${H}_{1}\left(\vec{u},r,t\right)\equiv 0\equiv {H}_{1}\left({\rm{| }}\upsilon {\rm{| }}=c,r,t\right)$ and equation (18) also imply equation of combination coefficients ${\tilde{b}}_{m\geqslant 2}.$ If expressing these combination coefficients in terms of power series of $\tilde{\vec{u}}\equiv \vec{u}/c:$ $\{{\tilde{\vec{u}}}^{0},{\tilde{\vec{u}}}^{1},...,{\tilde{\vec{u}}}^{n},...\},$ equations (14), (16) will imply that the ${H}_{1}{{\rm{| }}}_{t=0}$-profile naturally defines a matrix $M$ linking the set $\{...,{\tilde{b}}_{i\geqslant 2},...\}$ with the set $\{{\tilde{\vec{u}}}^{0},{\tilde{\vec{u}}}^{1},...,{\tilde{\vec{u}}}^{n},...\}$

$$\begin{eqnarray}&&{\tilde{b}}_{i\geqslant 2}{{\rm{| }}}_{t=0}=\displaystyle \sum _{j\geqslant 0}{M}_{i-2+1,j}{\tilde{\vec{u}}}^{j}{{\rm{| }}}_{t=0}.\end{eqnarray} \tag{ 19 }$$

Now that ${\sum }_{i\geqslant 2}{\tilde{b}}_{i}\left(\vec{u}\left(r,t=0\right)\right)* {\left[\upsilon -\vec{u}\left(r,t=0\right)\right]}^{i}$ can satisfy equations (14), (16), (17) at $t=0,$ it is easy to prove its counterpart, which is to replace $\vec{u}\left(r,0\right)$ with $\vec{u}\left(r,{\rm{t}}\right),$ to still satisfy equations (14), (16), (17) at $t={\rm{\Delta }}t.$ We term it as translation operation on the ${H}_{1}{{\rm{| }}}_{t=0}$-profile and denote it as $\hat{{TL}}\left({\rm{\Delta }}t\right){H}_{1}{{\rm{| }}}_{t=0}$ in which $\vec{u}\left(r,0\right)$ in ${H}_{1}{{\rm{| }}}_{t=0}$ is replaced by $\vec{u}\left(r,{\rm{\Delta }}t\right),$ but the matrix linking $\{...,{\tilde{b}}_{i\geqslant 2}{{\rm{| }}}_{t={\rm{\Delta }}t},...\}$ and $\{...,{\tilde{\vec{u}}}^{m}\left(r,{\rm{\Delta }}t\right),...\}$ is as same as that linking $\{...,{\tilde{b}}_{i\geqslant 2}{{\rm{| }}}_{t=0},...\}$ and $\{...,{\tilde{\vec{u}}}^{m}\left(r,0\right),...\}.$ Such an operation can be expressed via a formula

$$\begin{eqnarray}&&\hat{{TL}}\left({\rm{\Delta }}t\right)g\left(\upsilon ,\vec{u}\left(r,0\right)\right)=\lambda \left({\rm{\Delta }}t\right)* g\left(\upsilon ,\vec{u}\left(r,{\rm{\Delta }}t\right)\right),\end{eqnarray} \tag{ 20 }$$

and the function $\lambda \left({\rm{\Delta }}t\right)\equiv 1$ if $g$ represents ${H}_{1}$ (because $\int {H}_{1}{{\rm{| }}}_{t=0}{d}^{3}\upsilon \equiv 1$ and $\int \hat{{TL}}\left({\rm{\Delta }}t\right){H}_{1}{{\rm{| }}}_{t=0}{d}^{3}\upsilon \equiv 1$ do not allow $\lambda \left({\rm{\Delta }}t\right)\ne 1$ ), but $\lambda \left({\rm{\Delta }}t\right)$ can be not constant-valued if $g$ represents ${H}_{0}$ (because $\int {H}_{0}{{\rm{| }}}_{t=0}{d}^{3}\upsilon \equiv 0$ and $\int \hat{{TL}}\left({\rm{\Delta }}t\right){H}_{0}{{\rm{| }}}_{t=0}{d}^{3}\upsilon \equiv 0$ allow $\lambda \left({\rm{\Delta }}{\rm{t}}\right)\ne 1$ ).

The matrix linking the $\{{\tilde{b}}_{i\geqslant 2}\}$-set and the $\{{\tilde{\vec{u}}}^{m}\}$-set reflects geometric/toplogical nature/shape of the ${H}_{1}$-profile, and likewise the ${H}_{0}$-profile also has a matrix reflecting its geometric/toplogical nature/shape. If this matrix is unchanged during evolution, the related profile is viewed as keeping geometric/toplogical nature/shape. Thus, the phrase 'deformation' refers to this matrix varying with respect to $t$ and hence the geometric/toplogical nature/shape is not fixed relative to an observer in the $u$-frame. Namely, the ${H}_{0}$-profile can allow deformation, (relative to an observer in the $u$-frame), varying with respect to $t.$

Because of $\int {H}_{0}{d}^{3}\upsilon =0,$ ${H}_{0}$ is therefore not a function positivity-preserving, and hence has zeros, beside $\upsilon =\vec{u}$ and ${\rm{| }}\upsilon {{\rm{| }}}^{2}={c}^{2},$ in $\upsilon$-space. Thus, a typical form of ${H}_{0},$ or that of functions satisfying equation (17), is

$$\begin{eqnarray}&&{H}_{0}^{{typical}}=\displaystyle \prod _{m}{\left[\upsilon /c-{Zer}{o}_{m}\left(\vec{u}\left(r,t\right)\right)/c\right]}^{{Re}{p}_{m}}\cdot \left[{\upsilon }^{2}/{c}^{2}-1\right]{\rm{;}}\end{eqnarray} \tag{ 21 }$$

each root ${\rm{| }}{Zer}{o}_{m}/c{\rm{| }}\lt 1$ is mono-variable functional of $\vec{u}$ and each number ${Re}{p}_{m}$ is repetition index of the root ${Zer}{o}_{m}$ (and at least a root satisfies ${Zer}{o}_{n}\left(\vec{u}\right)\equiv \vec{u}$). Such a product ${\prod }_{m}$ represents a class of $\upsilon$-space geometric/toplogic nature/shape. Any $\upsilon$-independent function multiplying this typical form still satisfies equation (17). From

$$\begin{eqnarray}&&\int {H}_{0}^{{typical}}{d}^{3}\upsilon =0{\rm{;}}\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&\int \left[\upsilon -\vec{u}\right]{H}_{0}^{{typical}}{d}^{3}\upsilon =0{\rm{;}}\end{eqnarray} \tag{ 23 }$$

we can find

$$\begin{eqnarray}&&0={\partial }_{t}0={\partial }_{t}\int {H}_{0}^{{typical}}{d}^{3}\upsilon =\int {\partial }_{t}{H}_{0}^{{typical}}{d}^{3}\upsilon {\rm{;}}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray*}&&0={\partial }_{t}0={\partial }_{t}\int \left[\upsilon -\vec{u}\right]{H}_{0}^{{typical}}{d}^{3}\upsilon =\int \left[-{\partial }_{t}\vec{u}\right]{H}_{0}^{{typical}}{d}^{3}\upsilon \end{eqnarray*}$$

$$\begin{eqnarray*}&&+\int \left[\upsilon -\vec{u}\right]{\partial }_{t}{H}_{0}^{{typical}}{d}^{3}\upsilon =-{\partial }_{t}\vec{u}\int {H}_{0}^{{typical}}{d}^{3}\upsilon \end{eqnarray*}$$

$$\begin{eqnarray}&&+\int \left[\upsilon -\vec{u}\right]{\partial }_{t}{H}_{0}^{{typical}}{d}^{3}\upsilon =0+\int \left[\upsilon -\vec{u}\right]{\partial }_{t}{H}_{0}^{{typical}}{d}^{3}\upsilon {\rm{;}}\end{eqnarray} \tag{ 25 }$$

and (in similar procedure)

$$\begin{eqnarray}&&0={{\rm{\nabla }}}_{r}0={{\rm{\nabla }}}_{r}\int {H}_{0}^{{typical}}{d}^{3}\upsilon =\int {{\rm{\nabla }}}_{r}{H}_{0}^{{typical}}{d}^{3}\upsilon {\rm{;}}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray*}&&0={{\rm{\nabla }}}_{r}0={{\rm{\nabla }}}_{r}\int \left[\upsilon -\vec{u}\right]{H}_{0}^{{typical}}{d}^{3}\upsilon =\int \left[-{{\rm{\nabla }}}_{r}\vec{u}\right]{H}_{0}^{{typical}}{d}^{3}\upsilon \end{eqnarray*}$$

$$\begin{eqnarray*}&&+\int \left[\upsilon -\vec{u}\right]{{\rm{\nabla }}}_{r}{H}_{0}^{{typical}}{d}^{3}\upsilon =-{{\rm{\nabla }}}_{r}\vec{u}\int {H}_{0}^{{typical}}{d}^{3}\upsilon \end{eqnarray*}$$

$$\begin{eqnarray}&&+\int \left[\upsilon -\vec{u}\right]{{\rm{\nabla }}}_{r}{H}_{0}^{{typical}}{d}^{3}\upsilon =0+\int \left[\upsilon -\vec{u}\right]{{\rm{\nabla }}}_{r}{H}_{0}^{{typical}}{d}^{3}\upsilon .\end{eqnarray} \tag{ 27 }$$

Because of

$$\begin{eqnarray}&&{\partial }_{t}{H}_{0}^{{typical}}={\partial }_{\vec{u}}{H}_{0}^{{typical}}\cdot {\partial }_{t}\vec{u}{\rm{;}}\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{{\rm{\nabla }}}_{r}{H}_{0}^{{typical}}={\partial }_{\vec{u}}{H}_{0}^{{typical}}\cdot {{\rm{\nabla }}}_{r}\vec{u}{\rm{;}}\end{eqnarray} \tag{ 29 }$$

Equations (24)–(29) will yield

$$\begin{eqnarray}&&\int {\partial }_{\vec{u}}{H}_{0}^{{typical}}{d}^{3}\upsilon =0{\rm{;}}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\int \left[\upsilon -\vec{u}\right]{\partial }_{u}{H}_{0}^{{typical}}{d}^{3}\upsilon =0.\end{eqnarray} \tag{ 31 }$$

If we use the symbol '$\{{H}_{0}\}$' to represent the set of functions satisfying equation (17), equations (30), (31) imply that if ${H}_{0}^{{typical}}\in \{{H}_{0}\}$ exists, so does ${\partial }_{\vec{u}}{H}_{0}^{{typical}}\in \{{H}_{0}\}.$

The VE is indeed an expression of the local growth rate ${\partial }_{t}\left[{H}_{1}+{H}_{0}\right]$ [59]:

$$\begin{eqnarray}&&{\partial }_{t}\left[{H}_{1}+{H}_{0}\right]=\hat{{GR}}\left[{H}_{1}+{H}_{0}\right]{\rm{;}}\end{eqnarray} \tag{ 32 }$$

where the effect of the operator $\hat{{GR}}$ on a function $F$ obey following formula [59]

$$\begin{eqnarray}&&\hat{{GR}}F\equiv -\left[\upsilon \cdot {\rm{\nabla }}-e\left[\vec{E}+\upsilon \times \vec{B}\right]\cdot {\partial }_{p}\right]F+\frac{\left[{\partial }_{t}+\upsilon \cdot {\rm{\nabla }}\right]\left[n-{n}_{0}\right]}{n-{n}_{0}}F+K\left(\upsilon ,r,t\right){\rm{;}}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&K\left(\upsilon ,r,t\right)\equiv \frac{\left[\upsilon \cdot {\rm{\nabla }}u\right]{n}_{0}}{n-{n}_{0}}{\partial }_{\upsilon }\delta \left(\upsilon -\vec{u}\right).\end{eqnarray} \tag{ 34 }$$

Note that a fraction of the local growth rate

$$\begin{eqnarray}&&\hat{{TLGR}}F{{\rm{| }}}_{t=0}\equiv {li}{m}_{{\rm{\Delta }}t\to 0}\frac{\hat{{TL}}\left({\rm{\Delta }}t\right)F{{\rm{| }}}_{t=0}-F{{\rm{| }}}_{t=0}}{{\rm{\Delta }}t},\end{eqnarray} \tag{ 35 }$$

is responsible for $\hat{{TL}}\left(t\right)F{{\rm{| }}}_{t=0},$ and another fraction

$$\begin{eqnarray}&&\hat{{DEGR}}F{{\rm{| }}}_{t=0}\equiv \left[\hat{{GR}}-\hat{{TLGR}}\right]F{{\rm{| }}}_{t=0}\end{eqnarray} \tag{ 36 }$$

is responsible for the deformation of $F$ relative to an observer in the $u$-frame.

Because ${H}_{0}$ must satisfy equation (17), its general form reads $\lambda \left(t\right){H}_{0}^{{typical}}.$ Substituting this general form into the VE, we have

$$\begin{eqnarray*}&&{\partial }_{t}{H}_{0}={d}_{t}\lambda \cdot {H}_{0}^{{typical}}+\lambda {\partial }_{t}\vec{u}\cdot {d}_{\vec{u}}{H}_{0}^{{typical}}=\hat{{DEGR}}{H}_{1}\end{eqnarray*}$$

$$\begin{eqnarray*}&&+\lambda \hat{{GR}}{H}_{0}^{{typical}}=\hat{{DEGR}}{H}_{1}+\lambda \vec{W}\left(\upsilon ,r,t\right)\cdot {d}_{\vec{u}}{H}_{0}^{{typical}}\end{eqnarray*}$$

$$\begin{eqnarray}&&+\lambda S\cdot {H}_{0}^{{typical}}{\rm{;}}\end{eqnarray} \tag{ 37 }$$

or

$$\begin{eqnarray}&&\frac{{d}_{t}\lambda }{\lambda }{H}_{0}^{{typical}}+\vec{W}\cdot {d}_{\vec{u}}{H}_{0}^{{typical}}+S\cdot {H}_{0}^{{typical}}=\frac{1}{\lambda }\hat{{DEGR}}{H}_{1},\end{eqnarray} \tag{ 38 }$$

where

$$\begin{eqnarray}&&\vec{W}\left(\upsilon ,r,t\right)\equiv {\partial }_{t}\vec{u}-\frac{1}{{m}_{e}}\left[\upsilon \cdot {\rm{\nabla }}\vec{u}+\left[{F}_{\upsilon }{{\rm{\Gamma }}}_{\upsilon }^{-3}-{F}_{\vec{u}}{{\rm{\Gamma }}}_{\vec{u}}^{-3}\right]\right]{\rm{;}}\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{\vec{F}}_{\upsilon }\equiv e\left[\vec{E}+\upsilon \times \vec{B}\right]{\rm{;}}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\upsilon }\equiv \frac{1}{\sqrt{1-{\left(\upsilon /c\right)}^{2}}}{\rm{;}}\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&S\left(\upsilon ,r,t\right)\equiv \frac{\left[{\partial }_{t}+\upsilon \cdot {\rm{\nabla }}\right]\left[n-{n}_{0}\right]}{n-{n}_{0}}{\rm{;}}\end{eqnarray} \tag{ 42 }$$

and the functional $K$ defined by equation (34) has been contained in $\hat{{GR}}{H}_{1},$ as well as in $\hat{{DEGR}}{H}_{1}.$

Because of equations (16), (17), (32)–(36), it is easy to find that

$$\begin{eqnarray}&&\int \{\hat{{DEGR}}\left[{H}_{1}+{H}_{0}\right]\}{d}^{3}\upsilon =0{\rm{;}}\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&\int \left[\upsilon -\vec{u}\right]\{\hat{{DEGR}}\left[{H}_{1}+{H}_{0}\right]\}{d}^{3}\upsilon =0{\rm{;}}\end{eqnarray} \tag{ 44 }$$

, or $\hat{{DEGR}}\left[{H}_{1}+{H}_{0}\right]\in \{{H}_{0}\}.$ Therefore, the VE reflects a relation among these members of the set $\{{H}_{0}\},$ and has been converted into equation (38), a formal $1$st-order ordinary differential equation (ODE) of ${H}_{0}^{{typical}}$ with respect to variable $\vec{u}.$ Its formal strict solution

$$\begin{eqnarray}&&{H}_{0}^{{typical}}={\exp }\left[-\int \frac{\frac{{d}_{t}\lambda }{\lambda }+S}{\vec{W}}d\vec{u}\right]\int \frac{1}{\lambda }\Space{0ex}{0.25ex}{0ex}\{\frac{\hat{{DEGR}}{H}_{1}}{\vec{W}}{\exp }\left[\int \frac{\frac{{d}_{t}\lambda }{\lambda }+S}{\vec{W}}d\vec{u}\right]\Space{0ex}{0.25ex}{0ex}\}d\vec{u}{\rm{;}}\end{eqnarray} \tag{ 45 }$$

reflects a relation between ${H}_{0}^{{typical}},$ a member of $\{{H}_{0}\},$ and another member $\hat{{DEGR}}{H}_{1}.$

Asymptotic behavior of ${H}_{0}$ at $t=0$ are $0={H}_{0}{{\rm{| }}}_{t=0}=\{\lambda * {H}_{0}^{{typical}}\}{{\rm{| }}}_{t=0}$ and $0\ne {\partial }_{t}\{\lambda * {H}_{0}^{{typical}}\}{{\rm{| }}}_{t=0},$ and lead to

$$\begin{eqnarray}&&\lambda \left(t=0\right)=0{\rm{;}}\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{d}_{t}\lambda \left(t=0\right)=1,\end{eqnarray} \tag{ 47 }$$

which can easily lead to

$$\begin{eqnarray}&&\lambda =\frac{{\sin }\left(\nu t\right)}{\nu },\end{eqnarray} \tag{ 48 }$$

where $\nu$ is a constant corresponding to the minimum of all time scales (each represents the time for $\left[{H}_{1}+{H}_{0}\right]$ at a given coordinate $\left(\upsilon ,r\right)$ decreasing to $0:$

$$\begin{eqnarray}&&\nu ={Max}\Space{0ex}{0.25ex}{0ex}\{\Space{0ex}{0.25ex}{0ex}| \frac{\left[{\partial }_{t}{H}_{0}-\hat{{TLGR}}{H}_{0}\right]{{\rm{| }}}_{t=0}}{{H}_{1}{{\rm{| }}}_{t=0}+{H}_{0}{{\rm{| }}}_{t=0}-0}\Space{0ex}{0.25ex}{0ex}| \Space{0ex}{0.25ex}{0ex}\}={Max}\Space{0ex}{0.25ex}{0ex}\{\Space{0ex}{0.25ex}{0ex}| \frac{\left[{\partial }_{t}{H}_{0}\right]{{\rm{| }}}_{t=0}}{{H}_{1}{{\rm{| }}}_{t=0}}\Space{0ex}{0.25ex}{0ex}| \Space{0ex}{0.25ex}{0ex}\}.\end{eqnarray} \tag{ 49 }$$

Such a choice of $\nu$ is for warranting ${H}_{1}+{H}_{0}={H}_{1}+\lambda \left(t\right){H}_{0}^{{typical}}$ non-negative always. Thus, equations (45), (48), (49) reveal that the C-PDF ${H}_{1}+{H}_{0}$ corresponds to a time-dependent global oscillation of ${H}_{0}$ around ${H}_{0}\equiv 0\forall \upsilon ,r,t,$ or that of ${H}_{1}+{H}_{0}$ around $\hat{{TL}}{H}_{1}.$ Because these fixed points, which correspond to ${H}_{0}=0,$ exist, the collective mode exhibits as a 'standing-wave' in $6D$ phase space.

Considering the phrase 'Eulerian' is often for describing $f$ by its values at a mesh of the $6D$ phase space, we term above procedure as 'Eulerian-like' because it is based on the description of $f$ by specified functions defined over the $6D$ phase space. These specified functions reflect geometric/toplogial nature of initial $f$-profile. The well-known drawback of the Eulerian approach, i.e., too huge mount of mesh points and corresponding discrete values, is avoided by high-efficiently expressing them through collective mode, represented by the ${H}_{0}$-profile, in the $6D$ phase space. Namely, the Eulerian approach is in a one-by-one manner, to update a mesh ${f}_{{Ir},I\upsilon }^{n}$ to another mesh ${f}_{{Ir},{\rm{I}}\upsilon }^{n+1}$ through a mesh ${\partial }_{t}{f}_{{Ir},I\upsilon }^{n}.$ Each mesh has its own geometric/tpolgical nature. In contrast, our approach is to utilize respective geometric/topologic nature of three meshes and high-efficiently express the difference ${f}_{{Ir},I\upsilon }^{n+1}-{f}_{{Ir},I\upsilon }^{n}$ as a collective mode ${H}_{0}^{n+1}-{H}_{0}^{n},$ which is determined by ${H}_{1}^{n}.$ One-by-one mesh points calculation on individual $f$-values is replaced by calculating amplitude of the collective mode. Moreover, compared with various ansatz of $f$ [40] used in simulating the V-M system, $f=n* \left[{H}_{1}+{H}_{0}\right]$ is a strict expression. On the other hand, the Lagrangian approach is to express the initial-value problem of the V-M system or the N-M system as those of many initial mono-variable functions (of $t$), or Lagrangian particles, identified by different points on a $6D$ $\left({r}_{s},{\upsilon }_{s}\right)$-profile. By now, the calculation in the Lagrangian approach is rarely to take into account the geometrical/toplogical nature of this $6D$ initial-values-profile, and to choose a rough treatment in which all initial-value problems of these Lagrangian particles are uniformly/equally treated as having a same time-convergence. This might be an unavoidable difficulty for the Lagrangian approach.

Here, we present some examples of application of this new simulation method. It simulates an electron bunch (e-bunch) travelling within a space inhomogeneous DC-fields configuration. The DC-fields configuration is described as follows: ${\vec{E}}_{{DC}}={E}_{0}\mathop{{e}_{x}}\limits^{\leftarrow}$ and ${\vec{B}}_{{DC}}=\vec{B}\left(x\right)\mathop{{e}_{z}}\limits^{\leftarrow}={B}_{0}\frac{x}{w}\mathop{{e}_{z}}\limits^{\leftarrow}$ if ${\rm{| }}x{\rm{| }}\lt w$ and ${\vec{B}}_{{DC}}={B}_{0}{sign}(x)\mathop{{e}_{z}}\limits^{\leftarrow}$ elsewhere. Namely, ${\rm{| }}{\vec{B}}_{{DC}}$ drops from ${B}_{0}$ at $x=w$ to $-{B}_{0}$ at $x=-w.$ The initial PDF of the e-bunch is $f{{\rm{| }}}_{t=0}={\exp }\left(-\frac{{\left(x-{x}_{0}\right)}^{2}+{\left(y-{y}_{0}\right)}^{2}}{{R}_{T}^{2}}-\frac{{z}^{2}}{{R}_{L}^{2}}\right){\exp }\left(-\frac{1}{\sqrt{1-{\rm{| }}\upsilon /c{{\rm{| }}}^{2}}}+\beta \cdot \frac{{\upsilon }_{z}/c}{\sqrt{1-{\rm{| }}\upsilon /c{{\rm{| }}}^{2}}}\right),$ where ${R}_{T}=30\mu m$ and ${R}_{L}=10\mu m$ are transverse and longitudinal radius respectively and $\beta c=0.05c$ is the global moving velocity (along the $z$-direction) of the e-bunch, ${x}_{0}=5\mu m$ and ${y}_{0}=5\mu m$ are initial transverse coordinates of the symmetric axis of the e-bunch relative to the $z$-axis. Note that this $f{{\rm{| }}}_{t=0}$-profile corresponds to initial current profile having a vortex ${\rm{\nabla }}\times \vec{j}{{\rm{| }}}_{t=0}\ne 0.$ As described in last subsection, we can derive, from this initial PDF, $5$ derivative initial profiles and calculate their evolutions. Some examples are presented in figures 4, 5. Velocity subspace boundary conditions are ${H}_{0}{{\rm{| }}}_{{v}^{2}={c}^{2}}=0={H}_{1}{{\rm{| }}}_{{v}^{2}={c}^{2}},$ ${H}_{0}{{\rm{| }}}_{v=u}=0={H}_{1}{{\rm{| }}}_{v=u}.$ ${H}_{0}$ and ${H}_{1}$ are finite-valued and ${H}_{1}$ are always non-negatively-valued. These examples well illustrate that phase space details of the evolving $f$-profile can be well reflected by those of the evolving ${H}_{0}$-profile.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For comparison with other kinetic simulation schemes, we present, in figure 6, results of applying this method to laser wakefield, a topic attracting many PIC simulations [60, 61]. Electronic density wake of a laser with a triangular-function-type vector potential initial envelope $a\left(t,r\right)={a}_{i}{\cos }\left(0.5\pi t/{\tau }_{i}\right){\cos }\left(0.5\pi r/{\sigma }_{i}\right),$ which is displayed in figure 6(a), has a shape similar to that displayed in figure 1 of [60], For a laser with a Gaussian-function-type initial envelope $a={a}_{0}{\exp }\left(-{z}^{2}/2{L}^{2}\right){\exp }\left(-{x}^{2}/{r}_{0}^{2}\right),$ electronic density wake displayed in figure 6(c) has a shape similar to that displayed in figure 19 of [61].

Zoom In Zoom Out Reset image size

Because of incomplete information disclosure in these refs, such as missing detailed forms of the interpolation functions used in PIC simulations, it is only possible to make a qualitative comparison. Our results presented in figure 6 also display 'bubble structure' displayed in figure 19 of [61]. But the depth and longitudinal position of density valley differ from those in figure 19 of [61]. For example, in the upper-left panel, the region $5\lt z/{lambda}\lt 10,x/{lambda}0$ is a density valley, it is spaced by a density peak from another density valley at the region $15\lt z/{lambda}\lt 25,x/{lambda}0.$ Likewise, in the down-left panel, a large density valley is at the region $25\lt z/{lambda},$ another small valley is at the region $15\lt z/{lambda}\lt 16.$

But in-depth understanding the unavoidable difficulty in particle simulation and seeking for possible solution for it are necessary. A possible solution might arise from a subtle, dual role of an electron. Each electron can be both the source of a scalar potential and that of a vector potential. A notable fact is that scalar potentials from different source electrons will not offset mutually but vector potentials will. This fact implies that the microscopic vector field exerted by an electron can be divided into two parts: ${Vec}={Ve}{c}_{c}+{Ve}{c}_{o},$ where the subscript 'c' represents 'contributed' (to macroscopic vector field) and 'o' represents 'offset'. For all electrons, there are at least a scalar field and two vector fields defined over the initial phase space $\left({r}_{s},{\upsilon }_{s}\right)$ giving a description on the system. One vector field is from those ${Ve}{c}_{o}$ and the other from those ${Ve}{c}_{c}.$ For example, accelerating rate of each electron $a\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right)$ is a vector, and hence there should be

$$\begin{eqnarray}&&0=\int {a}_{o}\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right)d{r}_{s}d{\upsilon }_{s}=\int {a}_{o}\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right)d{\upsilon }_{s}{\rm{;}}\end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray}&&{a}_{o}\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right)=\frac{e}{{m}_{e}}\left[{d}_{t}r\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right)-u\left(r\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right),t\right)\right]\times B\left(r\left(t{\rm{;}}{r}_{s},{\upsilon }_{s}\right),t\right)\end{eqnarray} \tag{ 51 }$$

These unique properties of the ${Ve}{c}_{o}$-field can be utilized for solving difficulties mentioned before.

Note that the MEs do not couple with ${Ve}{c}_{o}.$ The uniform convergence requirement forces us to re-consider the role of those RNEs. It is appropriate to view them as a definition of respective ${a}_{o},$ or all those RNEs yield the expression of an ${a}_{o}$-field, which does not couple with the MEs. In other words, the particle simulation can be in a more efficient scheme if noticing the division ${Vec}={Ve}{c}_{c}+{Ve}{c}_{o}$ and the fact that the MEs do not couple with ${Ve}{{\rm{c}}}_{o}.$ Ignoring them will make the particle simulation entangled by exactly solving ${a}_{o}$ from those RNEs. Once the ${a}_{o}$ field cannot be exactly solved, the whole particle simulation is ruined. Actually, it is unnecessary to bind the particle simulation with this nearly impossible task (because of the uniform convergence requirement). The condition for $N$ expression formula of ${a}_{o}$ able to be simultaneously solved, or the fact that those ${Ve}{c}_{c}$ (such as ${\upsilon }_{c}\equiv u,$ ${E}_{{Lo}},$ ${E}_{{Tr}}$ and $B$) allow so diverse, even opposite, values of ${a}_{o},$ lead to an equation of these ${Ve}{c}_{c},$ and this equation couples with the MEs to form a natural closure of these ${Ve}{c}_{c}$ (because the MEs do not couple with ${Ve}{c}_{o}$) [56].