Improved Algorithms for Calculating the Space-Charge Field in Vacuum Devices

Abstract

The space-charge field (SCF) is a key factor in vacuum electronic devices, accelerators, free electron lasers and plasma systems, etc. The calculation of the SCF is very important since it has a great influence on the precision of numerical simulation results. However, calculating the SCF usually takes a lot of time, especially when the number of simulated particles is large. In this paper, we used a vectorization, parallelization and truncation method to optimize the calculation of the SCF based on the traditional calculation algorithms. To verify the validity of the optimized SCF calculation algorithm, it was applied in the performance simulation of a millimeter wave traveling wave tube. The results showed that the time cost was reduced by three orders compared with conventional treatment. The proposed algorithm also has great potential applications in free electron lasers, accelerators and plasma systems.

1. Introduction

The space-charge field (SCF) [1,2,3,4,5,6] is an important physical factor in electronic vacuum devices, particle accelerators [7,8,9,10,11,12], free electron lasers (FEL) [13] and plasma systems [14,15,16,17,18]. It is generated by electrostatic interaction between charged particles. In these systems, the SCF directly changes the distribution of the charged particles. Thus, understanding the SCF can improve the accuracy of beam–wave interaction simulations. Usually, the electrostatic interaction between any two macro particles needs to be calculated in the calculation of the SCF. This will cause the simulation to take lot of time as the number of macro particles increases. Therefore, it is necessary to optimize the calculation method to calculate the SCF so that calculation efficiency can be improved.

In this paper, we proposed an optimized SCF calculation algorithm for an E-band traveling wave tube (TWTs). A Lagrangian [19,20] description was applied to simulate the operation of TWTs since it is widely used for calculating the beam–wave interaction process. In this description, the electron beam is treated as a series of macroparticles that can be tracked. Each macroparticle is distinguished by the phase when it enters the interaction region. The electromagnetic field is excited by the current determined by the position of the macro particles. The SCF generated by the electrostatic interaction of the macro particles is described by the disc model proposed by Tien [21,22]. A group of nonlinear equations that describe the physical process of the electromagnetic wave excitation and the motion of the macro particles have been established and solved by numerical methods. The phase-space evolution of the macro particles, the output power and bunching factor can be obtained via these simulations.

During the process of solving these equations, the SCF calculation takes the most time. In addition, in the calculation of some characteristics of vacuum devices such as intermodulation and high harmonics, the number of macro particles will increase to several thousands and even more. This will unfortunately cause a long simulation time. Hence, the reduction in the SCF calculation time is an important issue in the device simulation. Based on the traditional algorithm for SCF calculation, we used the vectorization method, the parallelization method and the truncation method to optimize the SCF calculation so that time is saved. The results show that the time cost was reduced by three orders, which demonstrates a very good effectiveness.

This paper is organized as follows: In Section 2, we present the expression of the SCF under the Discs Model and its normalized expression. In addition, an analysis of the macro particle distribution is presented, which is important for the numerical calculation of the SCF. In Section 3, a traditional algorithm and three optimized algorithms are applied in calculating the SCF. The traditional algorithm uses a simple loop statement. The optimized algorithms include three different optimization approaches: the parallel algorithm, the vectorization algorithm and the truncation algorithm. As an example, these algorithms are applied in the simulation of a traveling wave tube and a performance comparison is made among them. Finally, the conclusion is drawn in Section 4.

2. Analysis of Space Charge Field

2.1. Expression of the SCF

The Disc Model is a one-dimensional macro particle model, which divides an electron beam into a group of rigid discs along the axial direction, and each disc can be regarded as a “macro particle”, as shown in Figure 1.

The charge is uniformly distributed on the disc and remains unchanged in the process of motion. To describe the transcendence phenomenon, it is assumed that the discs can pass through each other. Then, the SCF is a statistical field generated by any two charged disks in the cylindrical channel. It can be solved by Green’s function.

$$\begin{array}{c} E_{sc}=\frac{Q}{2\pi {\epsilon }_0{r}_{\mathrm{b}}{}^2}Z\\ Z=\mathrm{sign}\left(z-z^{\prime }\right){\displaystyle \sum _{n=1}^{\infty }\left[\frac{2{\mathrm{J}}_1\left({\mu }_{0n}\frac{r_{\mathrm{b}}}{r_{\mathrm{c}}}\right)}{{\mu }_{0n}{\mathrm{J}}_1\left({\mu }_{0n}\right)}\right]}\exp \left(-{\mu }_{0n}\frac{\left|z-z^{\prime }\right|}{r_{\mathrm{c}}}\right)\end{array}$$

(1)

where $Q$ is the charge carried by each disk, $r_{\mathrm{c}}$ and $r_{\mathrm{b}}$ are the radius of the channel and the radius of the beam, ${\mathrm{J}}_1$ is the first order Bessel function of the first kind, ${\mu }_{0n}$ is the $n$th null point of the zero-order Bessel function of the first kind and $z-z^{\prime }$ is the displacement difference between any two discs. In Equation (1), there is a need to compute the sum of the series and solve the null points of the Bessel function, which can cause considerable computational costs. There is an approximate function that can replace $Z$ to simplify the Equation (2), i.e.,

$$Z=\mathrm{sign}\left(z-z^{\prime }\right)\exp (-\frac{2\left|z-z^{\prime }\right|}{r_{\mathrm{b}}})$$

(2)

To verify the rationality of this approximate expression, we can compare the original expression and the approximate expression, as shown in Figure 2. The approximate expression is close to the curves whose $r_{\mathrm{b}}/r_{\mathrm{c}}$ (radius ratio between beam and cavity) lies between 0.5~0.9. Generally, $r_{\mathrm{b}}/r_{\mathrm{c}}$ is about 0.5~0.7 in TWTs with a cylindrical beam channel. So, the approximate expression can take the place of Equation (1) and be applied to the solution of the SCF in TWTs.

The SCF of all source particles interacting with the objective particle in the vacuum tube is described by the following Equation.

$$E_{\mathrm{sc}}=\frac{Q}{2\pi {\epsilon }_0{r}_{\mathrm{b}}{}^2}{\displaystyle \underset{-\infty }{\overset{+\infty }{\int }}\mathrm{sign}\left(z-z^{\prime }\right)e^{-\frac{2\left|z-z^{\prime }\right|}{r_{\mathrm{b}}}}}d{z}^{\prime }$$

(3)

where $z$ is the position of the objective particle affected by the SCF and $z^{\prime }$ is the position of the source particle which produces the SCF. The working equations of the Lagrangian description for a TWT are

$$\begin{array}{l}{C}_1\frac{d^2{F}_n\left(\theta \right)}{d{\theta }^2}-j2n\frac{d{F}_n\left(\theta \right)}{d\theta }+n^2(b_n+C_1{b}_n-j{C}_1{d}_n)F_n\left(\theta \right)\\ \begin{array}{cccc}& & & \begin{array}{cc}& \end{array}=j2n^3{\left(\frac{C_n}{C_1}\right)}^3{(1+C_{1}b-jC{d}_n)}^3{I}_n\left(\theta \right)\end{array}\end{array}$$

(4)

$$\frac{{\partial }^2\mathsf{\Phi }\left(\theta ,{\phi }_0\right)}{\partial {\theta }^2}=-{\left(\frac{\partial \mathsf{\Phi }\left(\theta ,{\phi }_0\right)}{\partial \theta }+\frac{1}{C_1}\right)}^3\mathrm{Re}\left({\displaystyle \sum }_n{F}_n\left(\theta \right)e^{jn\mathsf{\Phi }}+F_{\mathrm{sc}}\right)$$

(5)

$$I_n\left(\theta \right)=\frac{1}{\pi }\underset{0}{\overset{2\pi }{{\displaystyle \int }}}{e}^{-jn\mathsf{\Phi }}d{\phi }_0$$

(6)

Equation (4) is the wave excitation equation, which describes the different harmonic electromagnetic waves excited by the corresponding harmonic beam current. $F_{\mathrm{sc}}$ is the normalized electric field. $n$ is the wave harmonic number. $\theta$ is the normalized distance along the tube. $C_n$ is the Pierce gain parameter and $C_1$ is the fundamental component. $b_n$ and $d_n$ are the nonsynchronous parameters and loss constant of different harmonics, respectively. $I_n$ is the harmonic component of the beam current. Equation (5) is the motion equation of the macro particles. Every macro-particle has its own equation. $\mathsf{\Phi }\left(\theta ,{\phi }_0\right)$ is the normalized-phase position of a specific macro-particle. ${\phi }_0$ is the entry phase of that macro-particle. $F_{\mathrm{sc}}$ is the space-charge field of that macro-particle, which is independent of other macro-particles. Equation (6) describes the beam current with different harmonic numbers. Equation (3) is changed to Equation (7) according to the same normalized method as Equations (4)–(6).

$$\begin{array}{l}{F}_{\mathrm{sc}}=\frac{{\omega }_{\mathrm{p}}^2}{2C_1^2{\omega }^2}\underset{-\infty }{\overset{+\infty }{{\displaystyle \int }}}\mathrm{sign}\left(\mathsf{\Phi }\left(\theta ,{\phi }_0^{\prime }\right)-\mathsf{\Phi }\left(\theta ,{\phi }_0^{}\right)\right)\cdot \\ \exp \left(-\frac{2}{r_{\mathrm{b}}}\cdot \left|\frac{1}{{\beta }_{\mathrm{e}}}\frac{1}{1+C_1\frac{d\mathsf{\Phi }\left(\theta ,{\phi }_0^{\prime }\right)}{d\theta }}\left(\mathsf{\Phi }\left(\theta ,{\phi }_0^{\prime }\right)-\mathsf{\Phi }\left(\theta ,{\phi }_0^{}\right)\right)\right|\right)d{\phi }_0^{\prime }\end{array}$$

(7)

$\mathsf{\Phi }\left(\theta ,{\phi }_0^{}\right)$ and $\mathsf{\Phi }\left(\theta ,{\phi }_0^{\prime }\right)$ are the phase of the objective particle and source particle. The discrete form of the space-charge field equation is

$$\begin{array}{l}{F}_{\mathrm{sc}}\left(\theta ,n\right)={\displaystyle {\displaystyle \sum }_{n^{\prime }=1}^{\infty }}\frac{{\omega }_{\mathrm{p}}^2}{2C_1^2{\omega }^2}\mathrm{sign}\left(\mathsf{\Phi }\left(\theta ,n^{\prime }\right)-\mathsf{\Phi }\left(\theta ,n\right)\right)\cdot \\ \exp \left(-\frac{2}{r_{\mathrm{b}}}\cdot \left|\frac{1}{{\beta }_{\mathrm{e}}}\frac{1}{1+C_1\frac{d\mathsf{\Phi }\left(\theta ,n^{\prime }\right)}{d\theta }}\left(\mathsf{\Phi }\left(\theta ,n^{\prime }\right)-\mathsf{\Phi }\left(\theta ,n\right)\right)\right|\right)\mathrm{\Delta }{u}_0\end{array}$$

(8)

where $n$ and $n^{\prime }$ are the serial number of the objective particle and the source particle, respectively, and $\mathrm{\Delta }{u}_0$ is the initial phase interval between two particles adjacent to each other in the injection port.

2.2. Analysis of the SCF

Equation (8) shows the SCF acting on the objective particle is produced by all other particles in the TWT. However, if the phase difference between two particles is greater than $\pi$, the interaction between them would be negligible. Therefore, only the source particle whose phase difference with the objective particle is within $(-\pi ,+\pi )$ is considered.

The incident particles are uniformly distributed within $(0,2\pi )$ at the initial moment. As the beam–wave interaction goes on, the incident particles may escape from $(0,2\pi )$ and the maximum phase difference between the incident particles would be larger than $2\pi$, as is shown in Figure 3a. Nevertheless, considering the periodicity of the particle distribution, we can make a periodic extension from other incident particles to solve this problem, as is shown in Figure 3b. The mathematical expression for this operation is shown as

$$\mathsf{\Phi }\left(n\right)=\{\begin{array}{c}\mathsf{\Phi }\left(n\right)-2\pi ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathsf{\Phi }\left(n\right)>2\pi \\ \mathsf{\Phi }\left(n\right)+2\pi ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathsf{\Phi }\left(n\right)<2\pi \end{array},(n=1,2\dots ,\mathrm{N})$$

(9)

To calculate the SCF of a specific objective particle, all source particles whose phase difference is within $\left(-\pi ,+\pi \right)$ are chosen from the extended particles. N is the total number of macro particles. This process is shown in Figure 4, and the mathematical expression for this process is given by

$$\mathrm{\Delta }\mathsf{\Phi }\left(n^{\prime },n\right)=\{\begin{array}{c}\mathrm{\Delta }\mathsf{\Phi }\left(n^{\prime },n\right)+2\pi ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }\mathsf{\Phi }<-\pi \\ \mathrm{\Delta }\mathsf{\Phi }\left(n^{\prime },n\right)-2\pi ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Delta }\mathsf{\Phi }>\pi \end{array},\begin{array}{c}\\ \end{array}(n=1,2,\dots N;n^{\prime }=1,2,\dots N)$$

(10)

These chosen particles are be used for calculating the SCF of a given objective particle. Therefore, to calculate the SCF of all particles, we need to repeat the above operations to obtain the chosen particles for different objective particles.

3. Solution of the SCF

3.1. The Traditional Algorithm to Calculate the SCF

Based on the above analysis of the particle phase distribution and Equation (8), a traditional algorithm is used to calculate the SCF first, whose schematic diagram is shown in Figure 5. This algorithm uses a major loop and a minor loop. The minor loop is used to calculate the SCF of a single objective particle, and the SCF excited by source particles ($\mathrm{\Delta }\mathsf{\Phi }\in (-\pi ,\pi )$) is calculated in sequence. The major loop is used for calculating the SCF of all objective particles. Then, the calculated value of the SCF is substituted into Equation (4) for solving the beam–wave interaction.

3.2. The Vectorized Algorithm

Because two nested-loop statements are involved, the complexity and the time cost greatly increase in the traditional algorithm. Here, a vectorized algorithm based on the traditional algorithm is proposed. The minor loop of the traditional algorithm is eliminated by vectorizing the input parameters $\mathsf{\Phi }$ and $d\mathsf{\Phi }/d\theta$. As shown in Figure 6, ${[\mathsf{\Phi }]}_N$ and ${[{\mathsf{\Phi }}^{\prime }]}_N$ are two N-dimensional vectors, in which both contain the phase information of N macro particles. Then, the time cost of the vectorized algorithm is only 0.4% of that of the traditional algorithm (as N = 12,000). In addition, the curves of Figure 7 show that the complexity of the traditional algorithm is about $O(n^2)$, while the complexity of vectorized algorithm is only about $O(n)$.

3.3. Parallelized Algorithm

In the major loop in the traditional algorithm, the SCF of each objective particle is calculated separately. Because the interaction between the objective particle and the source particles is independent, the calculation of the major loop could be parallelized. Therefore, based on the vectorized algorithm presented above, we propose two parallelized algorithms, which include CPU parallel computing and GPU parallel computing.

As an example, we used a CPU, Intel Core i7-9700K, for testing the parallelized algorithm. The Intel Core i7-9700K has 8 cores, and the basic frequency is 3.6 GHz. In the CPU-parallelized algorithm, we split the minor loop into several subloops based on the number of logical processors in the CPU, and the subloops were calculated simultaneously by the CPU, as shown in Figure 8.

Next, the GPU was used for the trial, since the GPU has thousands of scales of cores and has a better performance when performing large-scale calculations. In this case, we used a GPU, NVIDIA GeForce RTX 2060 SUPER for the testing of the GPU-parallelized algorithm, whose number of CUDA Cores is 2176. To realize the calculation of the SCF, the loop statements need to be eliminated by matrixing the algorithm. We constructed input parameters into matrices such as ${[\mathsf{\Phi }]}_{N\times N}$ and ${[{\mathsf{\Phi }}^{\prime }]}_{N\times N}$ (where ${\mathsf{\Phi }}^{\prime }=d\mathsf{\Phi }/d\theta$). Then, we obtained the phase difference matrix ${[\mathrm{\Delta }\mathsf{\Phi }]}_{N\times N}$ by ${[\mathrm{\Delta }\mathsf{\Phi }]}_{N\times N}={[\mathsf{\Phi }]}_{N\times N}-{[\mathrm{\Delta }\mathsf{\Phi }]}_{N\times N}^{\mathrm{T}}$, as shown in Figure 9. Finally, ${[\mathrm{\Delta }\mathsf{\Phi }]}_{N\times N}$ was substituted into Equation (8) for solving the SCF by matrix operations.

Compared with the vectorized algorithm, the CPU-parallelized algorithm significantly reduced the time cost by 74% (as N = 12,000), as shown in Figure 10. And the time cost of the GPU-parallelized algorithm was only 30% (as N = 12,000) of the CPU-parallelized algorithm. In addition, the time cost was reduced by nearly 90% (as N = 12,000) compared with that of the vectorized algorithm. The calculation time reduction in the GPU-parallelized algorithm was not notable compared with that of the CPU-parallelized algorithm, though the core number of the GPU was much greater than that of the CPU. This was because the data transmission between one CPU and other CPUs consumed as much time as the GPU is used. In addition, the computation capability of the single-core of the GPU was weaker than that of the single-core of the CPU. Furthermore, the GPU memory limited the particle number in the simulation process. As the particle number increased, the GPU memory usage also increased. When the particle number reached 13,000, the GPU memory usage reached saturation in our simulation. Therefore, using a more advanced GPU could increase the data volume for calculating.

3.4. The Truncation Algorithm

Because the SCF decreased exponentially as the distance between the source particle and the objective particle increased, the calculation efficiency could be enhanced by truncating the SCF induced by macro particles with a large distance. Here, the curve of the normalized SCF with respect to the distance between any two particles is presented in Figure 11. As the particle phase difference of less than $0.5\pi$ was considered, the normalized SCF was 0.008 (about 0.1%), where the interaction between two particles is negligible.

In addition, we analyzed the beam–wave interaction of a TWT under different truncated ranges of phase difference. We obtained the curve of the relative error of the output power, the most important parameter, with respect to ${\theta }_{\lim }$, which is the truncation radian of the phase difference (the range of the phase difference is $\mathrm{\Delta }\mathsf{\Phi }\in (-{\theta }_{\lim },+{\theta }_{\lim })$), as shown in Figure 12a. The relative error (RE) is given by $\mathrm{RE}=\left|P_{out}-P^{\prime }\right|/P_{out}$, where $P_{out}$ is the output power without truncation and $P^{\prime }$ is the output power with truncation. Figure 12 shows that the RE tended to be a constant value (about 12.5%) as ${\theta }_{\lim }$ approached 0, which is the case when the SCF is absent. As ${\theta }_{\lim }=0.5\pi$, the RE was only 0.73%. Moreover, the time cost was reduced by 23.8% (as the number of particles was 12,000), which is illustrated in Figure 12b. In this case, ${\theta }_{\lim }=0.5\pi$ was a feasible truncation radian of the phase difference for solving the beam–wave interaction.

4. Conclusions

In this paper, we presented three fast algorithms to calculate the SCF, i.e., the vectorized algorithm, the parallelized algorithm and the truncated algorithm, which were based on the traditional calculation algorithms. The time cost was reduced by 99.6% (when the number of particles was 12,000) for the vectorized algorithm. Then, we proposed parallelized algorithms which were based on a CPU and GPU. Compared to the vectorized algorithm, the time cost of the CPU-parallelized algorithm was reduced by ~74% and that of the GPU-parallelized algorithm was reduced by ~90%. Although the GPU-parallelized algorithm performed better than the CPU-parallelized algorithm, it had a memory usage higher than the latter. As the truncation radian of the phase difference was limited to $(-0.5\pi ,+0.5\pi )$, the relative error of the output power was only 0.73% and the time cost was reduced by 23.8% (when the number of particles was 12,000). Combining all these three algorithms and applying them to the SCF calculation, the time cost was reduced by three orders relative to that of the traditional algorithm (reduced to 0.048%). These algorithms, although used in a one-dimensional disk model, could be extended to two- and three-dimensional SCF calculations in order to increase the accuracy of vacuum electron device simulations. In addition, these algorithms could be applied to solve other SCF calculation problems in FEL, accelerators and plasma systems, etc.