## 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^{7} 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^{7} 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*and a charge

_{e}*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 6*l* functions: {*x _{i}
*(

*t*),

*y*(

_{i}*t*),

*z*(

_{i}*t*),

*d*(

_{t}x_{i}*t*),

*d*(

_{t}y_{i}*t*),

*d*(

_{t}z_{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*,

_{t}x*y*,

*z*

_{2k }(

*t*= 0) −

*d*,

_{t}x*y*,

*z*

_{2k + 1}(

*t*= 0), is,

*d*,

_{t}x*y*,

*z*

_{2k }−

*d*,

_{t}x*y*,

*z*

_{2k + 1}is always taken as 0 or the 2

*k*th realistic particle is always bound with the 2

*k*+ 1th one (their relative motion is ignored completely). This is a distorted description. The 2

*k*th realistic particle and the 2

*k*+ 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*,

_{t}x*y*,

*z*

_{2k }(

*t*= 0) and

*d*,

_{t}x*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*and a fixed charge

_{e}*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

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

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 ∂*
_{t}E* 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)
*LF _{υ}
* =

*e*[

*E*(

*r*,

*t*) +

*υ*×

*B*(

*r*,

*t*)] represents the Lorentz force and

*p*(

*υ*) =

*υ*Γ(

*υ*) and

*fd*

^{3}

*rd*

^{3}

*υ*=

*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*
_{0}, *M*
_{1}), where *M _{i}
* = ∫

*υ*

^{i}fd^{3}

*υ*, and formally decouples with all

*M*

_{ i > = 2}. According to open fluid equations set

*, ∇*

_{t}M_{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 − ∫

*υ*· ∂

^{i}LF_{υ}

_{p}fd^{3}

*υ*, a term in the equation

*M*from

_{m}*m*=

*i*to

*m*= ∞, each equation

*M*from

_{m}*m*=

*i*to

*m*= ∞, and thus all equations

*i*= 1 to

*i*= ∞ will mean a ∞ × ∞ matrix describing the relation among all moments

*M*from

_{m}*m*= 1 to

*m*= ∞. Clearly, obtaining exact solutions of (

*E*,

*B*) from such an open equation set

Another open set {*D _{i}
*, 0 = <

*i*< ∞}, where

*M*-set. Clearly,

*D*

_{0}= 0 and

*D*

_{1}= 0 automatically exist. Each equation

*D*-set

and coefficients *A _{i}
*,

*B*,

_{i}*C*are known functionals of (

_{i}*E*,

*B*,

*M*

_{0},

*M*

_{1}). Starting from the

*i*= 1 case, we can formally obtain an expression of

*D*

_{2}in all terms

*D*

_{ i > = 3}and then substitute it into the

*i*= 2 case and formally obtain an expression of

*D*

_{3}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*. Namely, the open equation set

_{i}*E*,

*B*,

*M*

_{0},

*M*

_{1}).

According to MEs, (*E*, *B*) depends on (*M*
_{0}, *M*
_{1}) 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
*D _{i}
* and cannot have an effect on (

*E*,

*B*,

*M*

_{0},

*M*

_{1}).

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
*M _{i}
*, 0 = <

*i*< ∞} has a closed subset {

*M*

_{0},

*M*

_{1}}, 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
*δ* is Dirac function, and a thermal subsystem denoted by *f* − *F*(*f*). The cold subsystem inherits a part of each moment ∫ *υ ^{i}f*(

*r*,

*υ*,

*t*)

*d*

^{3}

*υ*, so does the thermal subsystem. Namely, the cold subsystem has its own self-consistent field (

*E*,

^{F}*B*), where ∇ ·

^{F}*E*~ ∫

^{F}*F*(

*r*,

*υ*,

*t*)

*d*

^{3}

*υ*and ∇ ×

*B*~ ∫

^{F}*υF*(

*r*,

*υ*,

*t*)

*d*

^{3}

*υ*. 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*) plays a role of external fields to the thermal subsystem and so does (

^{F}*E*

^{ f − F },

*B*

^{ f − F }) to the cold subsystem. In particular, two subsystems have a common fluid velocity:

*F*has a known dependence on

*υ*, its higher order moments ∫

*υ*

^{ i > 1}

*Fd*

^{3}

*υ*are of simple forms

*u*

^{ i > 1}∗ ∫ [

*f*∗

*δ*(

*υ*−

*u*)]

*d*

^{3}

*υ*. 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*

^{3}

*υ*still obey a set of equations in infinite number. The equations of these higher-order moments ∫

*υ*

^{ i > 1}

*Fd*

^{3}

*υ*in infinite number are found, as proven later, to be equivalent to an equation of

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
*F* as *n*
_{0}
*δ*(*υ* − *u*), where *n*
_{0} ≡ ∫ [*f*∗*δ*(*υ* − *u*)]*d*
^{3}
*υ* < ∫ *fd*
^{3}
*υ* = *M*
_{0}. From the definition of *F*, there is always
*d _{υ}δ*(

*υ*−

*u*)

*d*

^{3}

*υ*= 0, ∫ [

*υ*−

*u*]∗

*δ*(

*υ*−

*u*)

*d*

^{3}

*υ*= 0, ∫ [

*υ*−

*u*]∗[

*υ*· ∇ + (

*υ*· ∇

*u*)∂

*]*

_{υ}*δ*(

*υ*−

*u*)

*d*

^{3}

*υ*= 0, and ∫ [

*υ*−

*u*]∗

*W*(

*r*,

*υ*,

*t*)∗

*d*(

_{υ}δ*υ*−

*u*)

*d*

^{3}

*υ*= ∫ [

*υ*−

*u*]∗

*W*(

*r*,

*υ*,

*t*)∣

_{ υ = u }∗

*d*(

_{υ}δ*υ*−

*u*)

*d*

^{3}

*υ*, where

*W*(

*r*,

*υ*,

*t*) is the arbitrary function, the definition of Ω can be directly written as

which has a strict solution Ω = 0. Therefore, there is a theorem: For any V-M system, its *F*(*f*) meets Ω = 0. Namely,

Due to 0 = ∫ *υ*
^{0}Ω*d*
^{3}
*υ*, ∫ *υ*
^{1}Ω*d*
^{3}
*υ* = 0 will be be equivalent to [17, 18]

For equations in infinite number 0 = ∫ *υ*
^{
i > = 2}Ω*d*
^{3}
*υ*, 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
*D*
_{
i > = 2}. To determine these *D*
_{
i > = 2}, the whole set

The whole system is described by

Likewise, the cold subsystem by

and the thermal subsystem by

Two universal properties imply a universal scheme of exact dynamics calculation of any V-M system. That is, the VE
*υ ^{i}
*Ω

*d*

^{3}

*υ*= 0, 0 = <

*i*< ∞}, which implies a relation between

*E*,

*B*), which is also a relation between

*E*,

*B*). On the other hand, with the help of the CE

*D*

_{ i > = 2}in terms of (

*E*,

*B*,

*M*

_{0},

*M*

_{1}). After obtaining exact solutions of (

*E*,

*B*,

*M*

_{0},

*M*

_{1}) from Eq. (2) and MEs, all

*D*

_{ i > = 2}can be known through

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*}, we can always define two fields (in Lagrangian expression) [19]:

_{t}r_{i}which are fluid velocity field and relative velocity field, respectively. Especially, there is always a formula

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:

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

(13) |

where
*RV* field, of course including a common value *RV* = 0, is that the following equation is valid:

Its Euler expression reads

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* ≡ ∂*
_{t}p*(

*u*) +

*ZN*/3 and

_{i}r*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]:

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 ∂

*+*

_{t}p*ZN*/3 = −

_{i}r*λ*[

*u*×

*B*] [17, 18], where

*λ*≠ 1 is a constant, and then

*E*(

*r*,

*t*) from

*E*= ∂

*+*

_{t}p*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* and *B*,
*E* = − ∂*
_{t}B*. Here,

(17) |

It is easy to verify that any function of

where we have used the property

On the other hand, for more general ∇Φ, we can find that any monovariable function of

Therefore, for running-wave form,
*u*:

Likewise, same procedure exists for more general ∇Φ and

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
*ɡ*.

Detailed procedure of determining the function form of *ɡ* is described as follows [9, 11]: we can seek for special space position *R* in which
*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
*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*is known, the expansion coefficient

_{i}p^{i}*d*in

_{i}*ɡ*(

*K*) = ∑

**

_{i}*d*are also easy to be calculated.

_{i}K^{i}Clearly, BGK modes [23] are analytic exact solutions of the VE in *B* ≡ 0 case:

whose solutions are monovariable functions of
*ϕ* 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
*ɡ* 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*
_{1}(*r*, *t*)*δ*(*υ* − *u*
_{1}(*r*, *t*)) + *f*
_{2}(*r*, *t*)*δ*(*υ* − *u*
_{2}(*r*, *t*)) + …. This agrees with previous conclusion that functions of a general form *F*
_{1}(*r*, *t*)*δ*(*υ* − *F*
_{2}(*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.