Axisymmetric Charge-Conservative Electromagnetic Particle Simulation Algorithm on Unstructured Grids: Application to Microwave Vacuum Electronic Devices

Axisymmetric Charge-Conservative Electromagnetic Particle Simulation Algorithm on Unstructured Grids:
Application to Microwave Vacuum Electronic Devices

Dong-Yeop Na Yuri A. Omelchenko Haksu Moon Ben-Hur V. Borges Fernando L. Teixeira ElectroScience Laboratory and Department of Electrical and Computer Engineering, The Ohio State University, Columbus, OH 43212, USA Trinum Research Inc., San Diego, CA 92126, USA Intel Corporation, Hillsboro, OR 97124, USA Department of Electrical and Computer Engineering, University of Sao Paulo, Sao Carlos, SP 13566-590, Brazil

We present a charge-conservative electromagnetic particle-in-cell (EM-PIC) algorithm optimized for the analysis of vacuum electronic devices (VEDs) with cylindrical symmetry (axisymmetry). We exploit the axisymmetry present in the device geometry, fields, and sources to reduce the dimensionality of the problem from 3D to 2D. Further, we employ ‘transformation optics’ principles to map the original problem in polar coordinates with metric tensor to an equivalent problem on a Cartesian metric tensor with an effective (artificial) inhomogeneous medium introduced. The resulting problem in the meridian () plane is discretized using an unstructured 2D mesh considering -polarized fields. Electromagnetic field and source (node-based charges and edge-based currents) variables are expressed as differential forms of various degrees, and discretized using Whitney forms. Using leapfrog time integration, we obtain a mixed finite-element time-domain scheme for the full-discrete Maxwell’s equations. We achieve a local and explicit time update for the field equations by employing the sparse approximate inverse (SPAI) algorithm. Interpolating field values to particles’ positions for solving Newton-Lorentz equations of motion is also done via Whitney forms. Particles are advanced using the Boris algorithm with relativistic correction. A recently introduced charge-conserving scatter scheme tailored for 2D unstructured grids is used in the scatter step. The algorithm is validated considering cylindrical cavity and space-charge-limited cylindrical diode problems. We use the algorithm to investigate the physical performance of VEDs designed to harness particle bunching effects arising from the coherent (resonance) Cerenkov electron beam interactions within micro-machined slow wave structures.

Electromagnetic simulation; particle-in-cell; unstructured meshes; vacuum electronic devices; backward-wave oscillator; Cerenkov radiation.

1 Introduction

Historically the need for high-power electromagnetic (EM) radiation sources in the gigahertz and terahertz frequency ranges has triggered significant technical advances in vacuum electronic devices (VED) gold1997review ; cairns1997generation ; booske2011vacuum ; booske2008plasma ; li2013experimental , such as the gyrotron, free electron Laser, and traveling wave tube (TWT). These devices serve as a basis for a variety of applications in radar and communications systems, plasma heating for fusion, and radio-frequency (RF) accelerators gaponov1994applications ; schamiloglu2004high .

Amplification of RF signals is usually obtained by exploiting resonance Cerenkov interactions between an electron beam and the modal field supported by a slow-wave structure (SWS) cairns1997generation ; shiffler1990high ; johnson1955backward ; gunin1998relativistic ; case1984space . SWSs are often made by imposing periodic ripples on the conducting wall of cylindrically symmetric waveguides, as illustrated in Fig. (a)a, so that the phase velocity of the modal field becomes slower than the speed of light in vacuum due to Bragg scattering chew1995waves . According to the dispersion relations associated with the geometry of SWSs, resultant Cerenkov interactions can amplify forward or backward waves111Along the direction of the group velocity w. r. t. the beam velocity.. Similarly to plasma instabilities, the evolution of forward and backward waves can be characterized by convective instabilities that grow over time while traveling away from the location of initial disturbance and absolute instabilities that propagate a local initial disturbance throughout the whole device volume cairns1997generation . Traveling-wave tube amplifiers (TWTA) and backward-wave oscillators (BWO) are two practical examples utilizing convective and absolute instabilities, respectively.

Recent studies have shown that a particular SWS geometry may significantly enhance the system performance of TWTs. For example, nonuniform (locally periodic) ripples used in BWOs may improve mode conversion efficiency moreland1994efficiency ; chipengo2015novel , and tapering ripples may reduce reflections at the output of TWTA and prevent internal oscillations shiffler1990high ; shiffler1991high . More importantly, smooth device edges are preferred for high-output power applications in order to mitigate pulse shortening, which is a major bottleneck for increasing output powers beyond the gigawatts range agee1998evolution ; korovin2000pulsewidth . This is because extremely strong field singularities, which accumulate on the sharp edges, may create interfering plasmas that terminate the output signal at an earlier time. Sinusoidally corrugated slow wave structures (SCSWS) have been increasingly adopted in many modern high-power BWO systems to combat this problem vlasov2000overmoded . In addition, a variety of micro-machining fabrication techniques have been developed to enable better device performances by using much tighter tolerances. These technological advances have allowed the production of devices operating at higher frequencies, including the THz regime.

Figure 1: Schematics of two examples of axisymmetric vacuum electronic devices. (a) Backward-wave oscillator producing bunching effects on an electron beam. Wall ripples are designed to support slow-wave modes in the device. (b) Space-charge-limited cylindrical vacuum diode.

Computational experiments for VEDs employ electromagnetic particle-in-cell (EM-PIC) algorithms birdsall2004plasma ; hockney1988computer ; grigoryev2002numerical ; candel2010parallel , which numerically solve the Maxwell-Vlasov equations describing weakly coupled (collision-less) systems, where the collective behavior of charged particles prevails over their binary collisions birdsall2004plasma ; hockney1988computer ; dawson1983particle . A typical PIC algorithm tracks the temporal evolution of macro-particles seeded in a coarse-grained six-dimensional (6D) phase space222That is, a finite-size ensemble of physical particles with positions and momenta.. A typical PIC algorithm consists of four basic steps, viz. the field solver, the field gather, the particle push, and the particle charge/current scatter, which are repeated at every time iteration. This provides a self-consistent update of particles and field states in time.

As a field solver, most EM-PIC simulations employ the celebrated Yee’s finite-difference time-domain (FDTD) method for regular grids, due to its simplicity. There are a plethora of FDTD-based EM-PIC codes such as UNIPIC, MAGIC, TWOQUICK, KARAT, VORPAL, and others wang2009unipic ; wang2010three ; nieter2009application . However, the relatively poor grid-dispersion properties of this algorithm taflove2000computational causes spurious numerical Cerenkov radiation greenwood2004elimination . Moreover, in complex geometries such as those of modern VEDs, “staircase” (step-cell) effects present a critical challenge. Using FDTD for an accurate analysis of geometrically complex devices, which typically have curved boundaries or very fine geometrical features, may require excessive mesh refinement and therefore result in a waste of computational resources. Many studies have been done to mitigate the staircasing errors in finite-difference (FD) methods, in particular through using conformal FD discretizations meierbachtol2015conformal ; wang2016conformal .

On the other hand, the finite-element time-domain (FETD) method teixeira2008time ; lee1997time fundamentally eliminates the undesired staircase effects since it is naturally based on unstructured (irregular) grids, which can more easily be made conformal to complex geometries and can be augmented by powerful mesh refinement algorithms. Unfortunately, conventional FETD-based PIC algorithms have historically faced numerical challenges that result from a lack of exact charge conservation on unstructured grids. This gives rise to the accumulation of spurious charges which must be removed by applying costly a posteriori corrections eastwood1991virtual ; marder1987method . In addition, implicit time updates used in conventional FETD require repeated linear solves at each time-step teixeira2008time . Recently, a novel charge-conserving scatter scheme for unstructured grids, inspired by differential-geometric ideas and the exterior calculus of differential forms burke1985 ; flanders1989 ; teixeira1999lattice ; kotiuga2004 ; he2007differential ; teixeira2014lattice , has been proposed in moon2015exact . Other charge-conservative EM-PIC algorithms for unstructured grids were also developed under similar tenets in squire2012geometric ; pinto2014charge ; pinto2016divergence . In addition, a charge-conserving EM-PIC algorithm with explicit time-update that is both local (i.e. preserves sparsity) and obviates the need for linear solvers at each time step has been described in he2006sparse ; kim2011parallel ; na2016local based on the sparse approximate inverse (SPAI) of the discrete Hodge operator (the finite-element “mass” matrix).

These recent advances have made possible the present work, which is motivated by the demand to accurately capture realistic physics of beam-SWS interactions in complex geometry devices. In this paper we present a charge-conservative EM-PIC algorithm based on unstructured grids and optimized for the analysis and design of axisymmetric VEDs. Since conventional SWSs are cylindrically axisymmetric (invariant along ), SWS studies can be best done with algorithms that explore this symmetry, so that significant computational resources can be saved and the algorithm may be feasibly implemented as a forward engine in a design loop wang2016conformal . We show that under the assumption of cylindrical symmetry of fields and sources one can reduce the original 3D geometry to a 2.5D setup by introducing an artificial inhomogeneous medium333In a manner akin to the “transformation optics” (TO) technique teixeira1999lattice ; teixeira1999differential ; pendry2006controlling ; chen2010nature ., considering -polarized fields in the meridian () plane. The spatial discretization of EM variables, which are represented by differential forms of various degrees, is achieved through applying discrete Whitney forms on an unstructured mesh in the meridian plane. Using leapfrog time integration, we obtain space- and time-discretized Maxwell’s equations, a so-called mixed finite-element time-domain (FETD) scheme, for the field solver. A local and explicit time-update requiring no repeated linear solver is made possible by employing a SPAI algorithm na2016local . The interpolation of field values to particles’ positions for solving the Newton-Lorentz equations of motion is also done via Whitney forms. Relativistic particles are accelerated and pushed in space with a corrected Boris algorithm verboncoeur2005particle ; vay2008simulation . In the scatter step, we adapt an exact charge-conserving scheme developed for Cartesian unstructured grids moon2015exact . We validate our algorithm using analytical results for a cylindrical cavity and previously obtained results for a space-charge-limited (SCL) vacuum diode (see Fig. (b)b). We include a micro-machined SWS-based BWO example, designed to harness particle bunching effects from coherent Cerenkov beam-wave interaction, to demonstrate the advantages of utilizing unstructured grids without staircasing error to predict the device performance.

2 Spatial dimensionality reduction

In this section we describe a numerical model for 3D VEDs with cylindrical axisymmetry based on an equivalent 2D model discretized on an unstructured (irregular) grid in the meridian plane.

2.1 Exterior calculus representation of Maxwell’s equations

We represent Maxwell’s equations using the exterior calculus of differential forms flanders1989 ; kotiuga2004 ; teixeira2014lattice ; squire2012geometric ; FEMSTER2005 as


where and are 1-forms for the electric and magnetic field intensity, and are 2-forms for the electric and magnetic flux density, is 2-form for the electric current density, is 3-form for the electric charge density, and operator is the exterior derivative encompassing conventional curl and divergence operators warnick1997teaching ; sen2000geometric ; bossavit1988whitney ; teixeira2013differential . These 1-, 2-, and 3-forms can be expressed using a set of non-orthonormal-basis in a cylindrical coordinate system teixeira1999differential ; warnick1997teaching . For instance, (1-form) is expressed as ; then, its vector proxy can be written as . Similarly, (2-form) is given by (where is the exterior or wedge product burke1985 ; flanders1989 ; sen2000geometric ) and its vector proxy takes the form of  warnick1997teaching .

2.2 Cylindrical axisymmetry constraints

Consider a charged ring with constant density along azimuth that travels inside of a cylindrically axisymmetric tube bounded by a perfect electric conductor (PEC) with a radial boundary profile of where is the computational domain and is the wall radius which only depends on , as shown in Fig. (a)a.

Figure 2: A charged ring travels inside an axisymmetric object bounded by PEC: (a) a 3D view, (b) the meridian plane.

Cylindrical axisymmetry, used here, implies that there is no variation along () in the device geometry, fields, and sources. It should be noted that axisymmetric sources in the meridian plane are represented as charge rings (see Fig. (b)b). There exist two useful constraints that simplify the original 3D problem: (i) elimination of the term in the exterior derivative , viz. and (ii) retainment of only transverse magnetic (TM) eigenmodes with 444The index is used here to denote azimuthal harmonics.. The first constraint enables the same calculus in the meridian plane as in the 2D Cartesian coordinate system with the cylindrical metric factor embedded into the constitutive relations, as discussed below. The second constraint simplifies expressions for fields and sources as


From equations of (5), (6), (7), (8), it is straightforward to show that a 3D problem with cylindrical axisymmetry can be represented by a 2D problem describing -polarized fields on the meridian plane (see also Fig. 3).

2.3 Modified Hodge star operator

The Hodge star operator maps -forms into -forms in -dimensional space555The Hodge star operator can be understood geometrically as yielding the orthogonal complement of a given differential form to the volume form in -space. In the 3D case for example, = , = , = = and so forth warnick1997teaching . For an arbitrary -form in -space, we have , where is the volume element in -space. As such, provides the ( or“energy”) element norm of . he2007differential ; teixeira2014lattice ; warnick1997teaching ; teixeira2013differential ; tarhasaari1999some ; gillette2011dual . In our case, we have


The Hodge operators incorporate the metrical properties of the system, which in the cylindrical case are expressed by a metric tensor . For the magnetic field and flux density, substituting (8) into the left-hand side term of (9) gives


and substituting (6) into the right-hand side term of (9) yields


Note that only acts on the differentials such as , , and . By comparing (11) and (12) and introducing an artificial magnetic permeability, , we can extract the radial factor from the Hodge star operator. As a result, we can reuse simple Cartesian space calculus with metric tensor and a Hodge operator devoid of additional metric factors. For the electric field and flux density, substituting (7) into the left-hand side term of (10) yields


and substituting (5) into the right-hand side term of (10) gives


In a similar fashion, an artificial electrical permittivity takes the form of . Essentially, an original 3D problem with cylindrical axisymmetry is replaced by an equivalent 2D problem with -polarized fields immersed in Cartesian space with an inhomogeneous medium with artificial permittivity and permeability, as illustrated schematically in Fig. 3.

Figure 3: The original problem shown in Fig. 2 is replaced by an equivalent 2D problem in the meridian plane as depicted above, which considers -polarized EM fields on Cartesian space with an artificial inhomogeneous medium. The variable coloring serves to stress the dependency of the artificial medium parameters on .

As a result, the present field-solver borrowing the concept of TO is free from the axial singularity associated with the factor present in the cylindrical nabla operator.

3 An axisymmetric EM-PIC algorithm

EM-PIC algorithms primarily operate in time domain. Their computational cycle typically consists of the following steps: (1) EM field update (solution of Maxwell’s equations), (2) field gather (interpolating field values to the particle positions), (3) particle update (solution of the Newton-Lorentz equations of motion), and (4) particle scatter (computing grid currents and charges from the kinetic information stored in charged particles).

3.1 Field solver: Mixed FETD with local and explicit time update

In the field update, space- and time-discretized Maxwell’s equations are solved for all discrete EM fields in one time-step. For the 2D equivalent problem discussed in Sec. 2, we can define , , , and as a 1-, 2-, 1-, and 0-forms, respectively, and and as 1- and 0-forms 666Since and are 1- and 2-forms, respectively.. It should be stressed that the operator used here for the equivalent problem is defined in the Cartesian 2D space and hence devoid of the radial factor (which enters instead in the definition of the artificial medium properties). Because and are ordinary differential forms with internal orientations associated with the primal lattice teixeira2014lattice ; teixeira2013differential , we expand them by using Whitney 1- and 2-forms (cochains) residing in the primal mesh, as follows:


where is the number of -cells for (number of nodes, edges, facets, respectively, on the grid), is the Whitney -form associated with the -th -cell on the grid, and and are the (time-evolving) degrees of freedoms (DoFs) for and , respectively white2015mixed . On the other hand, , , , and are twisted differential forms with external orientations, associated with the dual lattice teixeira1999lattice ; teixeira2014lattice ; sen2000geometric . Fig. 4 shows schematically how the DoFs are associated to the primal and dual meshes elements in our 2D setting.

Figure 4: The (2+1) setup for fields and sources defined in (a) primal and (b) dual meshes.

To avoid the need for working with the dual mesh, we employ the Hodge duals and which can be expanded using Whitney 1- and 0-forms in the primal mesh, as follows:


where and are Dos of -th edge and node for and , respectively. Applying a pairing operation and generalized Stokes’ theorem to Faraday’s law in the primal mesh and Ampere’s law in the dual mesh teixeira1999lattice ; teixeira2014lattice ; he2005degrees , the spatially discretized Maxwell’s equations (see also A) are given by


with , , , and where and are incidence matrices encoding the discrete representation of the exterior derivative (or, equivalently, the curl operator distilled from the metric) on the primal and dual meshes, respectively teixeira1999lattice ; teixeira2014lattice ; hughes1981lagrangian ; guth1980existence . Rows of have a association to faces of the primal mesh and columns to edges. On a primal triangular mesh, each row of has three nonzero elements corresponding to the three edges of the triangular cell associated to that row (with entries or depending on their relative orientation). All the other elements of that row (corresponding to edges that are not part of the boundary of the triangular cell boundary) are zero. A similar relationship exists for so that both of these matrices take their elements in the set of as a consequence of their metric-free character. It can be shown that = teixeira1999lattice . In (19) and (20), and are discrete Hodge operators (see B), which are symmetric positive-definite matrices obtained by the Galerkin method kim2011parallel ; tarhasaari1999some ; gillette2011dual ; he2006geometric . Their elements are expressed as the area (volume in 3D) integrals


where is the intersection of the Whitney functions spatial support, respectively. Vector proxies of (21) and (22) are provided in B. Note that the discrete Hodge matrix for the electric current density becomes the identity matrix for the barycentric dual mesh (see C). We apply leapfrog time integration for time discretization with the electric field, flux density and charge density defined at integer time-steps and the magnetic flux density, field and current density defined at half integer time-steps . The resulting space- and time-discretized Maxwell’s equations form the so-called mixed finite-element time-domain (FETD) scheme moon2015exact ; kim2011parallel ; na2016local ; donderici2008mixed :


where superscripts denote time-step and is the time-step increment. Local explicit time-update is achieved using the SPAI approach to obtain a representation for and obviate the need for a (expensive) linear-solver at each time-step kim2011parallel ; na2016local . Since the discrete equations (23) and (24) are basically sparse linear systems, matrix multiplications with compressed sparse row format are adopted to reduce computational costs.

3.2 Field gather

The EM fields must be determined at the positions of -species particles in order to solve the (Newton-Lorentz) equations of motion for the latter. This is again done through using Whitney forms, which recover the necessary field values from their grid-defined representations as moon2015exact


where is a vector proxy of the Whitney -form for -th grid element and is the particle position vector.

A key component in the operation of BWOs is the beam focusing system (BFS) which axially confines an electron beam within an extremely small gyroradius in the polar () plane to avoid the bombardment of the SWS by the beam electrons. Beam focusing is usually achieved by applying a strong static axial magnetic field in the SWS region. Because a strong confinement of the electron beam, azimuthal particle bunching effects can be ignored.

3.3 Relativistic particle pusher

The positions and velocities of -species particles are updated by solving the Newton-Lorentz equations of motion as na2016local .


where is a velocity vector, is a particle mass, is a particle charge, and , and . The leapfrog time-discretized Newton-Lorentz equations of motion can be written as


Here, the position and velocity vectors of particles are defined at integer and half-integer time-steps, respectively. In order to solve the implicit equation (30), we adopt the Boris algorithm with a correction verboncoeur2005particle ; vay2008simulation .

Figure 5: The grid charge evolution at node during a single particle push. The charge states at time-steps of , , , and are illustrated in (a), (b), (c), and (d), respectively. Note that a colored triangular area, formed by the particle position, nodes , , corresponds to the grid charge of node marked as a blue-colored circle.

3.4 Charge-conserving particle scatter

The particle scatter transfers the kinetic information from the charged particles to the grid for its subsequent use by the field solver. As described before, the discrete (Hodge dual) charge and current densities are defined at nodal and edge elements of the primal mesh, respectively. As mentioned in Sec. 1, conventional interpolation techniques for irregular grids may violate charge conservation and cause spurious (nonphysical) charge deposition on the grid pinto2014charge ; na2016local . A charge-conserving scatter algorithm has been recently proposed moon2015exact based on the use of discrete Whitney forms as a consistent interpolation (“particle shape”) factor. For piecewise-linear particle trajectories, charges and currents associated with 0- and 1-forms are defined as follows moon2015exact :


where represents the -th grid element, and and denote indices of lower dimensional manifolds of -th element777For example, if -th element is an edge, and indicate the two nodes associated with the said edge.. The functions and are scalar and vector proxies of Whitney 0- and 1-forms, respectively, and is a barycentric coordinate for an -species particle at time-step , respectively. One can find the derivations of Eqs. (31) and (32) in more detail in moon2015exact .

Figure 6: The grid current evolution at edge is indicated by the thickness of the green line. Similar to the nodal charge, its equivalent amount is illustrated by a colored triangular area.

Fig. 5 and Fig. 6 provide the geometrical depiction on the node-based charges and edge-based currents, respectively, obtained by Eqs. (31) and (32). For example, assume a charged particle movement from time-step, to time-step, . As shown in Fig. 5, the grid charge (blue-colored circle) at node is calculated based on (31). The associated charge value is represented by the colored triangular area formed by the particle position and nodes and . Likewise, as shown in Fig. 6, the respective grid current at edge (32) is given by a gradient-colored triangular area formed by the particle trajectory and the node that does not belong to that edge.

3.5 Verification of charge conservation

3.5.1 Discrete continuity equation (DCE)

In order to verify charge conservation, we examine the discrete charge continuity equation (DCE) given by moon2015exact ; na2016local


where is an incidence matrix encoding the discrete representation of the exterior derivative in the dual mesh or, equivalently, the divergence operator distilled from metric information. Incidence matrices take their elements in the set of , see also the Appendices. Note that Hodge matrices for the discrete Hodge duals of charges and currents ( and ) become identity matrices for barycentric dual lattices. Therefore, and .

Figure 7: Schematics for verifying charge conservation. The charge rate at node exactly matches the net current flowing from/to node .

The first term in the left-hand side of (33) stands for the rate of change in the charge during one time-step. The second term in the left-hand side of (33) represents the net current flowing from/to the nodes. Consider a charged particle that moves from time-step of to time-step in a given triangular element (see Fig. 7). Then, the rate of change of the nodal charge at node is


The net current flowing from/to node is a sum of the edge currents at edges and as follows:


Thus, the rate of change of the nodal charge at node is equal to the net current flow. From the geometrical perspective, this is depicted in Fig. 7.

3.5.2 Discrete Gauss’ law (DGL)

Premultiplying both sides of (33) by , where is the incidence matrix representing the discrete divergence operator on the dual grid and is an exact inverse of obtained by using SPAI algorithm, and utilizing the identity  teixeira1999lattice ; teixeira2014lattice yields


which can be rearranged as


and, using (33), rewritten as


Eq. (38) implies that residuals of the discrete Gauss’ law (DGL) at any two successive time steps remain the same, in other words


and by induction,


for all . This means that if consistent initial conditions are set, then the DGL remains valid for all time steps.

4 Validation

In this section, we provide validation examples. First, we consider a metallic cylindrical cavity and compare the resonant frequencies of the cavity modes obtained by the present field solver with the exact (analytic) results. Second, we model a space-charge-limited cylindrical diode with a finite-length emitter and compare the maximum injection currents for divergent and convergent electron beam flows against previously published results.

4.1 Metallic hollow cylindrical cavity

Assume a hollow cylindrical cavity with metallic walls, radius m and height m (see Fig. 8). Two magnetic point sources and excited by Gaussian-modulated pulses of broadband spectrum are placed at arbitrary positions given as


with ns, ns, and rad/s. The source locations are denoted by Tx in Fig. 8. These sources excite resonant modes inside the cavity. The 2D meridian plane of the cylindrical cavity is discretized by an unstructured grid with nodes, edges , and faces. The three lateral metallic boundaries are assumed as perfect electric conductors (PEC). The remaining boundary is the -axis (axisymmetric boundary). The time step interval is chosen as ps, and the simulation runs over a total of time steps. Fig. 8 shows a snapshot for electric field distribution in the cavity at s. The RGB colormap and the white arrows indicate magnitudes and vectors of the electric fields, respectively.

Figure 8: Snapshots for electric field distribution at 2 s. Note that RGB colors and white arrows indicate magnitudes and vectors of the electric fields, respectively.

The time signal are detected at (Rx in Fig. 8) and a Fourier analysis is performed to obtain a spectrum of the signal. The resulting spectrum shows the resonant cavity modes from MHz to GHz in Fig. 9, where blue-solid lines are axisymmetric FETD field-solver results and red-dashed lines are analytic results. An excellent agreement can be observed. In addition, Table. 1 shows the resonant frequencies for cavity modes888, , and are associated with eigenmode orders along azimuthal (), radial (), and longitudinal directions (), respectively. and the normalized error defined as where and are numerical and analytic resonant frequencies, respectively. It is seen that all resonant cavity modes are axisymmetric () and the normalized errors are fairly low.

Figure 9: Spectrum for resonant cavity modes from MHz to GHz.
[GHz] [%] [GHz] [%]
0.941 0.016
Table 1: Resonant frequencies for axisymmetric cavity modes and normalized errors between numerical and analytic works.

4.2 Space-charge-limited (SCL) cylindrical diode

As a second validation example, we test the accuracy of present EM-PIC algorithm by modeling a SCL cylindrical diode with finite-length emitter. By applying an external voltage to the cathode, a fast rise in the number of electrons emitted from the cathode initially occurs; however, in the steady-state the injection current density eventually becomes saturated due to space charge effects. For an infinitely long cylindrical diode and electrodes, Langmuir-Blodgett’s law describes the SCL current per unit length as


where is the external voltage, and are charge and mass of electrons, is a radial coordinate, and


with and denotes the radius of the cathode. Here, we consider the emitting electrode with finite length . The solution domain is with mm, mm, and mm. The domain has the horizontal wall segments representing electrode surfaces (cathode or anode), as shown in Fig. 10. If the inner conductor is chosen as a cathode, the electron flow becomes divergent, on the other hand, it is convergent. The left and right boundaries of the domain are truncated by a perfectly matched layer (PML) teixeira1998general ; teixeira1999unified ; WangTPS2006 . The unstructured mesh has nodes, edges, and faces. We choose ps and the simulation runs up to a total time-steps.

Figure 10: Schematics for divergent and convergent flows in the cylindrical diode.

Fig. 10 illustrates snapshots for electron beam distribution with mm and kV. In order to determine , we first fix the superparticle scaling factor , indicating the number of electrons for each superparticle, and gradually increase the superparticle injection rate until the virtual cathode starts to form. We simulate a total of cases, including , , , and for each divergent and convergent electron flow, and compare the maximum injection current density without formation of a virtual cathode, to the previous results obtained by kostov2002space with KARAT, which is a FDTD-based EM-PIC algorithm. Fig. 11 shows versus for divergent and convergent electron flows. Red-solid (divergent) and -dashed (convergent) lines are KARAT results and blue-markers with upper ranges are obtained from present EM-PIC simulations. The upper ranges on the blue markers stand for an interval where an exact solutions for the current density may exist. The marker indicates maximum current density without the virtual cathode formation and upper horizontal line is minimum injection current density with virtual cathode formation. Two markers are used because of the stepwise increases in the current density by the assumed superparticle number in our EM-PIC model. Decreasing superparticle scaling factor or increasing superparticle injection rate can yield higher resolution. Fig. 11 shows very good agreement between the results of the present EM-PIC algorithm and those of KARAT, for both divergent and convergent flows. In addition, it is seen that, as increases, the current density converges to the limit of Langmuir-Blodgett’s law as expected.

Figure 11: Space-charge-limited current density for various and comparison between present EM-PIC simulations and KARAT by kostov2002space .

Fig. 12 shows the magnitude of the electrical self-field and external field at the instant of virtual cathode formation. Since the external field is stronger as decreases, a divergent flow produces more charges on the cathode surface so that their self-field may cancel the external field. As a result, the current density of divergent flows is larger than that of convergent flow.

Figure 12: Electric field intensity of self- and external fields at the instant of virtual cathode formation.

5 Numerical examples

In this section, we present simulations of a relativistic BWO device operating at -point by using the proposed EM-PIC algorithm. First we consider a SWS with sinusoidal ripples on a cylindrical waveguide section and determine its characteristics by performing a “cold” test (i.e. without the presence of an electron beam). Then, we perform “hot” tests (with electron beams) of BWO to check the reliability and validity of our axisymmetric EM-PIC algorithm. In particular, we are interested in investigating effects of staircasing errors on the predicted behavior of the BWO system.

Figure 13: Schematics of backward-wave oscillator with an instant particle distribution snapshots at  ns.
Figure 14: Electric potential distribution (contour plots) and corresponding electric fields (vector plots) between the cathode and the anode.

5.1 Relativistic backward-wave oscillator (BWO)

Consider a BWO system composed of the cathode-anode, the SCSWS region, the beam collector, the output port, and the beam focusing system. In order to produce the relativistic electron beam, we apply an external voltage difference between cathode and anode. We choose kV which produces the electron beam with mean axial velocity with the width of mm. Fig. 14 illustrates electric potential distribution (contour plots) and corresponding electric fields (vector plots) by solving Poisson equations. Each super-particle represents electrons so that the total injection current is around kA. On average, macro-particles are assigned to each grid cell. The Debye length in our simulations is  mm and much larger than the average grid (edge) length,  mm. This avoids artificial heating of the electron beam. The electron beam is emitted from the cathode and eventually absorbed at the collector. We consider a SCSWS with radial profile , where and are maximum and minimum radii, respectively, and is the corrugation period. The total number of corrugations along the structure is denoted as . Based on an eigenmode analysis 999Corresponding to the absence of an electron beam or a so-called ‘cold test’. (see the double-refined case in Fig. 24), the SCSWS was designed to have  mm,  mm,  mm, and for Ku band operation. We terminate the output ports of the BWO system by inserting a PML. All left vertical walls except for the cathode are truncated by PML to avoid spurious reflections. By using a PML with thickness equals to where is wavelength of the center frequency, the PML reflection coefficient, defined in donderici2008mixed , is as low as dB. In the beam focusing system, a static axial magnetic field is applied over the SWS region. We enforce axisymmetric fields at the -axis, by applying a perfect magnetic conductor (PMC) boundary condition there (note that only the component is present on axis). For the particles incident on the axis we use a perfectly reflecting boundary condition.

Figure 15: A zoomed-in region of four rightmost corrugations of Fig. 13 with RGB color scales reflecting particle velocities.

5.1.1 System performance

Fig. 13 illustrates a snapshot of the electron beam at  ns. In Fig. 15 we plot a zoomed-in beam picture that shows the four rightmost corrugations, where RGB colors indicate variations in normalized particle velocities . The periodic particle beam bunching is a result of particles being accelerated or decelerated, which means that the beam electrons synchronously lose and recover their kinetic energy. The net energy is transferred from the beam to the waves as seen in Fig. 16, which illustrates a phase space distribution of the particle beam at  ns. The particles decelerated from the initial velocity ( m/s) dominate the accelerated particles, so that a net loss of the beam kinetic energy leads to amplification of the mode.

Figure 16: Phase-space plot at ns.

Fig. 17 presents a vector plot of self-fields generated by the electron beam at  ns (in steady state). Coherent Cerenkov beam-wave interactions give rise to a strong mode that may be observed within the SWS region.

Figure 17: A snapshot of steady-state self-fields ( ns).

The self-fields at the output port were analyzed in time and frequency domains, as shown in Fig. (a)a and Fig. (b)b, respectively. It can be seen in Fig. (a)a that although initially there are no oscillations, the output signal starts to oscillate simultaneously with the onset of beam bunching. This signal keeps evolving and eventually approaches a steady state at about  ns. By performing the Fourier analysis of steady-state output signals, we show that their spectrum shows a good degree of single-mode purity at  GHz, as shown in Fig. (b)b.

Figure 18: Output signal analysis in (a) time and (b) frequency domains.
Figure 19: Verification of charge conservation at nodes along time (at time-steps of , , ) by testing NR levels of (a) DCE and (b) DGL.

In order to verify charge conservation, we plot the normalized residuals (NR) versus the nodal indices for the DCE and DGL na2016local in Figs. (a)a and Fig. (b)b, respectively. These residuals, evaluated for each time step or and node , are defined as


Fig. (a)a shows that NR levels for DCE remain fairly low at all nodes and very close to the double precision floor (below ). Sparsely observed peaks are due to the fact that the electron-beam edges, where electrons occasionally travel through extremely small cell fractions, generate numerical noise during the scatter step. However, these errors still remain well below . NR levels for DGL are also distributed around the double precision floor as shown in Fig. (b)b, which means that charge conservation is maintained within round-off errors.

Figure 20: 3D velocity plots for an electron beam with the BFS magnetic field of  T.

The transverse dynamics of an electron beam in the cross-section plane are sensitive to gyroradius value which depends on the BFS magnetic field strength. Fig. 20 shows 3D electron velocity plots for a weaker BFS magnetic field reduced from  T to  T. RGB colors again indicate particle velocity magnitudes. The weaker axial magnetic field yields a relatively larger in the polar () plane, which makes the electron beam to gradually expand radially.

5.1.2 Staircasing error analysis

We now examine staircasing errors resulting from a discrete approximation of the SCSWS boundary. We consider eight cases comprising the same BWO geometry modeled by: () a coarse mesh, () a double-refined mesh, () a quadruple-refined mesh, () a octuple-refined mesh, () a coarse mesh with staircased boundary, () a double-refined mesh with staircased boundary, () a quadruple-refined mesh with staircased boundary, and () a octuple-refined mesh with staircased boundary. Fig. 21 depicts one period of the SCSWS boundary as rendered by the different meshes.

Figure 21: SCSWS boundary profiles for all cases.

All cases employ unstructured meshes. Although unstructured meshes do not necessitate staircased boundaries, the latter are enforced to reproduce the staircasing that would be present on structured meshes. At the same time, this setup allows for a detailed study that effectively isolates staircasing error from grid-dispersion errors. The mesh information for all cases is shown in Table 2.

Case ave [mm]
coarse 4,564 12,978 8,415 2.50
double-refined 5,801 16,545 10,745 1.67
quadruple-refined 8,771 25,236 16,466 1.11
octuple-refined 15,333 44,542 29,210 0.74
coarse staircased 4,593 13,005 8,413 2.50
double-refined staircased 5,909 16,766 10,858 1.67
quadruple-refined staircased 8,939 25,547 16,610 1.11
octuple-refined staircased 15,543 44,904 29,362 0.74
Table 2: Mesh information for different SCSWS cases

We have performed hot tests for all cases, keeping the number of superparticle per each cell by on average.

Figure 22: Field signal at the output port in (a) SCSWS and (b) staircased SCSWS in the time domain.
Figure 23: Normalized spectral amplitude at the output port in SCSWS and staircased SCSWS.

Field output signals in time and oscillation frequencies are displayed in Fig. 22 and Fig. 23. The results based on meshes devoid of staircasing error presented in Fig. (a)a converge much faster than those of meshes with staircased boundaries in Fig. (b)b The oscillation frequency values in Fig. 23 also shows the fast convergent rate of present EM-PIC simulations. In particular, it is interesting to see that the results from the BWO with a staircased SCSWS are unable to capture much RF oscillation at all in the case of coarse and double-refined meshes (case (v) and (vi)). The RF oscillation is more visible on quadruple- or octuple-refined meshes (case (vii) and (viii)). This is because underestimation of the -point frequency101010-point denotes the solution of the modal field which have maximum frequency in the passband. in the modal dispersion causes a slow-charge mode driven by the electron beam benford2015high to falls into the forward wave region and, as a result, the system does not act as an oscillator anymore. Note that the BWO system here is designed to operate at -point of the modal fields. The underestimation of the -point frequency in the staircased boundary is shown by Fig. 24, which depicts the dispersion diagrams for the mode of each SWS. The analytic dispersion relation for the mode is obtained based on Floquet’s theory by considering harmonics up to order barroso2002cylindrical . It is clearly seen that the -point frequency decreases as the staircasing error become significant. On the other hand, the cases without staircasing errors rapidly converge to the analytic prediction.

Figure 24: Dispersion relations from “cold tests”.

6 Concluding Remarks

We introduced a new axisymmetric charge-conservative EM-PIC algorithm on unstructured grids for the analysis and design of micromachined VEDs with cylindrical axisymmetry. We demonstrated that cylindrical symmetry in device geometry, fields, and sources enables reduction of the 3D problem to a 2D one through considering -polarized fields in the meridian () plane and introducing an artificial inhomogeneous medium to account for the metric elements of the cylindrical coordinate system. As a result, computational resources are significantly reduced. The unstructured-grid spatial discretization is achieved by using Whitney forms in the meridian plane. Using leapfrog time integration, we obtain space- and time-discretized Maxwell’s equations in the form of a mixed FETD scheme. A local and explicit time update is made possible by employing the SPAI approach na2016local . Interpolation of the field values to the particles’ positions is also performed by Whitney forms and used for solving the Newton-Lorentz equations of motion of each particle. Relativistic particles are accelerated based on a relativistic Boris algorithm with correction. In the particle scatter step, we employ a recently introduced Cartesian charge-conserving scatter scheme for unstructured grids moon2015exact . The algorithm was validated by considering cylindrical cavity and space-charge-limited cylindrical diode problems. We also illustrated the advantages of the present algorithm in the analysis of a BWO system including a slow-wave waveguide structure with complex geometry.


This work was supported in part by U.S. NSF grant ECCS-1305838, OSC grants PAS-0061 and PAS-0110, FAPESP-OSU grant 2015/50268-5. H. Moon and D.-Y. Na were also supported by the OSU Presidential Fellowship program.

Appendix A Pairing operation and generalized Stokes theorem

One of the key properties of Whitney -forms is that they admit a natural “pairing” with the -cells of the mesh, where refers to the dimensionality i.e. refers to nodes, to edges, and so on teixeira1999lattice . Computationally, the pairing operation between an -th -cell of the grid and the Whitney form associated with the -th -cell is effected by the integral below and yields teixeira1999lattice ; teixeira2014lattice


for in the 3D space and where is the Kronecker delta. The generalized Stokes theorem teixeira1999lattice ; teixeira2014lattice ; cairns1936the states


where is the boundary operator. Using the above pairing operation and generalized Stokes’ theorem, we can obtain discrete Maxwell’s equation on a irregular lattice (unstructured grid). For example, applying pairing for -cells into Faraday’s law on the primal mesh gives


and applying generalized Stokes’ theorem into the left-hand side term of (49) yields


Substituting (5) and (6) into (50) and using where is an incidence matrix element taking a value in the set of  teixeira1999lattice ; teixeira2014lattice ; hughes1981lagrangian ; guth1980existence , we obtain


By using (47), (51) can be rewritten as


for . (52) represents the discrete representation of Faraday’s law as written in (19). Discrete Ampere’s law can be obtained by a similar procedure on the dual mesh.

Appendix B Discrete Hodge matrix

Using vector calculus proxies he2007differential ; kim2011parallel ; tarhasaari1999some ; gillette2011dual ; he2006geometric and considering a (doubly) artificial inhomogeneous medium at the meridian plane as done in (12) and (14), we can rewrite (21) and (22) as


where is a vector proxy for Whitney -form associated with the -th -cell. The vector proxies for Whitney 1- and 2-forms are written as moon2015exact