Phase Space Dynamics of Relativistic Particles

Abstract

By analyzing bottleneck of numerical study on six-dimensional (6D) phase space dynamics of electron beam, we present a universal and practical scheme of exactly simulating the 3D dynamics through available computer group condition. In this scheme, the exact 6D phase space dynamics is warranted by exact solutions of 3D self-consistent fields of electron beam.

1. Introduction

The dynamics of electron beam in various man-made apparatus is a kernel content of accelerator/beam physics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Such a beam can be described by well-known Vlasov-Maxwell (V-M) equations [1, 2, 3, 12, 13, 14]. Because usually the number of realistic particles (each has a mass m_e and a charge e) in an electron beam is at astronomical figure level, it is really a huge challenge to calculate realistic 3D dynamics of such a beam even through the most advanced computer groups available currently.

This can be well illustrated by following rough estimation. An electron needs to be described by six real-valued variables (for its three position components and three velocity components), and each real-valued variable will cost ~24 byte storage mount. For a realistic beam with a total charge at pC level, the number of realistic particles it contains will be about 10^{−12 + 19}/1.6~6.3 × 10⁷ and hence corresponds to 1.5 GB storage mount. If numerical experiment is conducted through a single-CPU computer, according to the most advanced personal computer currently available, 1.5 GB can be contained in the memory which is commercially available at 4 GB level. Thus, if all data can be stored in the memory, time-consuming data exchange between the memory, a high-speed storage medium, and the disk, a low-speed storage medium, can be avoided and the time cost of updating these data per round will be determined by the time cycle of CPU, which is at ns level. Because updating a real-valued variable will spend several cycles (and hence tens ns), the time cost for updating 1.5 GB data recording 6.3 × 10⁷ real-valued variables will be at minute level (if not taking into account the time cost for updating fields information). Because the time step is chosen to be sufficiently small even for explicit difference scheme, usually at least updating hundreds of times is needed for obtaining a snapshot of the beam at a specified time point. This implies at least hour-level time cost through the above-described most advanced computer condition currently. Of course, the larger the total charge is, the more the time cost is (if the data for the beam is still able to be contained in a memory).

Above discussion has clearly revealed that only for low-charge beam, the numerical experiment, which is based on explicit difference scheme, can be conducted on a single-CPU computer. Therefore, for more general cases with high-charge beam, people have to resort to computer group and parallel code. In such a situation, if data are not exchanged between high-speed memory and low-speed disk, the bottleneck of the experiment is not the time cycle of a CPU but that of the main board which might be several folds of that of the CPU and affect the dispatch of data bus responsible for data dispatch among CPUs. Sufficient CPUs can offset the effect of the time cycle of the main board. But the fact that the larger the total charge is, the more the time cost is still exists. Thus, for nC-level beam and stricter implicit-difference-scheme of updating fields, CPUs’ number in thousand levels is still unable to do as one would wish.

The contradiction between feasible computers’ condition and large amount of data describing realistic beams often force researchers to make some compromises. Merging N > 1 realistic particles (each is of a mass m_e and a charge e) into a macroparticle (of a mass Nm_e and a charge Ne) is a common compromise. Merging can cause less data mount for describing a beam and hence can speed up numerical experiment. Therefore, it is broadly adopted in popular scheme of simulating particle dynamics, such as particle-in-cell (PIC) simulation [15, 16].

Although merging can speed up the numerical experiment, it is at the cost of the reliability of the experiment. Let us analyze what will happen when N = 2. If there are l realistic particles in a beam, describing the beam will need 6l functions: {x_i (t), y_i (t), z_i (t), d_tx_i (t), d_ty_i (t), d_tz_i (t)} and 1 =  < i =  < l. Usually, the merging is based on the following principle: at t = 0, each realistic particle merges the nearest neighbor particle. Thus, no matter how large the difference between the initial velocities of any two merged particles, denoted as d_tx, y, z _2k (t = 0) − d_tx, y, z _2k + 1(t = 0), is, d_tx, y, z _2k  − d_tx, y, z _2k + 1 is always taken as 0 or the 2kth realistic particle is always bound with the 2k + 1th one (their relative motion is ignored completely). This is a distorted description. The 2kth realistic particle and the 2k + 1th one, if they have the same initial positions x, y, z _2k (t = 0) = x, y, z _2k + 1(t = 0), tend to move to different destinations at the next time point 0 + Δt (where Δt is the time step), because of their different initial velocities d_tx, y, z _2k (t = 0) and d_tx, y, z _2k + 1(t = 0). Namely, actual picture should be described as follows: at every time point t, a macroparticle will break into N equal pieces having their own destinations at t + Δt, rather than globally move into a destination at t + Δt or being equivalently expressed as: each macroparticle at t is a gathering of N equal pieces from different sources at t − Δt. Because of the distorted description, we briefly summarize the merging as the rigid-macroparticle approximation which refers to a macroparticle having a fixed mass Nm_e and a fixed charge Ne.

The direct result of the rigid-macroparticle approximation is that the charge density profile n of the beam at next time point t + Δt is distorted, and hence the self-consistent electric field E (meeting ∇ · E = ne) is inexact. If the dynamics calculation cannot be satisfactory even over a short time scale down to the time step Δt, it is hopeless to expect calculation over longer time scale which is multifold of Δt to be reliable.

Actually, even if computer condition is sufficiently advanced to directly calculate realistic particles, Newton equations (NEs) and Maxwell equations (MEs), complicated mathematical relation between E and B, which is represented by

$${\partial }_{t}E+\nabla \times B={\mu }_{0}j,$$ E1

leads to considerable difficulty in setting up difference scheme of each NE

$$d_{\mathit{tt}}{r}_i\left(t\right)=\frac{e}{m_e}\left[E\left(r_i\left(t\right),t\right)+d_t{r}_i\left(t\right)\times B\left(r_i\left(t\right),t\right)\right].$$ E2

This is because that at any time t, although E can be known from all particles’ information at t, B cannot. Due to the term ∂ _tE in Eq. (1), the information of B at t have dependence on all particles’ information at t + Δt. This implies that the NE is a complicated differential equation of r_i (t), and its difference scheme is difficult to be established. Namely, each difference version of Newton equation involves indeed all particles’ coordinates at t + Δt, i.e., {r_i (t + Δt); 1 =  < i =  < N}. Thus, N difference versions form a linear equation set of {r_i (t + Δt); 1 =  < i =  < N} whose solutions involve inevitably a N × N matrix. The larger N is, the more hopeless exact solution is.

Clearly, an exact or reliable dynamics calculation on any V-M system needs more in-depth exploring its universal properties, which can be favorable to avoid above-mentioned contradictions. In the following paragraph, we reveal some universal properties of the V-M system. Based on these universal properties, we establish a universal and practical scheme of exact dynamics calculation of any V-M system.

2. Universal properties of Vlasov-Maxwell system

It is well known that Vlasov equation (VE) $\widehat{L}f=0$ , where the operator is $\widehat{L}=\left[{\text{∂}}_t+\upsilon ·\text{∇}-L{F}_{\upsilon }·{\text{∂}}_p\right]$ and LF_υ  = e[E(r, t) + υ × B(r, t)] represents the Lorentz force and p(υ) = υΓ(υ) and $\mathrm{\Gamma }\left(\upsilon \right)=\frac{1}{\sqrt{1-\upsilon ·\upsilon }}$ correspond to the conservation law of total particle number ∫ fd ³ rd ³ υ = cons tan t. We should notice an important fact: if a beam is a particle number conserved system, its any portion or subsystem is not always also particle number conserved. For an unconserved subsystem, its distribution function g does not meet the VE. An unconserved subsystem means there are always realistic particles leaving it. A macroparticle is analogous to such an unconserved subsystem because it will break into equal pieces having different destinations at next time point t + Δt. Namely, there are always pieces leaving the macroparticle.

According to Maxwell equations (MEs), (E, B) couples with (M ₀, M ₁), where M_i  =  ∫ υⁱfd ³ υ, and formally decouples with all M _{i >  = 2}. According to open fluid equations set $\left\{\int {\upsilon }^{i>=0}\widehat{L}f{d}^3\upsilon =0,0=<i<\infty \right\}$ , each equation $\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0$ reflects a universal relation among ∂ _tM_i , ∇M _i + 1 and at least M _i − 1. This seems that exact solutions of all higher order moments M _{i >  = 2} are necessary for that of (E, B). On the other hand, because p(υ) and Γ(υ) are nonlinear functions of υ and hence − ∫ υⁱLF_υ  · ∂ _pfd ³ υ, a term in the equation $\int {\upsilon }^{i>=1}\widehat{L}f{d}^3\upsilon =0$ , is dependent on moments M_m from m = i to m = ∞, each equation $\int {\upsilon }^{i>=1}\widehat{L}f{d}^3\upsilon =0$ will imply a relation among moments M_m from m = i to m = ∞, and thus all equations $\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0$ from i = 1 to i = ∞ will mean a ∞ × ∞ matrix describing the relation among all moments M_m from m = 1 to m = ∞. Clearly, obtaining exact solutions of (E, B) from such an open equation set $\left\{\int {\upsilon }^{i>=0}\widehat{L}f{d}^3\upsilon =0,0=<i<\infty \right\}$ is impossible and, as discussed below, also unnecessary.

Another open set {D_i , 0 =  < i < ∞}, where $D_i={\left[\frac{M_i}{M_0}-{\left(\frac{M_1}{M_0}\right)}^i\right]}^{\ast }{M}_0$ , can be defined naturally through the M-set. Clearly, D ₀ = 0 and D ₁ = 0 automatically exist. Each equation $\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0$ can be expressed through the D-set

$$A_i{\partial }_t{D}_i+B_{i+1}\nabla D_{i+1}+\sum _{m>=i}‍C_m{D}_m=1,$$ E3

and coefficients A_i , B_i , C_i are known functionals of (E, B, M ₀, M ₁). Starting from the i = 1 case, we can formally obtain an expression of D ₂ in all terms D _{i >  = 3} and then substitute it into the i = 2 case and formally obtain an expression of D ₃ in all terms D _{i >  = 4}, etc. Finally, we will find that all D _{i >  = 2} are determined by D _∞ and all coefficients A_i , B_i , C_i . Namely, the open equation set $\left\{\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0,1=<i<\infty \right\}$ does not lead to a substantial constraint on (E, B, M ₀, M ₁).

According to MEs, (E, B) depends on (M ₀, M ₁) and is independent of the D-set. Therefore, the exact solutions of the D-set are not a necessary condition for those of (E, B). The open equation set $\left\{\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0,1=<i<\infty \right\}$ is only responsible for relations among those D_i and cannot have an effect on (E, B, M ₀, M ₁).

Two universal properties of any V-M system are favorable to achieve exact dynamics calculation. One is that the f of any V-M always has a thermal spread [18], or in stricter mathematical language, functions in a general form ɡ(r, t)∗δ(υ − u(r, t)), where ɡ and u are two functions of (r, t) and δ is Dirac function, cannot meet the VE or $\widehat{L}\left[ɡ\left(r,t\right)\ast \delta \left(\upsilon -u\left(r,t\right)\right)\right]\ne 0$ . The other is that the open set {M_i , 0 =  < i < ∞} has a closed subset {M ₀, M ₁}, which meets a closed equation set. Detailed proof of the second universal property is presented below.

Due to the first universal property, the whole system can be viewed as a superposition of a cold subsystem denoted by $F\left(f\right)=\delta \left(\upsilon -\frac{\int \mathit{f\upsilon }{d}^3\upsilon }{\int f{d}^3\upsilon }\right)\ast \int \left[f\ast \delta \left(\upsilon -\frac{\int \mathit{f\upsilon }{d}^3\upsilon }{\int f{d}^3\upsilon }\right)\right]d^3\upsilon$ , where δ is Dirac function, and a thermal subsystem denoted by f − F(f). The cold subsystem inherits a part of each moment ∫ υⁱf(r, υ, t)d ³ υ, so does the thermal subsystem. Namely, the cold subsystem has its own self-consistent field (E^F , B^F ), where ∇ · E^F ~ ∫ F(r, υ, t)d ³ υ and ∇ × B^F ~ ∫ υF(r, υ, t)d ³ υ. Likewise, the thermal subsystem has its own self-consistent field (E ^f − F , B ^f − F ), where E = E^F  + E ^f − F and B = B^F  + B ^f − F . (E^F , B^F ) plays a role of external fields to the thermal subsystem and so does (E ^f − F , B ^f − F ) to the cold subsystem. In particular, two subsystems have a common fluid velocity: $u=\frac{\int \mathit{\upsilon F}\left(r,\upsilon ,t\right)d^3\upsilon }{\int F\left(r,\upsilon ,t\right)d^3\upsilon }=\frac{\int \upsilon \left[f-F\right]\left(r,\upsilon ,t\right)d^3\upsilon }{\int \upsilon \left[f-F\right]\left(r,\upsilon ,t\right)d^3\upsilon }=\frac{\int \mathit{\upsilon f}\left(r,\upsilon ,t\right)d^3\upsilon }{\int \mathit{\upsilon f}\left(r,\upsilon ,t\right)d^3\upsilon }$ . Because F has a known dependence on υ, its higher order moments ∫ υ ^i > 1 Fd ³ υ are of simple forms u ^i > 1 ∗ ∫ [f∗δ(υ − u)]d ³ υ. The thermal subsystem inherits the thermal spread of f and hence all off-center moments, the dependence of f − F on υ is still unknown, and its higher-order moments ∫ υ ^i > 1[f − F]d ³ υ still obey a set of equations in infinite number. The equations of these higher-order moments ∫ υ ^i > 1 Fd ³ υ in infinite number are found, as proven later, to be equivalent to an equation of $u=\frac{\int \mathit{\upsilon F}\left(r,\upsilon ,t\right)d^3\upsilon }{\int F\left(r,\upsilon ,t\right)d^3\upsilon }$ .

Detailed strict proof has been made elsewhere [17, 18] (see the appendix of Ref. [17]).

Here, for convenience of readers, we briefly present the kernel of strict proof. For simplicity of symbols, we denote $\left[\widehat{L}+\left(\upsilon ·\text{∇}u\right){\text{∂}}_{\upsilon }\right]F$ as Ω and F as n ₀ δ(υ − u), where n ₀ ≡  ∫ [f∗δ(υ − u)]d ³ υ <  ∫ fd ³ υ = M ₀. From the definition of F, there is always $u=\frac{\int \mathit{\upsilon F}{d}^3\upsilon }{\int F{d}^3\upsilon }$ . Because of strict mathematical formulas ∫ d_υδ(υ − u)d ³ υ = 0, ∫ [υ − u]∗δ(υ − u)d ³ υ = 0, ∫ [υ − u]∗[υ · ∇ + (υ · ∇u)∂ _υ ]δ(υ − u)d ³ υ = 0, and ∫ [υ − u]∗W(r, υ, t)∗d_υδ(υ − u)d ³ υ =  ∫ [υ − u]∗W(r, υ, t)∣ _υ = u ∗d_υδ(υ − u)d ³ υ, where W(r, υ, t) is the arbitrary function, the definition of Ω can be directly written as

$$\begin{array}{l}\text{Ω}=\left[{\partial }_t+u·\nabla \right]n_0\ast \delta \left(\upsilon -u\right)-n_0\ast \left[{{\partial }_{t}u+\left(L{F}_{\upsilon }\ast {\partial }_p\upsilon \right)|}_{\upsilon =u}\right]\hfill \\ \ast {\delta }^{\prime }=\left\{{{\displaystyle \int }}^{\text{}}\text{Ω}{d}^3\upsilon \right\}\ast \delta -\left\{{{\displaystyle \int }}^{\text{}}\left[\upsilon -u\right]\ast \text{Ω}{d}^3\upsilon \right\}\ast {\delta }^{\prime },\hfill \end{array}$$ E4

which has a strict solution Ω = 0. Therefore, there is a theorem: For any V-M system, its F(f) meets Ω = 0. Namely, $\widehat{L}f=0$ means Ω = 0.

Due to 0 =  ∫ υ ⁰Ωd ³ υ, ∫ υ ¹Ωd ³ υ = 0 will be be equivalent to [17, 18]

$${\partial }_t\left[p\left(u\right)\right]+e\left[E+u\times B\right]=0.$$ E5

For equations in infinite number 0 =  ∫ υ ^{i >  = 2}Ωd ³ υ, the Dirac function dependence of F on υ makes all these equations indeed to be equivalent to a same equation, Eq. (5). This can be easily verified by simple algebra.

Few authors use Eq. (5) to calculate (E, B) in plasmas [17, 18, 19, 20, 21, 22]. As previously pointed out, the equation $\int \upsilon \widehat{L}f{d}^3\upsilon =0$ reflects a relation among all D _{i >  = 2}. To determine these D _{i >  = 2}, the whole set $\left\{\int {\upsilon }^{i>=0}\widehat{L}f{d}^3\upsilon =0,0=<i<\infty \right\}$ is required.

The whole system is described by

$$\widehat{L}f=0;$$ E6a

$${\partial }_{t}E=-\left[\nabla ·E-\mathit{Ze}{N}_i\right]\frac{\int \mathit{\upsilon f}{d}^3\upsilon }{\int f{d}^3\upsilon }+\nabla \times B;\nabla ·B=0$$ E6b

$$\nabla \times E=-{\partial }_{t}B.$$ E6c

Likewise, the cold subsystem by

$$\left[\widehat{L}+\left(\upsilon ·\nabla u\right){\partial }_{\upsilon }\right]F=0;$$ E7a

$${\partial }_{t}E=-\left[\nabla ·E-\mathit{Ze}{N}_i\right]\frac{\int \mathit{\upsilon F}{d}^3\upsilon }{\int F{d}^3\upsilon }+\nabla \times B;\nabla ·B=0$$ E7b

$$\nabla \times E=-{\partial }_{t}B,$$ E7c

and the thermal subsystem by

$$\left[\widehat{L}\left(f-F\right)-\left(\upsilon ·\nabla u\right){\partial }_{\upsilon }F\right]=0;$$ E8a

$${\partial }_{t}E=-\left[\nabla ·E-\mathit{Ze}{N}_i\right]\frac{\int \upsilon \left(f-F\right)d^3\upsilon }{\int \left(f-F\right)d^3\upsilon }+\nabla \times B;\nabla ·B=0$$ E8b

$$\nabla \times E=-{\partial }_{t}B.$$ E8c

Two universal properties imply a universal scheme of exact dynamics calculation of any V-M system. That is, the VE $\widehat{L}f=0$ implies Ω = 0, which implies {∫ υⁱ Ωd ³ υ = 0, 0 =  < i < ∞}, which implies a relation between $\frac{\int {\upsilon }^{1}F{d}^3\upsilon }{\int {\upsilon }^{0}F{d}^3\upsilon }$ and (E, B), which is also a relation between $\frac{\int {\upsilon }^{1}f{d}^3\upsilon }{\int {\upsilon }^{0}f{d}^3\upsilon }$ and (E, B). On the other hand, with the help of the CE $\int {\upsilon }^0\widehat{L}f{d}^3\upsilon =0$ , the open equation subset $\left\{\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0,1=<i<\infty \right\}$ will be equivalent to a ∞ × ∞ matrix equation which gives exact expressions of all higher order off-center moments D _{i >  = 2} in terms of (E, B, M ₀, M ₁). After obtaining exact solutions of (E, B, M ₀, M ₁) from Eq. (2) and MEs, all D _{i >  = 2} can be known through $\left\{\int {\upsilon }^i\widehat{L}f{d}^3\upsilon =0,1=<i<\infty \right\}$ .

This universal scheme gets rid of inevitable shortcomings of other methods. For example, well-known PIC simulation is to introduce macroparticle representation of initial distribution function and to calculate coupled evolution of these macroparticles and E and B. As analyzed above, such a macroparticle is unrealistic because it approximates that all contained realistic particles always stay in it and completely ignore the fact that realistic particles can leave it. Moreover, because the coupled evolution is calculated in an alternatively updating manner, this makes calculated E and B to be dependent of the grainness parameter which reflects how many realistic particles are contained in a macroparticle and hence are not exact solutions of E and B. Thus, even if f can be expressed through those macroparticles’ positions and velocities, inexact solutions of E and B make such an expression of f to be inexact.

Therefore, even though advanced computer group is applied to PIC simulation, its abovementioned inherent shortcoming cannot warrant more advanced computer condition corresponding to Exacter solutions of f (because E and B are inexact solutions). Namely, advanced computer condition does not get reasonable usage. In contrast, if the phase space dynamics is conducted according to above-described manner based on six-equation set, the power of the advanced computer condition can be fully utilized to get higher resolution of exact solutions of f (because E and B are exact solutions).

3. Universal properties of Newton-Maxwell system

Eq. (5) can also be derived from Newton-Maxwell (N-M) system, which contains N relativistic Newton equations (RNEs) and four MEs. Actually, for a group of electrons {r_i , d_tr_i }, we can always define two fields (in Lagrangian expression) [19]:

$$u\left(r_i\left(t\right),t\right)\equiv \frac{\sum _j‍d_t{r}_j\delta \left(r_i\left(t\right)-r_j\left(t\right)\right)}{\sum _j‍\delta \left(r_i\left(t\right)-r_j\left(t\right)\right)},$$ E9

$$\mathit{RV}\left(r_i\left(t\right),t\right)\equiv \frac{\sum _j‍\left[d_t{r}_j-d_t{r}_i\right]\delta \left(r_i\left(t\right)-r_j\left(t\right)\right)}{\sum _j‍\delta \left(r_i\left(t\right)-r_j\left(t\right)\right)},$$ E10

which are fluid velocity field and relative velocity field, respectively. Especially, there is always a formula

$$d_t{r}_i\equiv u\left(r_i\left(t\right),t\right)-\mathit{RV}\left(r_i\left(t\right),t\right),$$ E11

where ≡ means the formula is automatically valid for any r_i (t). If applying ∇ _{ri (t)} to this automatically valid formula, we will obtain another automatically valid formula:

$${\nabla }_{r_i\left(t\right)}\left[u\left(r_i\left(t\right),t\right)-\mathit{RV}\left(r_i\left(t\right),t\right)\right]\equiv {\nabla }_{r_i\left(t\right)}{d}_t{r}_i\equiv d_t{\nabla }_{r_i\left(t\right)}{r}_i\equiv d_{t}1\equiv 0.$$ E12

Because N electrons are described by four MEs and N RNEs, this automatically valid formula enables each RNE to be written as

$$\begin{array}{c}0=d_{t}p\left(d_t{r}_i\right)+eE\left(r_i\left(t\right)\text{,}t\right)+e{d}_t{r}_i\times B\left(r_i\left(t\right)\text{,}t\right)=d_{t}p\left(u\left(r_i\left(t\right)\text{,}t\right)-RV\left(r_i\left(t\right)\text{,}t\right)\right)\\ +eE\left(r_i\left(t\right)\text{,}t\right)+e\left[u\left(r_i\left(t\right)\text{,}t\right)-RV\left(r_i\left(t\right)\text{,}t\right)\right]\times B\left(r_i\left(t\right)\text{,}t\right)={\partial }_{t}p(u\left(r_i\left(t\right)\text{,}t\right)\\ -RV\left(r_i\left(t\right)\text{,}t\right))+d_t{r}_i\times {\nabla }_{r_i\left(t\right)}p\left[u\left(r_i\left(t\right)\text{,}t\right)-RV\left(r_i\left(t\right)\text{,}t\right)\right]+eE\left(r_i\left(t\right)\text{,}t\right)\\ +e\left[u\left(r_i\left(t\right)\text{,}t\right)-RV\left(r_i\left(t\right)\text{,}t\right)\right]\times B\left(r_i\left(t\right)\text{,}t\right)={\partial }_{t}p\left(u\left(r_i\left(t\right)\text{,}t\right)-RV\left(r_i\left(t\right)\text{,}t\right)\right)\\ +eE\left(r_i\left(t\right)\text{,}t\right)+e\left[u\left(r_i\left(t\right)\text{,}t\right)-RV\left(r_i\left(t\right)\text{,}t\right)\right]\times B\left(r_i\left(t\right)\text{,}t\right),\end{array}$$ E13

where $p\left(\mathit{var}\right)=\frac{\mathit{var}}{\sqrt{1-\mathit{\text{va}}\mathit{\text{r}}\ast \mathit{var}}}$ is relativistic momentum, because no matter what the relative field is, the RNE of every electron is always valid. Thus, the condition for the RNE being valid under arbitrary value of the RV field, of course including a common value RV = 0, is that the following equation is valid:

$$0={\partial }_{t}p\left(u\left(r_i\left(t\right),t\right)\right)+\mathit{eE}\left(r_i\left(t\right),t\right)+e\left[u\left(r_i\left(t\right),t\right)\right]\times B\left(r_i\left(t\right),t\right).$$ E14

Its Euler expression reads

$$0={\partial }_{t}p\left(u\left(r,t\right)\right)+\mathit{eE}\left(r,t\right)+e\left[u\left(r,t\right)\right]\times B\left(r,t\right),$$ E15

which is just Eq. (5). This motion equation of u(r, t) and four MEs form a closed five-equation set. No matter how large N is, it merely means more RV-field fluid elements being calculated. Finer description of the RV field will lead to Exacter information on thermal distribution, but the RV-field has no contribution to the electric current and hence E and B. Thus, the exact solutions of the N + 4 equation set (N RNEs and four MEs) can be obtained from an N + 5 equation set (N equations of RV(r_i (t), t) and a closed five-equation set of (u, E, B). Examples of electron density profile based on this closed five-equation set have been presented elsewhere [17, 18, 19, 20].

4. Examples of exact solution of phase space distribution

From Eq. (5) and MEs, we can find that W ≡ ∂ _tp(u) + ZN_ir/3 and u × B meeting a same equation [17, 18]. Further strict analysis reveals that this will lead to a self-consistent equation of u [17, 18]:

$${\partial }_{\mathit{tt}}p+\left[\left({\partial }_t\nabla ·p\right)\frac{p}{\sqrt{1+p^2}}\right]+\nabla \times \nabla \times p=-\left[\mathit{uZ}{N}_i\right]+\mathit{POT}$$ E16

where Z is ionic charge, N_i is ionic density, and POT is a constant vector determined by initial conditions. After solving u(r, t) from this equation, B(r, t) will be known from a formula ∂ _tp + ZN_ir/3 =  − λ[u × B] [17, 18], where λ ≠ 1 is a constant, and then E(r, t) from E = ∂ _tp + u × B [17, 18].

Here, we give some typical exact solutions of the distribution of f over the υ-space [17, 19, 20, 21]. If f is an exact solution of the VE, any monovariable function of f, or ɡ(f), is also an exact solution. Because Eq. (5) and MEs can have exact solutions of (E, B) in running-wave form, $E=E\left(r-\frac{1}{\eta }\mathit{ct}\right)$ and $B=B\left(r-\frac{1}{\eta }\mathit{ct}\right)$ , we should note a relation between E and B, $E=-\frac{1}{\eta }c\times B+\text{∇}\mathrm{\Phi }\left(r-\frac{1}{\eta }\mathit{ct}\right)$ , which arises from ∇ × E =  − ∂ _tB. Here, $\mathrm{\Phi }\left(r-\frac{1}{\eta }\mathit{ct}\right)$ is a scalar function but cannot be simply taken as electrostatic (ES) potential (because $-\frac{1}{\eta }c\times B$ also has divergence or $\text{∇}·\left(-\frac{1}{\eta }c\times B\right)=\frac{1}{\eta }c·\text{∇}\times B\ne 0$ ). In this case, the VE can be further written as

$$\begin{array}{c}0=\left[{\partial }_t+\upsilon ·\nabla -\left[E+\upsilon \times B\right]·{\partial }_{p\left(\upsilon \right)}\right]f=\left[{\partial }_t-E·{\partial }_{p\left(\upsilon \right)}\right]f+\upsilon ·\left[\nabla -B\times {\partial }_{p\left(\upsilon \right)}\right]\\ f=\left[{\partial }_t+\frac{1}{\eta }c\times B·{\partial }_{p\left(\upsilon \right)}\right]f+\upsilon ·\left[\nabla -B\times {\partial }_{p\left(\upsilon \right)}\right]f-\nabla \mathrm{\Phi }·{\partial }_{p\left(\upsilon \right)}\\ f=\left[{\partial }_t+\frac{1}{\eta }c·B\times {\partial }_{p\left(\upsilon \right)}\right]f+\upsilon ·\left[\nabla -B\times {\partial }_{p\left(\upsilon \right)}\right]\\ f-\nabla \mathrm{\Phi }·{\partial }_{p\left(\upsilon \right)}f=\left(\upsilon -\frac{1}{\eta }c\right)·\left[\nabla -B\times {\partial }_{p\left(\upsilon \right)}\right]f-\nabla \mathrm{\Phi }·{\partial }_{p\left(\upsilon \right)}f.\end{array}$$ E17

It is easy to verify that any function of $p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}$ will meet

$$\begin{array}{l}0=\left[{\partial }_t-E·{\partial }_{p\left(\upsilon \right)}\right]ɡ\left(p+{{\displaystyle \int }}^{\text{}}E\left(r-\frac{1}{\eta }ct\right)dt\right)\hfill \\ =\left[-\frac{1}{\eta }c·\nabla +\frac{1}{\eta }c·B\times {\partial }_{p\left(\upsilon \right)}\right]ɡ=-\frac{1}{\eta }c·\left[\nabla -B\times {\partial }_{p\left(\upsilon \right)}\right]ɡ,\hfill \end{array}$$ E18

where we have used the property ${\text{∂}}_t\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}=-\frac{1}{\eta }c·\text{∇}\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}$ . Thus, if ∇Φ ≡ 0, any mono-variable function of $p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}$ , or $ɡ\left(p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}\right)$ , will be an exact solution of the VE.

On the other hand, for more general ∇Φ, we can find that any monovariable function of $\sqrt{1+\frac{p^2}{c^2}}-\frac{1}{\eta }c·p+\mathrm{\Phi }$ or $ɡ\left(\sqrt{1+\frac{p^2}{c^2}}-\frac{1}{\eta }c·p+\mathrm{\Phi }\right)$ , is an exact solution of the VE. According to Eq. (14), ${\text{∂}}_p\left[\sqrt{1+\frac{p^2}{c^2}}-\frac{1}{\eta }c·p\right]$ will contribute a vector parallel to $\left(\upsilon -\frac{1}{\eta }c\right)$ and hence make the operator $\left(\upsilon -\frac{1}{\eta }c\right)·B\times {\text{∂}}_{p\left(\upsilon \right)}$ has zero contribution.

Therefore, for running-wave form, $E=E\left(r-\frac{1}{\eta }\mathit{ct}\right)$ and $B=B\left(r-\frac{1}{\eta }\mathit{ct}\right)$ , the phase space distribution, if ∇Φ ≡ 0, can be described by a positive-valued function of $p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}$ , for example, $exp\left[-{\left(p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}\right)}^2\right]$ , ${sin}^2\left(exp\left[-{\left(p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}\right)}^2\right]\right)$ , etc. We can further pick out reasonable solutions according to the running-wave form of u:

$$u=\frac{\int \frac{p}{\sqrt{1+p^2}}ɡ\left(p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}\right)d^{3}p}{\int ɡ\left(p+\int E\left(r-\frac{1}{\eta }\mathit{ct}\right)\mathit{dt}\right)d^{3}p}.$$ E19

Likewise, same procedure exists for more general ∇Φ and $ɡ\left(\sqrt{1+\frac{p^2}{c^2}}-\frac{1}{\eta }c·p+\mathrm{\Phi }\right)$ .

Actually, a set of macroscopic functions (E, B, u) can have multiple microscopic solutions of the VE. Therefore, usually we know the distribution of f over the υ-space from the initial condition f(r, p, t = 0). From the function dependence of f(r, p, t = 0) on p, we can obtain the function dependence of ɡ on $\sqrt{1+\frac{p^2}{c^2}}-\frac{1}{\eta }c·p$ and hence determine detailed function form of ɡ.

Detailed procedure of determining the function form of ɡ is described as follows [9, 11]: we can seek for special space position R in which $E\left(R,0\right)=-\frac{1}{\eta }c\times B\left(R,0\right),$ or ∇Φ(r, 0)| _r = R  = 0, exists. The initial distribution at R, i.e., f(R, p, t = 0), is thus a monovariable function p. At the same time, two expressions are equivalent, and hence there is f(R, p, t = 0) = ɡ(K + Φ(R, 0)) = ɡ(K), where $K=\sqrt{1+\frac{p^2}{c^2}}-\frac{1}{\eta }c·p$ and Φ(R, 0) = 0 (if ∇Φ(r, 0)| _r = R  = 0). Because of certain relations between p and K, once the expansion coefficient c_i in f(R, p, t = 0) = ∑ _i  ‍ c_ipⁱ is known, the expansion coefficient d_i in ɡ(K) = ∑ _i  ‍ d_iKⁱ are also easy to be calculated.

Clearly, BGK modes [23] are analytic exact solutions of the VE in B ≡ 0 case:

$$0=\left[{\partial }_t-E\left(r-\frac{1}{\eta }\mathit{ct}\right)·{\partial }_{p\left(\upsilon \right)}\right]f+\upsilon ·\nabla f,$$ E20

whose solutions are monovariable functions of $\varphi \left(r-\frac{1}{\eta }\mathit{ct}\right)+\sqrt{1+p^2}-\frac{1}{\eta }c·p$ , where ϕ is scalar potential and E =  − ∇ϕ. Moreover, there is a similar procedure of determining function form of ɡ.

We should note that K is a nonlinear function of υ and the maximum value of K or K _max is reached at $\upsilon =\frac{1}{\eta }c$ . Thus, even if ɡ is a Dirac function of K + Φ, ɡ cannot be a Dirac function of υ, i.e., ɡ~δ(υ − u(r, t)). The nonlinear function relation between K and υ determines that ɡ is at least a summation of two Dirac components: ɡ = f ₁(r, t)δ(υ − u ₁(r, t)) + f ₂(r, t)δ(υ − u ₂(r, t)) + …. This agrees with previous conclusion that functions of a general form F ₁(r, t)δ(υ − F ₂(r, t)) cannot meet the VE.

We should also note that, because of nonlinear function relation between ɡ and K + Φ, the maximum of ɡ, or ɡ _max, is usually reached at K + Φ ≠ K _max + Φ. Namely, if ɡ _max is reached at K = K _ɡmax, this K _ɡmax usually corresponds to two values of υ. In contrast, K _max merely corresponds to a value of υ. Thus, the contour plot of ɡ in the phase space will take on complicated structures, such as hole, island, etc. [17, 19].

5. Summary

By analyzing the bottleneck of numerical study the 6D phase space dynamics of a V-M system, we present, based on two universal properties, a universal and practical scheme of exact dynamics calculation of any V-M system. Our scheme is to harness, to the largest extent, available computer group condition to obtain as exact as possible solution. It gets rid of some inherent factors (of other methods), which yield incorrect physical results.