WU Peiyu，XIE Yongjun，*，NIU Liqiang，and JIANG Haolin

1. School of Electronic and Information Engineering，Beihang University，Beijing 100191，China;2. School of Information Science and Engineering，Southeast University，Nanjing 210096，China

Abstract: Based on the finite element method (FEM) in the frequency domain and particle-in-cell approach in the time domain，a hybrid domain multipactor threshold prediction algorithm is proposed in this paper. The proposed algorithm has the advantages of the frequency domain and the time domain algorithms at the same time in terms of high computational accuracy and considerable computational efficiency. In addition，the compute unified device architecture (CUDA) acceleration technique also can be employed to further enhance its simulation efficiency.Numerical examples are carried out to demonstrate the effectiveness of the proposed algorithm. The results indicate that the multipactor threshold can be accurately predicted and the computational efficiency can be improved.

Keywords: compute unified device architecture (CUDA)，finite element method (FEM)，hybrid domain，multipactor threshold prediction，particle-in-cell (PIC).

1. Introduction

As one of the most well-known radio frequency breakdown phenomenon in the high frequency and vacuum environments，the multipactor effect plays an important part in the satellite payloads in space [1-3]. Multipactor is a nonlinear process caused by the increment of the secondary electrons emission (SEE) phenomenon in a certain condition [4]. With the development of the near-future space sector，the integration level of the satellite payloads becomes higher than ever before [5]. Thus，devices with a higher power capacity are severely demanded to realize such a condition. Especially for the high-power microwave components，such as circulators and waveguide microwave devices，the occurrence of multipactor will lead to a series of risks resulting in noise production，power dissipation，performance degradation and component destruction [1-3，6-8]. Therefore，the prevention of the multipactor effect is regarded as the frontier science in the space industries and its application.

In order to prevent the multipactor effect from the microwave components，the multipactor breakdown threshold voltage is employed initially [9]. However，with the development of microwave technology，the multipactor breakdown threshold voltage is inefficient in a series of microwave components，such as multi-ports devices. To alleviate such a problem，the multipactor breakdown threshold power is employed since then [10]. To predict the multipactor breakdown threshold power，the parallel plate approximation and the Monte-Carlo kinematics motion are first theoretically carried out [11]. Based on these methodologies，the multipactor susceptibility curve is introduced as a standard by the European Space Agency [12].However，the above-mentioned theoretical methods show their unacceptable accuracy in the calculation of complex structures. To alleviate that problem，the numerical methods are introduced during the simulation [13]. In 2005，the three-dimensional full-wave electromagnetic simulation tool (FEST 3D) was developed based on the velocity-verlet algorithm [13]. Through the experimental results and the practical applications，it can be verified that the FEST 3D is merely efficient by employing huge amount of simulation steps. Recently，the particle-in-cell(PIC) approach has been proposed based on the first principle of electron dynamics to minimize the computational time [14]. Combining the finite difference time domain (FDTD) with the PIC approaches，the PIC-FDTD approach shows the advantages of both the PIC and the FDTD approaches in terms of tracking the electrons’ trajectories and observing the waveform of the fields [15].As an explicit method，the computational efficiency of the PIC-FDTD approach is severely limited by the courant limit. The time step is severely limited by the mesh sizes. To further enhance the efficiency，the compu-te unified device architecture (CUDA) acceleration is carried out combined with the PIC-FDTD approach [16-17].

The magnetized ferrite with the unique non-reciprocity characteristic has been widely employed in the microwave components [18]. The field distribution obtained by the time-domain algorithm is inaccurate during the calculation especially in complex structures. To improve the computational accuracy，the frequency domain finite element method (FEM) is regarded as one of the most powerful methods to simulate the field distribution and electric properties of the magnetized ferrite devices [19].

Here，the hybrid domain algorithm based on the frequency domain FEM and time domain PIC-FDTD approach is proposed to predict the multipactor threshold.The computational efficiency can be further enhanced by adopting the CUDA acceleration. Numerical examples are carried out to validate the effectiveness and the efficiency. It is demonstrated through the results that combing the field distribution obtained by the FEM with the PIC-FDTD approach，the proposal shows its effectiveness.By employing the CUDA acceleration，the computation efficiency can be increased by dozens of times.

2. Simulation methodologies

During the simulation，the field distribution is calculated and simulated by employing the FEM in the frequency domain at the first step. At the second time step，the mesh information and the field distribution is exported according to meshes and phases of the simulation. In the frequency domain FEM，the tetrahedral division is employed in the calculation. To this point of view，the field distribution can be obtained in discrete values. To obtain the field component in the whole computational domain，the cubic spline interpolation method is employed between each mesh grid and the field distribution. Finally，the frequency domain field resultant is employed in the PIC-FDTD approach and its CUDA acceleration in the time domain. The flow diagram of the proposed algorithm is shown in Fig. 1.

Fig. 1 Flow diagram of the hybrid domain multipactor threshold prediction algorithm with CUDA acceleration

2.1 PIC approach

The PIC-FDTD approach is employed to simulate the movement of the particle by employing the mesh and field information which is calculated by the FEM. In the PIC-FDTD approach，the particles are simulated by the macroparticles.

Firstly，the macroparticles generate randomly in the component with the certain charge and mass. Then，the particles are driven by the Lorentz force，its equation can be given as

whereEis the electric flux density，Bis the magnetic flux density，mis the mass of the particles，qis the quantity of the particles，vis the speed of the particles which can be given as

wherexη(η =x，y，z) is the movement distance of the particles driven by the Lorentz force in η-direction. By employing the central difference in the conventional FDTD method，the update equation can be written as

wherenis the number of time steps，and the coefficients can be written asa1=1/Δtanda2=qBn/(2m). By employing (3) and (4)，the speed and the distance of the particles’ movement can be updated explicitly in the PICFDTD approach.

At the boundaries of the simulation，the judgement processing is introduced to demonstrate the SEE phenomenon. Three situations are included during the calculation which can be expressed as follows:

(i) The particles move out of the computational domain;

(ii) The particles bounce at the boundary with the same quantity and mass;

(iii) The SEE phenomenon occurs during the collision.

The flow diagram of the PIC-FDTD approach is shown in Fig. 2.

2.2 Modeling of secondary emission

In the high power microwave devices，the SEE phenomenon happens when the particles impact with the mediums inside devices. Meanwhile，the secondary electronics are produced through such interactions processing.

Fig. 2 Flow diagram of the PIC-FDTD approach

Thus，the modelling of the secondary emission in the simulation is of vital importance. The secondary emission yield (SEY) is deemed as one of the most important parameters during the procedures of the secondary electronics. The empirical models which include the Vaughan and Furman models are carried out to describe the SEY.Among the several models，the Furman model has been widely employed. The yield of the secondary electron expressed can be written as

where δ is the yield of the secondary electron，δmaxis the SEE maximum produce value，θis the angle of incidence electron，ksis the surface smoothness factor andwcan be given as

whereEtis the threshold power of incident electronics andEiis the power of incident electronics. The relationship betweenwandkcan be given as

The power of outgoing electronics can be given as

whereEmaxis the power under the δmaxsituation.

Fig. 3 and Fig. 4 show the SEY of silver and aluminum obtained by employing the Furman model.

Fig. 3 SEY curves of silver obtained by the Furman model

Fig. 4 SEY curves of aluminum obtained by the Furman model

2.3 CUDA acceleration

Owing to its unique architectural features，the CUDA programming model is especially suitable for tackling problems that can be viewed as data-parallel tasks. On the other hand，owing to its nature as a massively parallel algorithm，the FDTD simulation can be treated as the dataparallel task via mapping per graphics processing units(GPU) thread to a Yee cell’s computation [20-21].

To implement the parallel FDTD method based on the CUDA technique efficiently，the entire code running on GPU can be divided into two parts. To be specific，one kernel is used to update the electritic field，while the oth-er one is used for calculating the magnetic field. Moreover，instead of using the global memory singly，shared memory and constant memory are both implemented in our proposed algorithm to further enhance its simulation efficiency.

3. Numerical results

To demonstrate the effectiveness of the proposed approaches，the parallel metal plate is simulated and calculated firstly. Then，the ridge waveguide model is carried out to further demonstrate the effectiveness and efficiency. Finally，a more complex anisotropic magnetized ferrite circulator structure is introduced. In this paper，a PC with Intel CoreTMi7-6700k 3.20 GHz，64 GB (DDR4，2 666 MHz)and a GPU with Nvidia GeForce RTX 2 080 Ti (11 GHz)is employed.

3.1 Effectiveness of the PIC approach and its CUDA acceleration

The effectiveness of the PIC-FDTD approach and its CUDA acceleration，denoted as PIC-CUDA-FDTD，is testified through a metal parallel plate which has the analytic solution in theory.

As shown in Fig. 5，supposing the length of the parallel copper plate is 5 cm inx-direction. Particles generate randomly on the metal plates with the parameters ofks=1，δmax=1.3 andErf0=600 eV. The distance between the lower boundary and the upper boundary can be described as the functiond. In this simulation，d=1 cm is chosen. The field between the parallel metal plates can be written as

whereErf0is the maximum value of the field.vr，vrxandvryare the speed of the particle，the speed of the particle inx-direction and the speed of the particle iny-direction，respectively.

Fig. 5 Model of the PIC-FDTD and PIC-CUDA-FDTD approaches simulation in the parallel copper plate

Fig. 6 shows the number of secondary electronics whenErf0=520V obtained by employing the PIC-FDTD and the PIC-CUDA-FDTD approaches. It can be observed that the curves obtained by different approaches and the theoretical solution are almost overlapped. Meanwhile，it can be observed that the accuracy of the PIC-CUDAFDTD approach is inferior compared with that of the PICFDTD approach.

Fig. 6 Number of the secondary electronics obtained by the PICFDTD，PIC-CUDA-FDTD approaches and the theoretical solution

To further investigate the accuracy of the calculation，the relative local error is introduced for comparison which can be written as

whereNtis the secondary electronics obtained by employing the proposed algorithms andNris the reference solution which is the theoretical solution of the equations.Fig. 7 shows the relative local error between the different algorithms and the theoretical solution. It can be observed that the computational accuracy of the PIC-CUDAFDTD decreases significantly during the simulation.However，the PIC-CUDA-FDTD approach still maintains a considerable performance with the highest relative local error of -47 dB (0.01%). As shown in Table 1，the computational efficiency can be increased by employing PIC-CUDA-FDTD significantly which can be increased by 36% compared with that of the PIC-FDTD.

As shown from the numerical example in the parallel metal plate，the effectiveness and the efficiency of the PIC-CUDA-FDTD approach have been demonstrated.The results show that the proposed approaches can not only obtain a considerable accuracy but also receive a higher computational efficiency.

Fig. 7 Relative local error obtained by the PIC-FDTD and PICCUDA-FDTD approaches

Table 1 Running time of different algorithms

3.2 Multipactor prediction for the ridge waveguide model

To further demonstrate the effectiveness of the proposed algorithm，a ridge waveguide model is carried out.As shown in Fig. 8，the waveguide is with the size of 60mm×20mm×10mminx-，y- andz-directions，respectively. The stages are with the size of12mm×20mminx- andy-directions and 2 0mm×1mm iny- andz-directions，respectively. The plane wave is a Gaussian pulse with the maximum frequency of 10 GHz. The walls of the waveguide structure can be deemed as the perfect E condition during the calculation. Electronics are formed randomly on the walls to investigate the multipactor threshold.

Fig. 8 Model of the ridge waveguide structure

The field of the waveguide structure is simulated and calculated by employing the FEM in the frequency domain. The multipactor threshold is simulated by employing the PIC-FDTD and PIC-CUDA-FDTD approaches，respectively. Table 2 shows the multipactor threshold power，the running time and the time reduction. The multipactor threshold is predicted by performing the proposal approaches for several times. As the number of particles reaches a constant in the late-time during the simulation，the input power can be deemed as the multipactor threshold.

Table 2 Time and multipactor threshold prediction of different algorithms

As shown in Table 2，compared with the commercial software and the PIC-FDTD approach，the PIC-CUDAFDTD approach holds the best computational efficiency.The number of particles versus time with different powers obtained by employing the PIC-CUDA-FDTD approach is shown in Fig. 9. It can be observed that the number of particles becomes a constant in the late-time during the calculation.

Fig. 9 Number of particles obtained by employing the PIC-CUDAFDTD approach at different input powers

3.3 Multipactor prediction for the magnetized ferrite circulator structure

To illustrate the proposed algorithm in the microwave components with complex structures and several mediums，the ferrite circulator structure is employed to investigate the multipactor threshold.

Due to the complex field distribution of the anisotropic magnetized ferrite medium，large calculation error will occur by employing the conventional time domain algorithms directly. To alleviate this problem，the field distribution of the circulator is simulated by employing the FEM in the frequency domain. By this means，the accuracy of the computation will be improved compared with that using the time domain algorithms.

The model of the circulator structure is shown in Fig. 10.Three waveguide ports are employed with the standard of BJ40. The center of the circulator is composed of five layers in the vertical directions. The top and the bottom layers are equilateral triangles with the side length of 40 mm which are made up of aluminum. The thickness of the layers is 0.5 mm. The cylinder in the center is made up of three layers including silver，ferrite and silver which are with the same thickness of 10 mm. It can be easily observed that the multipactor sensitive region is between the top and the bottom aluminum layers. The front view of the multipactor sensitive region in the circulator structure is shown in Fig. 11.

Fig. 10 Model of the ferrite circulator with three waveguide ports

Fig. 11 Model of front view of the multipactor sensitive region of ferrite circulator

A plane wave which is a Gaussian pulse with the maximum frequency of 10 GHz incidents from port 1 of the circulator structure. The mesh size can be obtained as 0.15 mm. The time step is 3 fs. In order to simulate different input powers during the calculation，different amplifies are employed to investigate the multipactor threshold. The prediction of the multipactor threshold are valued through taking the average value of the input power by several times. The electritic fields are calculated by the FEM in the frequency domain. Fig. 12 shows the electritic field distribution of the ferrite circulator structure from the top view at the multipactor threshold power. Fig. 13 shows the number of particles and its fitting curves versus different input powers obtained by employing the PICCUDA-FDTD approach.

Fig. 12 Field distribution of the ferrite circulator structure from the top view at the multipactor threshold power

Fig. 13 Number of particles and its fitting curves versus different input powers obtained by PIC-CUDA-FDTD approach

As shown in Fig. 13，it can be easily predicted that the multipactor threshold is 11 000 W. In order to further demonstrate the computational accuracy of the PICCUDA- FDTD approach，the multipactor threshold prediction obtained by PIC-FDTD and PIC-CUDA-FDTD is shown in Fig. 14. The vertical axis is the electronic numbers which indicates the multipactor threshold. It can be observed that the above-mentioned approaches can hold the considerable accuracy. Table 3 shows the multipactor threshold，the running time and the time reduction obtained by different algorithms.

Fig. 14 Fitting curves of particles obtained by PIC-FDTD and PICCUDA-FDTD approaches at multipactor threshold power

Table 3 Time and multipactor threshold prediction of different algorithms in the multipactor structure

As can be concluded from the results，the PIC-CUDAFDTD approach shows its good performance and considerable computational efficiency.

4. Conclusions

Based upon the frequency domain FEM and the time domain PIC-FDTD approach，the hybrid domain multipactor prediction algorithm is proposed to simulate the field distribution，calculate the number of particles and predict the multipactor threshold. The computational efficiency of the proposed algorithm can be enhanced by employing the CUDA acceleration. As is demonstrated in the parallel plane plate，the waveguide model and the circulator structure，the proposed algorithm shows the advantages of the time domain algorithm，the frequency domain algorithm and the CUDA acceleration in terms of improving the computational accuracy in complex structures，enhancing the computational efficiency and directly obtaining the field and particles distribution.