A Genetically Based Algorithm to Improve Execution Speed in Multipactor Simulations in Parallel-Plate Waveguides

This article presents an approach to enhance the convergence speed of Monte Carlo simulations for multipactor phenomena in parallel-plate geometries. Multipactor, a self-sustained electron discharge, may cause significantly harmful effects in electronic device operation. An adaptive genetic modification method is introduced to explore the parameter space of electron trajectories efficiently, leading to rapid identification of the critical combinations that induce population growth or decay. The proposed approach involves modifying key parameters, such as the initial phases and kinetic energies of individual particles, based on their performance in multipactor simulations. The reported results demonstrate a substantial improvement in convergence speed, achieving accurate results with fewer tracked particles. The adaptive genetic modification proves to be a promising technique for fast multipactor threshold predictions in parallel-plate configurations and has potential applications when analyzing electron-induced phenomena in a wide range of electronic systems.


A Genetically Based Algorithm to Improve Execution Speed in Multipactor Simulations
in Parallel-Plate Waveguides

I. INTRODUCTION
M ULTIPACTOR is an undesired electron avalanche-like discharge [1], [2], [3], [4], [5], which occurs when free electrons inside a device operating under high vacuum conditions move synchronously, driven by the RF fields within the microwave component under analysis.When the electrons impact against the metallic walls of the component, with kinetic energy able to release enough secondary electrons from the surface of the walls, the electron population can increase exponentially and, eventually, an electronic discharge occurs.Multipactor discharges can lead to several effects that are extremely harmful to the component performance [6].Some of these effects are: power dissipation through heating of the device walls, out-gassing, increase of reflected signal and noise floor, and detuning of resonant cavities.The consequences of the combination of these effects can even lead to physical damage in the device.
Multipactor can occur in a wide variety of systems.This article focuses on RF and microwave components designed to operate in space [6], [7], [8], but it can also appear inside klystrons [9], magnetrons [4], [10], particle accelerators [11] and many other devices, for which the techniques presented in this work can be also extended to.
Nowadays, the trend in satellite communications is to improve the data rate throughput by means of an increase of the transmitted power, as the coding and modulation speeds operate near the theoretical maximums established by the Shannon theorem [12].In this context, estimating the multipactor threshold is of paramount importance in order to ensure a proper operation for space applications.Usually, the power threshold can be assessed through experimental verification [13].Nevertheless, considering that experimental tests are time-consuming and expensive, simulation tools have become increasingly important, especially in the design phase of RF components [14], [15], [16].
Achieving statistically convergent results, however, can be challenging due to the presence of multiple random variables such as the emission energy and angle of secondary electrons.As a result, a large particle population must be tracked, leading to time-consuming simulations.To alleviate this situation, an adaptive algorithm that makes use of genetically modified particles is proposed in this work.By dynamically modifying the properties of each particle (genetic modification), we aim to demonstrate the potential for speed-up when compared to traditional simulation methods.The proposed approach has been applied to parallel-plate waveguides, as a proof of concept, since particle trajectories can be expressed analytically because of the simplicity of the corresponding electromagnetic (EM) fields.This research, therefore, aims to explore the use of genetic algorithms to enhance the computational efficiency of multipactor simulations.
In Section II, the theory behind the particle dynamics within a parallel-plate waveguide under RF excitation is revisited.Next, in Section III, the fundamentals of the proposed genetic algorithm for multipactor analysis are presented.Section IV focuses on the comparison of the computational performance of the traditional Monte Carlo method and the proposed approach.Finally, the study concludes in Section V, highlighting the potential effectiveness of genetic algorithms to enhance the convergence speed of multipactor analysis techniques.

II. MONTE CARLO FOR MULTIPACTOR SIMULATION
Consider a 2-D parallel-plate waveguide geometry, as depicted in Fig. 1.The two metal surfaces are separated by a distance d and are subjected to excitation from an RF voltage of the form v(t) = V 0 cos(ωt + α), where V 0 represents the voltage amplitude, ω = 2π f is the angular frequency, and α denotes the initial phase at time t = 0 s.This waveguide supports a fundamental TEM mode solution.For this mode, the electric field is perpendicular to both metal conductors, and is given by a time-varying expression of the form ⃗ E = −(V 0 /d) cos(ωt + α) ŷ = E 0 cos(ωt + α) ŷ.A low-energy electron, located inside the vacuum gap, is accelerated due to the uniform RF field toward one of the metal surfaces according to the electric field polarity.Upon collision, depending on the impact energy, emission of secondary electrons into the vacuum medium, or absorption of the primary electron by the metal wall, may occur [5], [34].The mean number of released electrons after each collision is determined by the secondary electron yield (SEY) δ of the surface.A value of δ > 1 indicates a net emission of electrons, and therefore an increase in the electron population, whereas δ < 1 implies a net electron absorption.Crucial parameters defining the SEY curve are the first crossover energy W 1 and W max , which are associated with the impact energies where δ = 1 and δ = δ max , respectively.For this work, the SEY function employs the modified Vaughan model [35] using the parameters described in the ECSS standard [36] for aluminum (W 1 = 17 eV, W max = 276 eV, δ max = 2.92, and δ low = 0.8).It is assumed that for those collisions with impact energies below W 0 = 10 eV, the electrons are elastically reflected with a SEY value of δ = δ low .
The algorithm tracks the trajectory of particles or effective electrons, being N i the number of electrons composing the ith particle.The dynamics of each particle across the plates can be derived from the nonrelativistic Lorentz force expression where y i represents the y-axis position of the ith effective electron, and e and m denote the charge and mass of an electron, respectively.Note that the magnetic field has been omitted in (1) for its negligible effect with respect to the electric field |⃗ v × ⃗ B| ≪ | ⃗ E| , as typically conducted in legacy multipactor analysis methods [3], [4], [37].
By solving (1), one can derive the analytical expressions for the position y i and y-axis velocity v y,i of the ith particle being y 0i and v y,0i the y-axis position and initial velocity just after the time instant t 0i of the last impact.
After each particle collision, the number of electrons corresponding to the particle is updated according to , where t corresponds to the time resolution of the simulation [20].The distribution of secondary electron departure kinetic energy, denoted as W s , is modeled using a Rayleigh probability density function (pdf) where W g represents the characteristic parameter of the Rayleigh distribution, which has been fixed to a value of W g = 3 eV.With regard to the release angle (with respect to the normal of the surface), φ, it is assumed to be a random variable following a cosine law distribution [35], [38], which is calculated as follows: where x is a random variable distributed uniformly between 0 and 1.
By considering the kinetic energy and the release angle, the initial velocity components v x,0i and v y,0i after the impact are evaluated.Note that for the x-axis, a constant-velocity movement is obtained as the electric field of the TEM mode does not have x-axis component.However, the variation of the SEY with the incident angle is taken into account using the usual Vaughan model [35], [39].
As the secondary electron emission is governed by a stochastic process, related to the random variables of energy and angle, a Monte Carlo simulation is usually implemented in order to attain high statistical accuracy in the multipactor threshold estimation.Therefore, the algorithm follows the independent trajectories as expressed in (2) for a sufficiently large number n of effective electrons.This is indeed a timeconsuming task, due to the large number of particles that must be tracked to obtain precise predictions.A multipactor discharge is thus identified if after k RF cycles, the equivalent electron population, calculated as the sum of N i (t 0 + kT ), is greater than the initial one [i.e., the sum of N i (t 0 )].In that case, an exponential growth can be expected due to the secondary emission process.Nevertheless, if the population decreases, the absorption process is assumed to dominate and consequently, no multipactor discharge is predicted to occur.

III. GENETIC METHOD IMPLEMENTATION
When a Monte Carlo simulation detects an electron population growth, typically only a few particles contribute significantly to this growth, while the rest exhibit a negligible mass relative to these prominent contributors.This occurs because the parameters of the remaining particles do not meet the resonance conditions required for the ignition of a multipactor discharge.The genetic method proposed in this work offers the capability to modify specific parameters of each particle along the simulation.By doing so, an accelerated convergence can be achieved by using only a limited number of particles, thus significantly reducing the computational effort over a traditional Monte Carlo multipactor simulation.The proposed approach allows to finely control the behavior of individual effective electrons during the simulation process, enhancing the overall computational performance.
Let S be the mathematical space formed by the vectors containing the information of the different n particles, which can also be referred as genetic information where R * ,n + represents the n-dimensional space of nonzero positive real numbers, whereas R n + depicts all positive real numbers.
Within this space S, there exists a subspace P (P ⊆ S) in which specific combinations of random parameters lead to the same convergent results than the ones obtained by using a sufficiently large number of tracked particles.The understanding of the process by which the subspace P is obtained, and the characterization of its properties, is crucial to understand the stochastic dynamics.Therefore, the main idea of this method is to identify, during the simulation, the most Fig. 2. Representation of space S containing subspace P. Initially, particles are distributed throughout the entire space.After applying the algorithm, all these particles eventually converge and end up within the subspace P.
critical particles that lead to the same multipactor prediction than a simulation in the full space, and then replicate its parameters for other tracked particles, as illustrated in Fig. 2. Thus, the simulation converges rapidly to the desired behavior as can be seen in the expression where g n and α i ( f, d, V 0 ) represents the temporal evolution of the electron population composed by n particles, and the population growth factor rate [16], [40] for a given f × d product and V 0 , respectively.Since P is a subset of S, g p is evaluated in a smaller search space and, therefore, has fewer possible parameter combinations compared to f p : S × R → R.
The optimized threshold prediction algorithm comprises the following three main tasks.
1) Start particle simulation to identify the most multipactorprone particle, that is, the one with the highest secondary emission yield value associated in the last set of N imp impacts (N imp = 100 in this work), and obtain their genetic information.2) Assign its characteristics (the departure kinetic energy W s and release angle φ after the impact providing a highest SEY value for the last set of N imp impacts, and the time instant of the last impact) to the less multipactorprone ones.A proportion β (out of one) of all particles are modified, thereby their genetic information is similar to the one from the most multipactor-prone particle.So, β is the first tuning parameter of the algorithm.
3) The assignment process is characterized by a parameter κ, which controls how similar the genetic information of the less multipactor-prone particles (proportion β) is with respect to the genetic information of the most multipactor-prone one.Therefore, κ is the second tuning or adjustment parameter.The search for the subspace P is initiated by populating the entire space S with seeding particles composed initially of only one electron.As the simulation progresses, some particles may approach the vicinity of the subspace P. The particles with largest effective mass m N i are the most multipactor-prone ones.The genetic information of these particles is stored, so that, the genetic information of the "worst" one is selectively replicated to the proportion β of effective electrons with the lowest mass, after a simulation interval T s (corresponding to N imp particle impacts) involving several RF cycles.
It is worth pointing out that this genetic information is not replicated equally.Instead, a uniform statistical distribution Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
within a range centered around the extracted value (µ) of the genetic characteristic, but constrained by a κ proportional value, is used (i.e., the value is in the interval [µ − κµ, µ + κµ]).This precise calibration is crucial, as it prevents immediate convergence to a local minimum, and facilitates a balanced transfer of information among the particles.
Importantly, the total effective mass (and charge) of the population is redistributed evenly among the particles, so that each particle in the system carries an equal share of the total effective mass at the beginning of each new simulation interval.This equitable distribution ensures that the impact of each individual particle on the overall behavior of the system is balanced, leading to a more stable and accurate representation of the collective behavior.Additionally, this equal sharing of effective mass plays a crucial role in achieving statistical convergence more rapidly, as it prevents certain particles from dominating the simulation while others remain underrepresented.As a result, the whole process can progress efficiently with a more uniform contribution from all particles, facilitating a better representation of the system dynamics and enhancing the accuracy of the predictions.
By carefully controlling the range κ and proportion β of the genetic information replicated, we ensure a more controlled exploration of the parameter space, enabling a faster identification of the critical parameter combinations that lead to the long-term evolution of the particle population.As a result, the number of tracked particles required to obtain accurate multipactor threshold predictions can be reduced, thus leading to an improvement in the computational efficiency of multipactor simulations with respect to the traditional Monte Carlo approach.

A. Convergence Analysis
In the current section, a convergence study for the proposed method is carried out with two different goals.The first one is to confirm that the new approach converges to the right multipactor threshold.The second goal is to compare the computational effort with respect to traditional Monte Carlo simulations.
The proposed multipactor threshold determination algorithm operates by considering a range of voltages (starting from an initial value), then checking if these fulfill the multipactor criterion.For the initial voltage value, the curve representing the time evolution of the electron population is obtained, and the associated growth factor is defined as the natural logarithm of the ratio between the final and initial overall particle mass [40].If the electron population count has increased (positive growth factor), a discharge is identified.Otherwise, the voltage is increased by a factor of (2) 1/2 (i.e., a scale factor of 2 in terms of power) until an increase in the particle population is observed.Next, at each step, the interval between the highest voltage without discharge and the lowest voltage where multipactor is predicted is halved (in logarithmic terms), by simulating the case corresponding to the geometric mean of the voltages of the interval extremes.The iterative method terminates when the difference between the current voltage value and the previous one is below 0.1 dB (precision), taking the geometric mean of the last two voltage values as the multipactor voltage threshold (V th ).This process is repeated k = 10 times for the same number n of tracked particles, obtaining the voltage threshold as the average threshold value of these k iterations.As a metric to measure the convergence, the root mean squared error (RMSE) is also used, as it quantifies the quadratic mean of the differences between the reference threshold value V th,ref and the predicted ones V th,i .The expression of RMSE is given by the following equation: First, a convergence analysis is performed at a frequency of 2.5 GHz with a gap of 1 mm, since for this f × d value (2.5 GHz•mm) a larger variability is expected due to the 1/2-to 3/2-mode transition in the multipactor susceptibility chart for ECSS aluminum surfaces.The reference threshold voltage has been obtained by a SPARK3D simulation with 1000 seeding electrons (one electron per particle) using the EM fields of an ideal parallel-plate waveguide with lateral magnetic walls, leading to V th,ref = 128 V.The effect of the proportion of particles with replicated information has been analyzed, considering the β values of 10%, 25% and 50% in the following set of results.The interval search proportion κ has been set to 10%.The check time T s considered by the genetic algorithm to perform a redistribution of the particle masses and charges is equal to the time required for N imp = 100 overall impacts.It is worth pointing out that this check time shows a dependence with the multipactor order due to the particle flight time.Lastly, after a convergence analysis, a very conservative resolution of 20 000 points per cycle has been used throughout the simulations (i.e., a time resolution lower than 1 ps).All the simulations are carried out in a computer server with an Intel Core i9-9900k at 3.6 GHz, with eight cores and 128 GB RAM.
The results from Fig. 3(a)-(d) clearly point out that the proposed genetic algorithm requires a much lower number of tracked particles than a traditional Monte Carlo analysis to obtain convergent results.For a proportion of altered particles β of 10%, around 30 particles are required with an average computational effort of 45.49 s, whereas the Monte Carlo method still exhibits a relevant variability [and therefore RMSE error, according to (8)] with 200 tracked particles and an average CPU time of 211.29 s, although both methods tend to converge to the reference voltage computed by SPARK3D.As a consequence, a significant reduction in the computational effort for the same statistical variability measured in terms of RMSE [see Fig. 4(a) and (b)] is achieved.
It was observed, however, that the update time T s and the proportion of particles with replicated information β play a key role in the outcome of the genetic algorithm [as can be seen in Fig. 3(d)].In case that T s is not large enough to properly identify the most multipactor-prone particle (wrongly taken, for instance, a particle with a lower order multipactor resonance which suffers more impacts per time unit), a high value of β would statistically enhance the mistake since a too large proportion of particles inherits the characteristics of the   "worst" one identified so far.This may result in an increase in the multipactor threshold, similar to a simulated annealing algorithm cooling too fast [41].This was the situation in Fig. 3(d) when the number of particles increases, resulting in a lower effective T s value.Note that T s should also be linked to the electron flight time, and therefore with the multipactor order.As pointed out in Table I and in Figs.4(a) and 5, a value of β around 10% (considering N imp = 100 impacts) seems to provide the best performance, as the method converges to the right multipactor threshold with the lowest number of tracked particles.Finally, it is worth mentioning that the CPU times in Fig. 4(b) are the ones corresponding to the complete set of simulations with different voltage levels required to obtain the final multipactor threshold.Anyway, the benefit of using a reduced number of tracked particles is evident, when considering that the classic Monte Carlo method would require much more than 60 particles to reduce the RMSE error and obtain a convergent multipactor threshold.

B. Algorithm Effectiveness
To further validate the proposed algorithm, the impact energy pdf for the genetic technique with a reduced number of Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.particles (n = 20) is compared with the one obtained with the traditional Monte-Carlo method using 1000 tracked particles.Besides, the results from CST Particle Studio simulations are included to compare the impact energy distribution with commercial particle-in-cell (PIC) codes.The CST simulation has been performed in the same reference structure with 1000 seeding particles.A Maxwellian particle source is used, whereas one electron per particle is considered.Moreover, it is worth pointing out that the modified Vaughan's SEY curve (aluminum ECSS [36]) is imported, and a temperature of 3 eV is used for the energy distribution of the secondaries [16].A 2.5 GHz sinusoidal signal (with a duration of 100 ns) is used to excite the structure.Fig. 6(a) and (b) shows very similar results comparing the traditional Monte Carlo simulation and the genetic one with the CST simulation (such results are nearly the same for different β values).Therefore, the simulation with a limited number of tracked particles corresponding to the genetic algorithm is fully representative of a traditional Monte Carlo analysis involving a larger number of particles.
Additionally, a susceptibility multipactor chart is generated to assess the effectiveness of the genetic algorithm.Fig. 7 the mean effective for simulation runs of 100 RF cycles (with 20 000 points per cycle), for resolution of 200 × 200 simulations the voltage and frequency-gap product domains.The threshold boundary (corresponding a mean effective SEY of 1) is denoted with a black is compared with the values obtained using the ECSS Multipactor Tool v.2 and SPARK3D simulators.An excellent agreement is indeed observed.Only for the low frequency-gap products small discrepancies arise, as the proposed algorithm provides a more conservative threshold (namely, the genetic method identifies a set of particle conditions allowing a sustained increase in the electron population at a lower voltage).Note, however, that this is indeed the most sensitive region for multipactor threshold prediction, and that both the ECSS Multipactor tool and SPARK3D are founded on nonstationary statistical theory.It must be recalled that the total time required to obtain the results in Fig. 7 (involving 40 000 particle simulations) was close to 29 h.Given these extensive computational demands, it is fully unfeasible to replicate this graph using the traditional Monte Carlo method with 1000 seeding particles.

V. CONCLUSION
In this article, a technique for speeding up multipactor particle simulations has been proposed.The use of an adaptive genetic modification algorithm in Monte Carlo simulations allows to reduce the number of particles to be tracked, whilst maintaining the same statistical performance.Thus, a reduction of about one order of magnitude in the minimum number of tracked particles required to attain a reduced threshold estimation error is achieved.This, in turn, leads to a relevant improvement in terms of computational cost.The algorithm has been successfully applied for a 2-D parallel plate geometry, which is the reference case typically used in multipactor standards.By optimizing the method, convergence is sped up and computational effort is reduced, with respect to traditional multipactor Monte Carlo simulations, through adaptive particle properties modification.
Thanks to its general scope, this technique could be generalized to more involved geometries such as those found in RF components, alleviating the design effort of complex devices with high-power requirements.However, this extension is not straightforward, as the electron loss rate (due to inhomogeneous and fringing fields) must be preserved to avoid an artificial reduction in the multipactor threshold.

Fig. 3 .
Fig.3.Evolution of the voltage threshold prediction (average value in V ) in terms of the number of tracked particles using (a) traditional Monte Carlo approach and the proposed genetic algorithm for β values of (b) 10%, (c) 25%, and (d) 50%.

Fig. 5 .
Fig. 5.Boxplot comparison for the threshold voltage between the traditional Monte Carlo approach and genetic algorithm with β values of 10% and 25%.

Fig. 6 .
Fig. 6.Comparison of impact energy pdfs for a voltage amplitude of (a) 100 V and (b) 160 V. Solid black and blue lines represent the energy distribution of traditional Monte Carlo and genetic algorithms, respectively, while dashed blue lines depict the energy distribution obtained from CST simulations.

TABLE I MULTIPACTOR
VOLTAGE THRESHOLD AVERAGE AND RMSE USING THE TRADITIONAL MONTE CARLO APPROACH IN COMPARISON WITH THE PROPOSED GENETIC ALGORITHM FOR β VALUES OF 10%, 25% AND 50% AND n TRACKED PARTICLES.A κ VALUE OF 10% IS CONSIDERED