Open access peer-reviewed chapter

Phase Space Dynamics of Relativistic Particles

By Hai Lin

Submitted: April 17th 2017Reviewed: September 13th 2017Published: December 20th 2017

DOI: 10.5772/intechopen.70957

Downloaded: 1129


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.


  • Vlasov-Maxwell system
  • electron beam
  • phase space dynamics
  • microscopic distribution function
  • self-consistent fields

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 meand 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 pClevel, the number of realistic particles it contains will be about 10−12 + 19/1.6~6.3 × 107 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 nslevel. 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 × 107 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 meand a charge e) into a macroparticle (of a mass Nmeand 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 lrealistic particles in a beam, describing the beam will need 6lfunctions: {xi(t), yi(t), zi(t), dtxi(t), dtyi(t), dtzi(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 dtx, y, z2k(t = 0) − dtx, y, z2k + 1(t = 0), is, dtx, y, z2k − dtx, y, z2k + 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, z2k(t = 0) = x, y, z2k + 1(t = 0), tend to move to different destinations at the next time point 0 + Δt(where Δtis the time step), because of their different initial velocities dtx, y, z2k(t = 0) and dtx, y, z2k + 1(t = 0). Namely, actual picture should be described as follows: at every time point t, a macroparticle will break into Nequal pieces having their own destinations at t + Δt, rather than globally move into a destination at t + Δtor being equivalently expressed as: each macroparticle at tis a gathering of Nequal 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 Nmeand a fixed charge Ne.

The direct result of the rigid-macroparticle approximation is that the charge density profile nof the beam at next time point t + Δtis 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 Δtto 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 Eand 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 Ecan be known from all particles’ information at t, Bcannot. Due to the term ∂tEin Eq. (1), the information of Bat thave dependence on all particles’ information at t + Δt. This implies that the NE is a complicated differential equation of ri(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., {ri(t + Δt); 1 =  < i =  < N}. Thus, Ndifference versions form a linear equation set of {ri(t + Δt); 1 =  < i =  < N} whose solutions involve inevitably a N × Nmatrix. The larger Nis, 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) L̂f=0, where the operator is L̂=t+υ·LFυ·pand LFυ = e[E(r, t) + υ × B(r, t)] represents the Lorentz force and p(υ) = υΓ(υ) and Γυ=11υ·υcorrespond to the conservation law of total particle number ∫ fd3 rd3 υ = 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 gdoes 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 (M0, M1), where Mi =  ∫ υifd3 υ, and formally decouples with all Mi >  = 2. According to open fluid equations set υi>=0L̂fd3υ=00=<i<, each equation υiL̂fd3υ=0reflects a universal relation among ∂tMi, ∇Mi + 1 and at least Mi − 1. This seems that exact solutions of all higher order moments Mi >  = 2 are necessary for that of (E, B). On the other hand, because p(υ) and Γ(υ) are nonlinear functions of υand hence − ∫ υiLFυ · ∂pfd3 υ, a term in the equation υi>=1L̂fd3υ=0, is dependent on moments Mmfrom m = ito m = ∞, each equation υi>=1L̂fd3υ=0will imply a relation among moments Mmfrom m = ito m = ∞, and thus all equations υiL̂fd3υ=0from i = 1 to i = ∞ will mean a ∞ × ∞ matrix describing the relation among all moments Mmfrom m = 1 to m = ∞. Clearly, obtaining exact solutions of (E, B) from such an open equation set υi>=0L̂fd3υ=00=<i<is impossible and, as discussed below, also unnecessary.

Another open set {Di, 0 =  < i < ∞}, where Di=MiM0M1M0iM0, can be defined naturally through the M-set. Clearly, D0 = 0 and D1 = 0 automatically exist. Each equation υiL̂fd3υ=0can be expressed through the D-set


and coefficients Ai, Bi, Ciare known functionals of (E, B, M0, M1). Starting from the i = 1 case, we can formally obtain an expression of D2 in all terms Di >  = 3 and then substitute it into the i = 2 case and formally obtain an expression of D3 in all terms Di >  = 4, etc. Finally, we will find that all Di >  = 2 are determined by D and all coefficients Ai, Bi, Ci. Namely, the open equation set υiL̂fd3υ=01=<i<does not lead to a substantial constraint on (E, B, M0, M1).

According to MEs, (E, B) depends on (M0, M1) 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 υiL̂fd3υ=01=<i<is only responsible for relations among those Diand cannot have an effect on (E, B, M0, M1).

Two universal properties of any V-M system are favorable to achieve exact dynamics calculation. One is that the fof 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 uare two functions of (r, t) and δis Dirac function, cannot meet the VE or L̂ɡrtδυurt0. The other is that the open set {Mi, 0 =  < i < ∞} has a closed subset {M0, M1}, 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 Ff=δυd3υfd3υfδυd3υfd3υd3υ, where δis Dirac function, and a thermal subsystem denoted by f − F(f). The cold subsystem inherits a part of each moment ∫ υif(r, υ, t)d3 υ, so does the thermal subsystem. Namely, the cold subsystem has its own self-consistent field (EF, BF), where ∇ · EF~ ∫ F(r, υ, t)d3 υand ∇ × BF~ ∫ υF(r, υ, t)d3 υ. Likewise, the thermal subsystem has its own self-consistent field (Ef − F, Bf − F), where E = EF + Ef − Fand B = BF + Bf − F. (EF, BF) plays a role of external fields to the thermal subsystem and so does (Ef − F, Bf − F) to the cold subsystem. In particular, two subsystems have a common fluid velocity: u=υFrυtd3υFrυtd3υ=υfFrυtd3υυfFrυtd3υ=υfrυtd3υυfrυtd3υ. Because Fhas a known dependence on υ, its higher order moments ∫ υi > 1 Fd3 υare of simple forms ui > 1 ∗ ∫ [fδ(υ − u)]d3 υ. The thermal subsystem inherits the thermal spread of fand hence all off-center moments, the dependence of f − Fon υis still unknown, and its higher-order moments ∫ υi > 1[f − F]d3 υstill obey a set of equations in infinite number. The equations of these higher-order moments ∫ υi > 1 Fd3 υin infinite number are found, as proven later, to be equivalent to an equation of u=υFrυtd3υFrυtd3υ.

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 L̂+υ·uυFas Ω and Fas n0 δ(υ − u), where n0 ≡  ∫ [fδ(υ − u)]d3 υ <  ∫ fd3 υ = M0. From the definition of F, there is always u=υFd3υFd3υ. Because of strict mathematical formulas ∫ dυδ(υ − u)d3 υ = 0, ∫ [υ − u]∗δ(υ − u)d3 υ = 0, ∫ [υ − u]∗[υ · ∇ + (υ · ∇u)∂υ]δ(υ − u)d3 υ = 0, and ∫ [υ − u]∗W(r, υ, t)∗dυδ(υ − u)d3 υ =  ∫ [υ − u]∗W(r, υ, t)∣ υ = udυδ(υ − u)d3 υ, 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, L̂f=0means Ω = 0.

Due to 0 =  ∫ υ0Ωd3 υ, ∫ υ1Ωd3 υ = 0 will be be equivalent to [17, 18]


For equations in infinite number 0 =  ∫ υi >  = 2Ωd3 υ, the Dirac function dependence of Fon υ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 υL̂fd3υ=0reflects a relation among all Di >  = 2. To determine these Di >  = 2, the whole set υi>=0L̂fd3υ=00=<i<is required.

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 L̂f=0implies Ω = 0, which implies {∫ υiΩd3 υ = 0, 0 =  < i < ∞}, which implies a relation between υ1Fd3υυ0Fd3υand (E, B), which is also a relation between υ1fd3υυ0fd3υand (E, B). On the other hand, with the help of the CE υ0L̂fd3υ=0, the open equation subset υiL̂fd3υ=01=<i<will be equivalent to a ∞ × ∞ matrix equation which gives exact expressions of all higher order off-center moments Di >  = 2 in terms of (E, B, M0, M1). After obtaining exact solutions of (E, B, M0, M1) from Eq. (2) and MEs, all Di >  = 2 can be known through υiL̂fd3υ=01=<i<.

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 Eand 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 Eand Bto be dependent of the grainness parameter which reflects how many realistic particles are contained in a macroparticle and hence are not exact solutions of Eand B. Thus, even if fcan be expressed through those macroparticles’ positions and velocities, inexact solutions of Eand Bmake such an expression of fto 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 Eand Bare 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 Eand Bare exact solutions).


3. Universal properties of Newton-Maxwell system

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


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 ri(t). If applying ∇ ri(t) to this automatically valid formula, we will obtain another automatically valid formula:


Because Nelectrons are described by four MEs and NRNEs, this automatically valid formula enables each RNE to be written as


where pvar=var1varvaris 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 RVfield, 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 Nis, it merely means more RV-field fluid elements being calculated. Finer description of the RVfield will lead to Exacter information on thermal distribution, but the RV-field has no contribution to the electric current and hence Eand B. Thus, the exact solutions of the N + 4 equation set (NRNEs and four MEs) can be obtained from an N + 5 equation set (Nequations of RV(ri(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) + ZNir/3 and u × Bmeeting a same equation [17, 18]. Further strict analysis reveals that this will lead to a self-consistent equation of u[17, 18]:


where Zis ionic charge, Niis ionic density, and POTis 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 + ZNir/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 fover the υ-space [17, 19, 20, 21]. If fis 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=Er1ηctand B=Br1ηct, we should note a relation between Eand B, E=1ηc×B+Φr1ηct, which arises from ∇ × E =  − ∂tB. Here, Φr1ηctis a scalar function but cannot be simply taken as electrostatic (ES) potential (because 1ηc×Balso has divergence or ·1ηc×B=1ηc·×B0). In this case, the VE can be further written as


It is easy to verify that any function of p+Er1ηctdtwill meet


where we have used the property tEr1ηctdt=1ηc·Er1ηctdt. Thus, if ∇Φ ≡ 0, any mono-variable function of p+Er1ηctdt, or ɡp+Er1ηctdt, will be an exact solution of the VE.

On the other hand, for more general ∇Φ, we can find that any monovariable function of 1+p2c21ηc·p+Φor ɡ1+p2c21ηc·p+Φ, is an exact solution of the VE. According to Eq. (14), p1+p2c21ηc·pwill contribute a vector parallel to υ1ηcand hence make the operator υ1ηc·B×pυhas zero contribution.

Therefore, for running-wave form, E=Er1ηctand B=Br1ηct, the phase space distribution, if ∇Φ ≡ 0, can be described by a positive-valued function of p+Er1ηctdt, for example, expp+Er1ηctdt2, sin2expp+Er1ηctdt2, etc. We can further pick out reasonable solutions according to the running-wave form of u:


Likewise, same procedure exists for more general ∇Φ and ɡ1+p2c21ηc·p+Φ.

Actually, a set of macroscopic functions (E, B, u) can have multiple microscopic solutions of the VE. Therefore, usually we know the distribution of fover 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 1+p2c21ηc·pand 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 Rin which ER0=1ηc×BR0,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=1+p2c21ηc·pand Φ(R, 0) = 0 (if ∇Φ(r, 0)| r = R = 0). Because of certain relations between pand K, once the expansion coefficient ciin f(R, p, t = 0) = ∑i ‍ cipiis known, the expansion coefficient diin ɡ(K) = ∑i ‍ diKiare also easy to be calculated.

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


whose solutions are monovariable functions of ϕr1ηct+1+p21ηc·p, where ϕis scalar potential and E =  − ∇ϕ. Moreover, there is a similar procedure of determining function form of ɡ.

We should note that Kis a nonlinear function of υand the maximum value of Kor Kmax is reached at υ=1η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 Kand υdetermines that ɡis at least a summation of two Dirac components: ɡ = f1(r, t)δ(υ − u1(r, t)) + f2(r, t)δ(υ − u2(r, t)) + …. This agrees with previous conclusion that functions of a general form F1(r, t)δ(υ − F2(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 + Φ ≠ Kmax + Φ. Namely, if ɡmax is reached at K = Kɡmax, this Kɡmax usually corresponds to two values of υ. In contrast, Kmax 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.

© 2017 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Hai Lin (December 20th 2017). Phase Space Dynamics of Relativistic Particles, Accelerator Physics - Radiation Safety and Applications, Ishaq Ahmad and Maaza Malek, IntechOpen, DOI: 10.5772/intechopen.70957. Available from:

chapter statistics

1129total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Motion of Electrons in Planar Ideal Undulator

By Nikolay Smolyakov

Related Book

First chapter

Introductory Chapter: Introduction to Ion Implantation

By Ishaq Ahmad and Waheed Akram

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us