Multipactor radiation analysis within a waveguide region based on a frequency-domain representation of the dynamics of charged particles

A technique for the accurate computation of the electromagnetic fields radiated by a charged particle moving within a parallel-plate waveguide is presented. Based on a transformation of the time-varying current density of the particle into a time-harmonic current density, this technique allows the evaluation of the radiated electromagnetic fields both in the frequency and time domains, as well as in the nearand far-field regions. For this purpose, several accelerated versions of the parallel-plate Green’s function in the frequency domain have been considered. The theory has been successfully applied to the multipactor discharge occurring within a two metal-plates region. The proposed formulation has been tested with a particle-in-cell code based on the finitedifference time-domain method, obtaining good agreement.


I. INTRODUCTION
The radiation of a charged particle moving within a certain metallic and dielectric waveguide structure has deserved the attention of many researchers in different fields of electromagnetic theory. Several authors have investigated the radiation of charges and beams in RF accelerating structures. Concretely, the study of the radiated fields left behind by a charged particle, which are called the wake fields, are particularly important because they influence the motion of the charged particles that follow them. The rigorous analysis of the wake fields is a nontrivial issue, and critically depends on the geometry and materials of the accelerator structure. In most of the current models found in the technical literature, the particles move linearly at constant velocity ͓1-8͔.
Another physical phenomenon directly related with the radiation of electrons moving within an RF waveguide passive component is the multipactor effect. Multipactor discharge is a nonlinear effect that may occur in high power microwave devices at very low pressures, especially in satellite subsystems ͓9-22͔ and particle accelerators ͓23-33͔.
Basically, multipactor is an electron avalanche in synchronism with the existing RF field caused by secondary-emission multiplication on the device walls: primary electrons accelerated by the RF fields impact on a surface releasing a larger number of secondary electrons, which may in turn be accelerated by the RF fields striking successively against the other surface, then extracting more and more electrons, and so on. The discharge, which is normally undesired, can occur for a wide range of frequencies ͑from megahertz to tens of gigahertz͒ and in a wide set of geometries. In microwave devices for space communications, its effects range from the signal degradation ͑increasing noise levels͒ to the complete destruction of the component. In accelerating structures, multipactor causes serious problems during accelerator operation ͑increasing the reflection coeffi-cient͒, and can generate dielectric window breakdown on the vacuum side.
The radiated power spectrum of a multipactor discharge between two parallel plates has already been investigated by the authors ͓34͔. In such work, the multipactor discharge was first approached as a spatially uniform electric current. However, in this paper we propose to transform the time-varying current density of the particle into a time-harmonic current density. Such transformation allows the accurate computation of the harmonic fields radiated by the charged particle, both in the near-and the far-field regions. Thus, the main goal of this contribution is the evaluation of the electric and magnetic fields emitted by a multipactor discharge, as well as * benito.gimeno@uv.es the radiated power spectrum, using an efficient technique which transforms the time-varying electron motion into frequency-domain equivalent currents. For such purpose, several accelerated versions for the frequency-domain Green's function of the considered metallic enclosure are used. In particular, the theory has been applied to the multipactor analysis in a parallel-plate waveguide following the well-known Sombrin's model ͓35͔, where a single effective electron is moving in one dimension under the influence of an external RF harmonic electric field.
The paper has been organized in three sections. In the next section, the frequency-domain representation of a charged particle motion is developed in detail. Afterward, we apply this theory to the multipactor analysis in a classical parallel-plate scenario, where the Green's function is expressed as an infinite series using the classical image theory, whose computation is efficiently accelerated through the Ewald's and Kummer's techniques. Next, numerical results obtained with our theory are successfully compared with results from a particle-in-cell code, thus validating the presented formulation. Finally, conclusions and future work are summarized.

II. FREQUENCY-DOMAIN REPRESENTATION OF A CHARGED PARTICLE MOTION
Let us consider a charged particle arbitrarily moving in one dimension ͑z axis͒ under the influence of an external force. From an electromagnetic point of view, this charge can be modeled as an electric density of charge and a current density J z , which depend on position ͑z͒ and time ͑t͒, ͑z,t͒ = q␦"z − z q ͑t͒…, where ␦ is the Dirac delta distribution, q is the particle charge, and z q ͑t͒ and v q ͑t͒ are its position and velocity, respectively, which obviously satisfy v q = ż q ͑t͒; note that dot indicates time derivative. It is well known that both singular densities are related to each other by the continuity equation. In this formulation we suppose that both trajectory z q = z q ͑t͒, and velocity v q = v q ͑t͒, of the particle are known functions.
As it is mentioned in the Introduction, our aim is to transform the singular electric current J z into a time-harmonic function depending on the position coordinate z. To proceed, we need to express the Dirac delta distribution in terms of a time-dependent function. First of all, we express the argument of the Dirac delta as an auxiliary function defined as follows: z ͑t͒ϵz − z q ͑t͒. Next, we expand the Dirac delta of such auxiliary function as a finite summation of timedependent Dirac elta distributions using the change of variable operation described in ͓36͔ ͑see Sec. 1.7͒, where t i , i =1,2, ... M, corresponds to the zeros of the function z ͑t͒, i.e., z ͑t i ͒ = 0, which occur when the particle is located at the position z: z q ͑t i ͒ = z. The following condition must also be satisfied for a consistent definition of Eq. ͑1͒: ͉ z ͑t i ͉͒ 0 for i =1,2, ... M. As a consequence, ͉ z ͑t i ͉͒ = ͉v q ͑t i ͉͒, which allows to rewrite the current density in this form.

͑2͒
For multipactor applications under perfect resonant conditions, we can consider that the particle velocity motion is a periodic function of period T q ͑which is directly related to the multipactor discharge order and the period of the harmonic external force͒. Therefore, the electric current density can be expanded in terms of the classical Fourier series as follows:

͑3͒
where a m and b m are unknown coefficients, and the angular frequency defined as q =2 / T q has been used. It is evident that if the motion of the particle is not a periodic function, the Fourier transform integral should be used instead of the Fourier series. Taking into account Eq. ͑2͒, these coefficients are trivially calculated, resulting in where the sign function defined as sgn͑x͒ϵx / ͉x͉ has been introduced. We emphasize that both coefficients a m and b m depend on the z coordinate because the zeros of z ͑t͒ directly depend on the considered z point. In order to use complex phasor forms, it is more convenient to define another In practice, it is more convenient to represent the current density using complex phasor forms defined as i m ͑z͒ where Re denotes real part. As a conclusion, in this section we have developed a time-harmonic representation of the charge motion in terms of harmonic currents i m ͑z͒, which only depend on the geometrical variable z.

A. Multipactor model
In order to implement the previous theory, we will apply it to the analysis of a multipactor discharge occurring within a parallel-plate region of gap d, as shown in Fig. 1. Several approaches have been developed in the technical literature to describe the multipactor effect, as we detailed in the Introduction. As a first simple model of this double-surface multipactor problem we will use Sombrin's model ͓35͔, which is based on the consideration of a unique effective electron which accumulates the total charge and mass of the generated electrons. Based on the direct integration of the Lorentz force, this model allows an analytical description of the motion equations and the resonance conditions, which will be key in the development of the present derivation. Other relevant parameters, such as the RF breakdown threshold voltage and the susceptibility charts, can be easily obtained using Sombrin's model. In order to obtain an analytical description of the physical phenomena, Coulombian repulsion of the electron cloud is not considered, and the formulation has been derived for one single electron, which imposes some limits in the application of the model. According to multipactor theory, during the discharge the charge grows exponentially as where n e ͑t͒ is the number of elementary charges, n e ͑0͒ is the number of charges at t = 0, and is the effective secondaryemission yield ͑SEY͒. Nevertheless, for multipactor regions close to the breakdown level ͑ =1͒, or high multipactor orders ͑high N͒, the charge grows slowly enough to be considered constant, n e ͑t͒Ϸn e ͑0͒, during a limited period of time. Therefore the theory developed hereafter for one single electron ͓n e ͑0͒ = 1 without loss of generality͔ remains valid within the proximity of the breakdown regions, which are indeed of the greatest interest for most applications. To proceed, we will suppose an electron with charge q =−e under the action of a z-oriented RF harmonic electric field given by E z ͑t͒ = E 0 cos͑t͒, f = / ͑2͒ being the RF frequency, T =1/ f its period, and E 0 the field amplitude, related with the voltage amplitude V 0 by E 0 = V 0 / d; obviously the driven voltage is V͑t͒ = V 0 cos͑t͒. In our model, an effective electron is launched from the lower plate ͑located at z = 0 plane͒ at time t = t ␣ with initial velocity v q 0 = v q ͑t ␣ ͒. The Lorentz force law for this simple case is given by the following equation: F z = qE z ͑t͒ = m q a z , m q being the electron mass and a z its acceleration. Thus, its first integration gives the velocity v q = v q ͑t͒, and the second one provides the electron position z q = z q ͑t͒. Under multipactor perfect resonance conditions, the transit time to the second surface ͑located at z = d plate͒ is N half RF cycles ͑N odd͒ ͓11͔. Under these considerations, the electron velocity is a periodic function with period T q = NT, given by the following expression: In Sombrin's model it is assumed that the magnitude of the initial velocity ͉v q 0 ͉ is constant in all launches ͓35͔.
As an example, we have simulated two multipactor discharges for orders N = 1 and N = 7, respectively; their more relevant parameters are summarized in Table I  both multipactor orders N = 1 and N = 7, respectively. As we observe in these graphics, the current patterns present more oscillations for higher harmonics, as predicted by Fourier analysis. It is worth noticing that these patterns are similar to the current distribution in a dipole antenna, of electrical length l e = m / 2, where is the free-space wavelength of the RF signal. For instance, the current for the m = 1 coefficient is similar to the current induced in a half wavelength dipole. The current for the m = 2 coefficient is similar to the current induced in a dipole one wavelength long.
We want to remark that the zeros of the z ͑t͒ function have been numerically calculated by means of the bisection technique ͓37͔.

C. Radiated fields by the harmonic currents
Our final objective is to evaluate to radiated fields of each harmonic current i m in the parallel-plate waveguide region. In order to proceed, we will use the classical radiation theory for time-harmonic fields, based on the use of the electric and magnetic auxiliary vector potentials, as detailed in the technical literature ͓38͔. For an electric current oriented in the z axis, only a z component of the magnetic vector potential exist, which satisfies the inhomogeneous Helmholtz equation, where A z m is the magnetic vector potential of the mth harmonic, k m is its wave number, and J z m is the z-oriented cur-rent density. The wave number of each harmonic is expressed as k m = m͑ q / c͒ = m͑k / N͒, k being the free-space RF wave number given by k = / c. Note that the term exp͑im q t͒ has been omitted throughout the formulation, and conventional cylindrical coordinates ͑r , , z͒ are used. The solution of Eq. ͑9͒ will be performed using the Green's function technique. The unknown vector potential of the problem can be expressed in terms of the following integral ͓39͔, The next step in the formulation deals with the evaluation of the Green's function of the problem g͑r ជ , r ជ Ј͒, which satisfies the scalar Helmholtz equation, with the Neumann boundary condition. The physical configuration of our problem consists of a time-harmonic z-directed electric point current element of value 1 A / m 2 , located at the source point r ជ Ј= zЈẑ, as depicted in Fig. 6. The vector potential generated by such current element ͓g͑r ជ , r ជ Ј͔͒ obviously satisfies Eq. ͑11͒. The method of images ͓40͔ provides the solution of this problem as an infinite spatial image series,

͑12͒
as shown in Fig. 6. Each of these elemental point currents radiate in free space; thus, the total field ͑in the region of interest͒ is the superposition of the radiated field by each infinitesimal element, where R sp is the magnitude of the relative vector position R ជ sp = r ជ − r ជ sp Ј , i.e., R sp = ʈR ជ sp ʈ. It should be emphasized that the singularity in this series is contained in the term s = 0 and p = + 1. Therefore, the magnetic vector potential can be calculated by means of the following integral: where we have taken J z m = i m , which is the case of our wire currents. Finally, the time-harmonic magnetic field can be evaluated through the well-known expression whereas the electric field is calculated using E ជ m = E m r r + E m z ẑ =−ٌ m − im q A ជ m , m being the electric scalar potential, resulting in The Lorentz gauge ͓39͔ has been used above to obtain the electric scalar potential: In order to evaluate numerically the integrals involved in Eqs. ͑15͒-͑17͒, the Gauss-Legendre quadrature algorithm has been used ͓37͔. The computation of the Green's function and its spatial derivatives is a difficult task, as the infinite series of spatial images shown in Eq. ͑13͒ converges very slowly. In most cases, a large number of terms is needed for achieving acceptable accuracy in the value of g͑r ជ , zЈẑ͒. Therefore, the evaluation of A z m and its derivatives for several hundreds of points becomes impractical. To solve this problem, some acceleration techniques have been applied to these series. First, by means of the Poisson formula, it is possible to transform Eq. ͑13͒ into an equivalent series written in the spectral domain, In Eq. ͑18͒, H 0 ͑2͒ is the Hankel function of zeroth order and second kind, which is defined as a combination of the Bessel functions. The Hankel function rapidly vanishes for large values of the spatial distance of the transversal coordinates ͑r͒, allowing therefore a very fast convergence for this new form of the series. However, if ͑r͒ tends to zero, the singularity of the Bessel functions appears, and the convergence rates are strongly affected in a negative way. For this situation, the authors have applied the Ewald method ͓41͔, which is a very efficient technique for accelerating this kind of series, even under singularity conditions ͓42-44͔. For our particular Green's function, the Ewald technique needs to be reformulated, thus obtaining the following expressions: g͑r ជ,zЈẑ͒ = g 1 + g 2 , ͑19a͒ 2d ͚ n=0 +ϱ ⑀ n I n ͑E͒cos͑k zn z͒cos͑k zn zЈ͒, ͑19c͒ In the previous expressions, erfc is the complementary error function. It can be noticed that this technique decomposes the original series into two new series, each one related to the spatial and the spectral domains. The value E is the so-called splitting parameter of the Ewald method, which controls the convergence of the two series. An optimal value of this parameter must be obtained for the correct behavior of the whole method, and several works have been developed concerning this matter ͓43,45͔. For our purposes, this value will be set as E = ͱ / d. The presence of the erfc functions ensures Gaussian convergence for the g 1 series, once this splitting parameter is properly selected. On the other hand, the term g 2 also converges very fast, provided that the integral I n is efficiently evaluated. The propagation constant ␥ n plays a very important role in the calculation of this integral. If k zn Ͼ k, which is the case of nonpropagative modes of the parallel plates, ␥ n 2 has a negative value, and the whole integral can be evaluated numerically without effort. On the contrary, for the complementary case of propagative modes, it is very important to ensure that Re͓ ␥ n 2 s 2 ͔ Ͻ 0 for the integral to be convergent ͓43͔. It is possible to fulfill this requirement if an adequate integration path is chosen, which should satisfy the following condition in the complex plane ͓43͔: where ␥ n and s denote the phases of ␥ n and the complex integration variable s, respectively. Figure 7 shows with dashed lines a generic path that fulfills condition ͑20͒. The authors have chosen a simple strategy based on a triangular trajectory, also shown in the same figure. This path allows splitting the integral I n into two linear integrals, which are very easy to evaluate numerically. Once I n is computed for all cases, the g 2 series rapidly converges, and the total Green's functions are obtained employing very few terms, even for the critical case of r → 0. Finally, in order to further reduce the computational effort, the authors propose the following algorithm based on the combination of several techniques for the efficient calculation of the Green's functions. First, the Ewald form of the series is only used in the near-singularity cases, and it is switched to the spectral form of the series for any other case. This is convenient since the evaluation of Hankel functions requires less CPU effort than the computation of Erfc functions. By choosing to evaluate Eqs. ͑18͒ or ͑19͒, depending on the value of the spatial distance, the calculation of the Green's functions is performed very efficiently. The authors have found that good performance is obtained by switching to the Ewald approach if r Ͻ 0.1. For the computation of the total electric and magnetic fields, the corresponding spatial derivatives are also needed. They can be obtained by directly operating with expressions ͑18͒ and ͑19͒, obtaining the corresponding accelerated form of the series for each case. Once all the different Green's functions have been evaluated properly, the electromagnetic fields can be computed very efficiently following this approach.

D. Numerical analysis of a multipactor discharge
We proceed to examine the field distributions of a multipactor discharge occurring in a parallel-plate scenario for the cases described in Table I the z component of the electric field E z m is plotted as a function of the normalized radial component r / , for z = d / 2 and multipactor orders N = 1 and N = 7, respectively. The results corresponding to the first Fourier harmonics are displayed. It can be seen that the electric field tends rapidly to zero for a distance of several wavelengths. Similar results have been observed for the imaginary part, which are omitted for brevity. We want to emphasize that the electromagnetic fields related with the even Fourier coefficients ͑m =2,4,6,...͒ in Eqs. ͑16͒ and ͑17͒ vanish, even though such coefficients are nonzero ͑see Figs. 4 and 5͒, which is in agreement with previously published simulations and experiments ͓46,47͔.
Second, the time evolution of the previous multipactor discharges has been compared with results obtained with an electromagnetic particle-in-cell ͑PIC͒ code based on the finite-difference time-domain technique ͓47͔. In order to compare with the analytical model, the numerical simulations were carried out with a constant charge scheme. Figure  10 depicts the z component of the electric field as a function of time, showing a good agreement between our simulations and the results provided by the PIC code.
Finally, the power generated by the discharge has been evaluated as a function of the frequency. In order to proceed, the real part of the complex Poynting vector ͑defined as N ជ = ͑1 / 2͒E ជ ϫ H ជ ‫ء‬ , ‫ء‬ denotes complex-conjugate operation͒ has been integrated over a cylindrical surface of height d and radius r = ; such cylinder is centered on the z axis. Figures 11 and 12 show the total power spectrum radiated by the multipactor discharges for the two cases of Table I, with multipactor order N = 1 and N = 7, respectively. Comparison with numerical results, computed with the cited PIC code ͓47͔ fully validates the presented theory. As stated before, only the odd harmonics, located at frequencies mf / N, carry some power. It can be seen that, for both test cases, the maximum emitted power corresponds to the harmonic located at the driven RF frequency, which is m = N, and agrees with previously published results ͓34͔. For the rest of harmonics the power level decreases as 1 / m.

IV. CONCLUSIONS
In this paper, we have proposed a technique to compute the electromagnetic radiated fields of a charged particle moving within a metallic enclosure. Such method is based on the transformation of the time-domain charged particle movement into an equivalent point-line time-harmonic current, which radiates into the considered environment. In particular, this formulation has been applied to study a multipactor discharge ͑simulated with Sombrin's approach͒ occurring in a parallel-plate waveguide region. A standard Fourier analysis has allowed the transformation of the periodic electron movement into equivalent spatial-dependent harmonic currents. The Green's function in frequency domain has been used to analyze the radiation of such currents within the parallel-plate guide, thus obtaining the total radiated fields both in the frequency and time domains. Comparisons with a PIC code fully validate the proposed theory.
Even though this formulation has been derived for a onedimension movement, it can be extended to the general case of a charged particle moving into a more realistic threedimensional scenario. In this way, this theory might be used to calculate the radiated fields in other more complex contexts, such as wake fields in RF accelerating structures, multipacting in microwave waveguides and other similar problems.