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

TL;DR

This paper introduces a charge-conservative electromagnetic particle-in-cell algorithm optimized for axisymmetric vacuum electronic devices, utilizing unstructured grids and transformation optics principles to improve simulation accuracy and efficiency.

Contribution

The paper presents a novel 2.5D EM-PIC algorithm that combines unstructured mesh discretization, transformation optics, and explicit time-stepping for efficient simulation of axisymmetric VEDs.

Findings

Validated with simulations of VEDs with particle bunching effects

Achieved accurate charge conservation on unstructured grids

Demonstrated effective modeling of resonance Cerenkov interactions

Abstract

We present a 2.5-dimensional charge-conservative electromagnetic particle-in-cell (EM-PIC) algorithm optimized for the analysis of vacuum electronic devices (VED) with cylindrical symmetry (axisymmetry). We explore the axisymmetry present in the device geometry, fields, and sources to reduce the dimensionality of the problem from 3D to 2D. Further, we explore `transformation optics' principles to map the original problem in polar coordinates to an equivalent problem on Cartesian coordinates with an effective (artificial) inhomogeneous medium introduced. The resulting problem in the meridian plane is discretized using an unstructured 2D mesh considering TE-polarized fields and properly scaled charges. EM field and source variables (node-based charges and edge-based currents) 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 a relativistic correction. In the scatter step, we apply a radial scaling factor on top of a charge-conserving scatter scheme tailored for 2-dimensional unstructured grids. As validation examples, we demonstrate simulations that investigate the physical performance of VEDs designed to harness particle bunching effects arising from the coherent (resonance) Cerenkov electron beam interactions within micromachined slow-wave structures.

Tables

Table 1. Table 1: Resonant frequencies for axisymmetric cavity modes and normalized errors between numerical and analytic works.

$(m,n,p)$	$f_a$ [GHz]	$\frac{\left|f_a-f_s\right|}{f_a}\times 100$ [%]	$(m,n,p)$	$f_a$ [GHz]	$\frac{\left|f_a-f_s\right|}{f_a}\times 100$ [%]
$(0,1,0)$	$0.230$	$0.083$	$(0,1,5)$	$0.784$	$0.0035$
$(0,1,1)$	$0.274$	$0.040$	$(0,2,4)$	$0.799$	$0.0048$
$(0,1,2)$	$0.378$	$0.040$	$(0,3,0)$	$0.826$	$0.0014$
$(0,1,3)$	$0.505$	$0.047$	$(0,3,1)$	$0.840$	$0.019$
$(0,2,0)$	$0.527$	$0.042$	$(0,3,1)$	$0.879$	$0.017$
$(0,2,1)$	$0.548$	$0.031$	$(0,3,1)$	$0.917$	$0.027$
$(0,2,2)$	$0.607$	$0.0030$	$(0,2,5)$	$0.929$	$0.027$
$(0,1,4)$	$0.642$	$0.023$	$(0,3,3)$	0.941	0.016
$(0,2,3)$	$0.693$	$0.034$	$-$	$-$	$-$

Table 2. Table 2: Mesh information for different SCSWS cases

Case	$N_0$	$N_1$	$N_2$	ave$\left(l_{\text{edge}}\right)$ [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

Equations

$$d\mathcal{E}$$

$$d\mathcal{H}$$

$$d\mathcal{D}$$

$$d\mathcal{B}$$

$$\mathcal{E}$$

$$\mathcal{B}$$

$$\mathcal{D}$$

$$\mathcal{H}$$

$$\mathcal{D}$$

$$\mathcal{H}=H_{\varphi }d\varphi ,$$

$${\mu }_0^{-1}\star \mathcal{B}$$

$$\mathcal{D}=D_{\rho }d\varphi \wedge dz$$

$${ϵ}_0\star \mathcal{E}$$

$$\mathcal{E}(t)=\sum _{i=1}^{N_1}{\mathbb{E}}_i(t)w_i^1,$$

$$\mathcal{B}(t)=\sum _{i=1}^{N_2}{\mathbb{B}}_i(t)w_i^2,$$

$${\mathcal{J}}_{\star }(t)=\sum _{i=1}^{N_1}{\mathbb{J}}_{\star ,i}(t)w_i^1,$$

$${\mathcal{Q}}_{\star }(t)=\sum _{i=1}^{N_0}{\mathbb{Q}}_{\star ,i}(t)w_i^0$$

$$\left[{\mathcal{D}}_{\text{curl}}\right]\cdot \left[\mathbb{E}\right]$$

$$\left[{\tilde{\mathcal{D}}}_{\text{curl}}\right]\cdot \left[{\star }_{{\mu }^{-1}}\right]\cdot \left[\mathbb{B}\right]$$

$${\left[{\star }_{ϵ}\right]}_{i,j}$$

$${\left[{\star }_{{\mu }^{-1}}\right]}_{i,j}$$

$${\left[\mathbb{B}\right]}^{n+\frac{1}{2}}$$

$${\left[\mathbb{E}\right]}^{n+1}$$

$${\vec{E}}^{\text{EM}}\left({\vec{r}}_{\alpha },n\Delta t\right)$$

$${\vec{B}}^{\text{EM}}\left({\vec{r}}_{\alpha },n\Delta t\right)={\vec{B}}_{\alpha }^{\text{EM},n}$$

$$\frac{d{\vec{r}}_{\alpha }}{dt}$$

$$\frac{d\left(m_{\alpha }{\gamma }_{\alpha }{\vec{v}}_{\alpha }\right)}{dt}$$

$$\frac{{\vec{r}}_{\alpha }^{n+1}-{\vec{r}}_{\alpha }^n}{{\Delta }_t}$$

$$\frac{{\vec{u}}_{\alpha }^{n+\frac{1}{2}}-{\vec{u}}_{\alpha }^{n-\frac{1}{2}}}{\Delta t}$$

$${\mathbb{Q}}_{\star ,i}^n$$

$${\mathbb{J}}_{\star ,i}^{n+\frac{1}{2}}$$

$$\frac{{\left[\mathbb{Q}\right]}^{n+1}-{\left[\mathbb{Q}\right]}^n}{\Delta t}+\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[\mathbb{J}\right]}^{n+\frac{1}{2}}=0$$

$$\frac{{\mathbb{Q}}_{\star ,1,\alpha }^{n+1}-{\mathbb{Q}}_{\star ,1,\alpha }^n}{\Delta t}=\frac{q_{\alpha }}{\Delta t}\left({\lambda }_{1,\alpha }^{n+1}-{\lambda }_{1,\alpha }^n\right).$$

$${\mathbb{J}}_{\star ,1,\alpha }^{n+\frac{1}{2}}+{\mathbb{J}}_{\star ,2,\alpha }^{n+\frac{1}{2}}$$

$$\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[{\star }_{ϵ}\right]}_a\cdot {\left[\mathbb{E}\right]}^{n+1}=\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[{\star }_{ϵ}\right]}_a\cdot {\left[\mathbb{E}\right]}^n-\Delta t\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[\mathbb{J}\right]}^{n+\frac{1}{2}}.$$

$$\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[{\star }_{ϵ}\right]}_a\cdot \left(\frac{{\left[\mathbb{E}\right]}^{n+1}-{\left[\mathbb{E}\right]}^n}{\Delta t}\right)=-\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[\mathbb{J}\right]}^{n+\frac{1}{2}},$$

$$\left[{\tilde{\mathcal{D}}}_{\text{div}}\right]\cdot {\left[{\star }_{ϵ}\right]}_a\cdot \left(\frac{{\left[\mathbb{E}\right]}^{n+1}-{\left[\mathbb{E}\right]}^n}{\Delta t}\right)=\frac{{\left[\mathbb{Q}\right]}^{n+1}-{\left[\mathbb{Q}\right]}^n}{\Delta t}.$$

Full text

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) [1, 2, 3, 4, 5], 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 [6, 7].

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) [2, 8, 9, 10, 11]. SWSs are often made by imposing periodic ripples on the conducting wall of cylindrically symmetric waveguides, as illustrated in Fig. 1a, so that the phase velocity of the modal field becomes slower than the speed of light in vacuum due to Bragg scattering [12]. 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 [2]. 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 [13, 14], and tapering ripples may reduce reflections at the output of TWTA and prevent internal oscillations [8, 15]. 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 [16, 17]. 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 [18]. 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.

Computational experiments for VEDs employ electromagnetic particle-in-cell (EM-PIC) algorithms [19, 20, 21, 22], 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 [19, 20, 23]. 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 [24, 25, 26]. However, the relatively poor grid-dispersion properties of this algorithm [27] causes spurious numerical Cerenkov radiation [28]. 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 [29, 30].

On the other hand, the finite-element time-domain (FETD) method [31, 32] 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 [33, 34]. In addition, implicit time updates used in conventional FETD require repeated linear solves at each time-step [31]. Recently, a novel charge-conserving scatter scheme for unstructured grids, inspired by differential-geometric ideas and the exterior calculus of differential forms [35, 36, 37, 38, 39, 40], has been proposed in [41]. Other charge-conservative EM-PIC algorithms for unstructured grids were also developed under similar tenets in [42, 43, 44]. 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 [45, 46, 47] 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 $\phi$), 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 [30]. 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 [37, 48, 49, 50]., considering $\text{TE}^{\phi}$-polarized fields in the meridian ($\rho z$) 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 $\mathcal{E}-\mathcal{B}$ 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 [47]. 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 [51, 52]. In the scatter step, we adapt an exact charge-conserving scheme developed for Cartesian unstructured grids [41]. 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. 1b). 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 [36, 38, 40, 42, 53] as

[TABLE]

where $\mathcal{E}$ and $\mathcal{H}$ are 1-forms for the electric and magnetic field intensity, $\mathcal{D}$ and $\mathcal{B}$ are 2-forms for the electric and magnetic flux density, $\mathcal{J}$ is 2-form for the electric current density, $\mathcal{Q}$ is 3-form for the electric charge density, and operator $d$ is the exterior derivative encompassing conventional curl and divergence operators [54, 55, 56, 57]. These 1-, 2-, and 3-forms can be expressed using a set of non-orthonormal-basis in a cylindrical coordinate system $\left(d\rho,d\phi,dz\right)$ [48, 54]. For instance, $\mathcal{E}$ (1-form) is expressed as $\mathcal{E}=E_{\rho}d\rho+E_{\phi}d\phi+E_{z}dz$; then, its vector proxy can be written as $\vec{E}=E_{\rho}\hat{\rho}+\frac{E_{\phi}}{\rho}\hat{\phi}+E_{z}\hat{z}$. Similarly, $\mathcal{B}$ (2-form) is given by $\mathcal{B}=B_{\rho}d\phi\wedge dz+B_{\phi}dz\wedge d\rho+B_{z}d\rho\wedge d\phi$ (where $\wedge$ is the exterior or wedge product [35, 36, 55]) and its vector proxy takes the form of $\vec{B}=\frac{B_{\rho}}{\rho}\hat{\rho}+B_{\phi}\hat{\phi}+\frac{B_{z}}{\rho}\hat{z}$ [54].

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 $\partial\Omega=R\left(z\right)$ where $\Omega$ is the computational domain and $R$ is the wall radius which only depends on $z$, as shown in Fig. 2a.

Cylindrical axisymmetry, used here, implies that there is no variation along $\phi$ (${\partial}/{\partial\phi}=0$) 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. 2b). There exist two useful constraints that simplify the original 3D problem: (i) elimination of the $d\phi\frac{\partial}{\partial\phi}$ term in the exterior derivative $d$, viz. $d=d\rho\frac{\partial}{\partial\rho}+dz\frac{\partial}{\partial z}$ and (ii) retainment of only transverse magnetic (TM) eigenmodes with $m=0$444The index $m$ 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

[TABLE]

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 $\text{TE}^{\phi}$-polarized fields on the meridian plane (see also Fig. 3).

2.3 Modified Hodge star operator

The Hodge star operator $\star$ maps $p$-forms into $\left(n-p\right)$-forms in $n$-dimensional space555The Hodge star operator can be understood geometrically as yielding the orthogonal complement of a given differential form to the volume form in $n$-space. In the 3D case for example, $\star(dx)$ = $dy\wedge dz$, $\star(dz\wedge d\rho)$ = $\rho d\phi$, $\star(1)$ = $dx\wedge dy\wedge dz$ = $\rho\,d\rho\wedge d\phi\wedge dz$ and so forth [54]. For an arbitrary $p$-form ${\cal A}$ in $n$-space, we have ${\cal A}\wedge(\star{\cal A})=|A|^{2}dV$, where $dV$ is the volume element in $n$-space. As such, ${\cal A}\wedge(\star{\cal A})$ provides the $L^{2}$ ( or“energy”) element norm of ${\cal A}$. [39, 40, 54, 57, 58, 59]. In our case, we have

[TABLE]

The Hodge operators incorporate the metrical properties of the system, which in the cylindrical case are expressed by a metric tensor $\mathrm{diag}(1,\rho^{2},1)$. For the magnetic field and flux density, substituting (8) into the left-hand side term of (9) gives

[TABLE]

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

[TABLE]

Note that $\star$ only acts on the differentials such as $dz$, $d\rho$, and $d\phi$. By comparing (11) and (12) and introducing an artificial magnetic permeability, $\mu\left(\rho\right)=\mu_{0}\rho^{-1}$, we can extract the radial factor $\rho$ from the Hodge star operator. As a result, we can reuse simple Cartesian space calculus with metric tensor $\mathrm{diag}(1,1,1)$ 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

[TABLE]

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

[TABLE]

In a similar fashion, an artificial electrical permittivity takes the form of $\epsilon\left(\rho\right)=\epsilon_{0}\rho$. Essentially, an original 3D problem with cylindrical axisymmetry is replaced by an equivalent 2D problem with $\text{TE}^{\phi}$-polarized fields immersed in Cartesian space with an inhomogeneous medium with artificial permittivity and permeability, as illustrated schematically in Fig. 3.

As a result, the present field-solver borrowing the concept of TO is free from the axial singularity associated with the $1/\rho$ 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 $\mathcal{E}-\mathcal{B}$ 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 $\mathcal{E}$, $\mathcal{B}$, $\mathcal{D}$, and $\mathcal{H}$ as a 1-, 2-, 1-, and 0-forms, respectively, and $\mathcal{J}_{\star}$ and $\mathcal{Q}_{\star}$ as 1- and 0-forms 666Since $\mathcal{J}$ and $\mathcal{Q}$ are 1- and 2-forms, respectively.. It should be stressed that the operator $\star$ used here for the equivalent problem is defined in the Cartesian 2D space and hence devoid of the radial factor $\rho$ (which enters instead in the definition of the artificial medium properties). Because $\mathcal{E}$ and $\mathcal{B}$ are ordinary differential forms with internal orientations associated with the primal lattice [40, 57], we expand them by using Whitney 1- and 2-forms (cochains) residing in the primal mesh, as follows:

[TABLE]

where $N_{p}$ is the number of $p$-cells for $p=0,1,2$ (number of nodes, edges, facets, respectively, on the grid), $w_{i}^{p}$ is the Whitney $p$-form associated with the $i$-th $p$-cell on the grid, and $\mathbb{E}_{i}(t)$ and $\mathbb{B}_{i}(t)$ are the (time-evolving) degrees of freedoms (DoFs) for $\mathcal{E}$ and $\mathcal{B}$, respectively [60]. On the other hand, $\mathcal{D}$, $\mathcal{H}$, $\mathcal{J}$, and $\mathcal{Q}$ are twisted differential forms with external orientations, associated with the dual lattice [37, 40, 55]. Fig. 4 shows schematically how the DoFs are associated to the primal and dual meshes elements in our 2D setting.

To avoid the need for working with the dual mesh, we employ the Hodge duals ${\mathcal{J}}_{\star}$ and $\mathcal{Q}_{\star}$ which can be expanded using Whitney 1- and 0-forms in the primal mesh, as follows:

[TABLE]

where $\mathbb{J}_{\star,i}$ and $\mathbb{Q}_{\star,i}$ are Dos of $i$-th edge and node for ${\mathcal{J}}_{\star}$ and $\mathcal{Q}_{\star}$, 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 [37, 40, 61], the spatially discretized Maxwell’s equations (see also A) are given by

[TABLE]

with $\left[\mathbb{E}\right]=\left\{\mathbb{E}_{1},...,\mathbb{E}_{N_{1}}\right\}^{T}$, $\left[\mathbb{B}\right]=\left\{\mathbb{B}_{1},...,\mathbb{B}_{N_{2}}\right\}^{T}$, $\left[\mathbb{J}_{\star}\right]=\left\{\mathbb{J}_{\star,1},...,\mathbb{J}_{\star,N_{1}}\right\}^{T}$, and where $\left[\mathcal{D}_{\text{curl}}\right]$ and $\left[\mathcal{\tilde{D}}_{\text{curl}}\right]$ 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 [37, 40, 62, 63]. Rows of $\left[\mathcal{D}_{\text{curl}}\right]$ have a $1:1$ association to faces of the primal mesh and columns to edges. On a primal triangular mesh, each row of $\left[\mathcal{D}_{\text{curl}}\right]$ has three nonzero elements corresponding to the three edges of the triangular cell associated to that row (with entries $1$ or $-1$ 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 $\left[\mathcal{\tilde{D}}_{\text{curl}}\right]$ so that both of these matrices take their elements in the set of $\left\{-1,0,1\right\}$ as a consequence of their metric-free character. It can be shown that $\left[\mathcal{\tilde{D}}_{\text{curl}}\right]$=$\left[\mathcal{D}_{\text{curl}}\right]^{T}$ [37]. In (19) and (20), $\left[\star_{\epsilon}\right]$ and $\left[\star_{\mu^{-1}}\right]$ are discrete Hodge operators (see B), which are symmetric positive-definite matrices obtained by the Galerkin method [46, 58, 59, 64]. Their elements are expressed as the area (volume in 3D) integrals

[TABLE]

where $K$ 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 $n$ and the magnetic flux density, field and current density defined at half integer time-steps $n\pm\frac{1}{2}$. The resulting space- and time-discretized Maxwell’s equations form the so-called mixed $\mathcal{E}-\mathcal{B}$ finite-element time-domain (FETD) scheme [41, 46, 47, 65]:

[TABLE]

where superscripts $n$ denote time-step and $\Delta t$ is the time-step increment. Local explicit time-update is achieved using the SPAI approach to obtain a representation for $\left[\star_{\epsilon}\right]^{-1}$ and obviate the need for a (expensive) linear-solver at each time-step [46, 47]. 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 $\alpha$-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 [41]

[TABLE]

where $\vec{W}_{i}^{p}$ is a vector proxy of the Whitney $p$-form for $i$-th grid element and $\vec{r}_{\alpha}$ 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 ($\rho\phi$) 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 $\alpha$-species particles are updated by solving the Newton-Lorentz equations of motion as [47].

[TABLE]

where $\vec{v}$ is a velocity vector, $m_{\alpha}$ is a particle mass, $q_{\alpha}$ is a particle charge, and $\gamma=\frac{1}{\sqrt{1-\left(\frac{\left|\vec{v}\right|}{c}\right)^{2}}}$, and $\vec{u}=\gamma\vec{v}$. The leapfrog time-discretized Newton-Lorentz equations of motion can be written as

[TABLE]

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 [51, 52].

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 [43, 47]. A charge-conserving scatter algorithm has been recently proposed [41] 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 [41]:

[TABLE]

where $i$ represents the $i$-th grid element, and $i_{a}$ and $i_{b}$ denote indices of lower dimensional manifolds of $i$-th element777For example, if $i$-th element is an edge, $i_{a}$ and $i_{b}$ indicate the two nodes associated with the said edge.. The functions $W_{i}^{0}$ and $\vec{W}_{i}^{1}$ are scalar and vector proxies of Whitney 0- and 1-forms, respectively, and $\lambda_{i,\alpha}^{n}$ is a barycentric coordinate for an $\alpha$-species particle at time-step $n$, respectively. One can find the derivations of Eqs. (31) and (32) in more detail in [41].

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, $n-2$ to time-step, $n+1$. As shown in Fig. 5, the grid charge (blue-colored circle) at node $1$ is calculated based on (31). The associated charge value is represented by the colored triangular area formed by the particle position and nodes $2$ and $3$. Likewise, as shown in Fig. 6, the respective grid current at edge $1$ (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 [41, 47]

[TABLE]

where $\left[\tilde{\mathcal{D}}_{\text{div}}\right]$ 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 $\left\{-1,0,1\right\}$, see also the Appendices. Note that Hodge matrices for the discrete Hodge duals of charges and currents ($\left[\mathbb{Q}_{\star}\right]$ and $\left[\mathbb{J}_{\star}\right]$) become identity matrices for barycentric dual lattices. Therefore, $\left[\mathbb{Q}\right]=\left[\mathbb{Q}_{\star}\right]$ and $\left[\mathbb{J}\right]=\left[\mathbb{J}_{\star}\right]$.

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 $n$ to time-step $n+1$ in a given triangular element (see Fig. 7). Then, the rate of change of the nodal charge at node $1$ is

[TABLE]

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

[TABLE]

Thus, the rate of change of the nodal charge at node $1$ 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 $\left[\tilde{\mathcal{D}}_{\text{div}}\right]\cdot\left[\star_{\epsilon}\right]_{\text{a}}$, where $\left[\tilde{\mathcal{D}}_{\text{div}}\right]$ is the incidence matrix representing the discrete divergence operator on the dual grid and $\left[\star_{\epsilon}\right]_{\text{a}}$ is an exact inverse of $\left[\star_{\epsilon}\right]]_{\text{a}}^{-1}$ obtained by using SPAI algorithm, and utilizing the identity $\left[\tilde{\mathcal{D}}_{\text{div}}\right]\cdot\left[\tilde{\mathcal{D}}_{\text{curl}}\right]=0$ [37, 40] yields

[TABLE]

which can be rearranged as

[TABLE]

and, using (33), rewritten as

[TABLE]

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

[TABLE]

and by induction,

[TABLE]

for all $n$. This means that if consistent initial conditions $\left[\tilde{\mathcal{D}}_{\text{div}}\right]\cdot\left[\star_{\epsilon}\right]_{a}\cdot\left[\mathbb{E}\right]^{0}=\left[\mathbb{Q}\right]^{0}$ 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 $\text{TM}_{0np}$ 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 $a=0.5$ m and height $h=1$ m (see Fig. 8). Two magnetic point sources $M_{1}\left(\rho,z,t\right)$ and $M_{2}\left(\rho,z,t\right)$ excited by Gaussian-modulated pulses of broadband spectrum are placed at arbitrary positions given as

[TABLE]

with $t_{g}=40$ ns, $\sigma_{g}=1$ ns, and $\omega_{g}=\pi\times 10^{9}$ 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 $8,095$ nodes, $23,928$ edges , and $15,834$ faces. The three lateral metallic boundaries are assumed as perfect electric conductors (PEC). The remaining boundary is the $z$-axis (axisymmetric boundary). The time step interval $\Delta{t}$ is chosen as $1$ ps, and the simulation runs over a total of $2\times 10^{6}$ time steps. Fig. 8 shows a snapshot for electric field distribution in the cavity at $2$ $\mu$s. The RGB colormap and the white arrows indicate magnitudes and vectors of the electric fields, respectively.

The time signal are detected at $\left(\rho=0.71,z=0.39\right)$ (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 $1$ MHz to $1$ 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 $\text{TM}_{mnp}$ cavity modes888$m$, $n$, and $p$ are associated with eigenmode orders along azimuthal ($\phi$), radial ($\rho$), and longitudinal directions ($z$), respectively. and the normalized error defined as $\frac{\left|f_{s}-f_{a}\right|}{f_{a}}$ where $f_{s}$ and $f_{a}$ are numerical and analytic resonant frequencies, respectively. It is seen that all resonant cavity modes are axisymmetric ($m=0$) and the normalized errors are fairly low.

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 $J_{1D,LB}$ as

[TABLE]

where $V$ is the external voltage, $e$ and $m$ are charge and mass of electrons, $\rho$ is a radial coordinate, and

[TABLE]

with $\mu=\ln\left(\rho/\rho_{c}\right)$ and $\rho_{c}$ denotes the radius of the cathode. Here, we consider the emitting electrode with finite length $L_{e}$. The solution domain is $\Omega=\left\{\left(\rho,z\right)\in\left[\rho_{i},\rho_{o}\right]\times\left[0,L_{z}\right]\right\}$ with $\rho_{i}=5$ mm, $\rho_{o}=15$ mm, and $L_{z}=100$ 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) [66, 67, 68]. The unstructured mesh has $7,313$ nodes, $21,334$ edges, and $14,022$ faces. We choose $\Delta t=0.15$ ps and the simulation runs up to a total $30,000$ time-steps.

Fig. 10 illustrates snapshots for electron beam distribution with $L_{e}=10$ mm and $V=1$ kV. In order to determine $J_{2D,\text{max injection}}$, we first fix the superparticle scaling factor $C_{sp}$, 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 $8$ cases, including $L_{e}/\rho_{o}=0.4$, $1$, $2.425$, and $4$ for each divergent and convergent electron flow, and compare the maximum injection current density without formation of a virtual cathode, $J_{2D,\text{max injection}}$ to the previous results obtained by [69] with KARAT, which is a FDTD-based EM-PIC algorithm. Fig. 11 shows $J_{2D,\text{max injection}}$ versus $L_{z}/\rho_{o}$ 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 $L_{e}/\rho_{o}$ increases, the current density converges to the limit of Langmuir-Blodgett’s law as expected.

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 $\rho$ 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.

5 Numerical examples

In this section, we present simulations of a relativistic BWO device operating at $\pi$-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.

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 $V_{0}$ between cathode and anode. We choose $-550$ kV which produces the electron beam with mean axial velocity $v_{\text{beam}}=0.877c$ with the width of $2$ mm. Fig. 14 illustrates electric potential distribution (contour plots) and corresponding electric fields (vector plots) by solving Poisson equations. Each super-particle represents $1.495\times 10^{7}$ electrons so that the total injection current is around $1.5$ kA. On average, $102$ macro-particles are assigned to each grid cell. The Debye length $\lambda_{D}$ in our simulations is $10.7258$ mm and much larger than the average grid (edge) length, $l_{av}=1.1123$ 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 $R(z)=\frac{1}{2}(A+B)+\frac{1}{2}(A-B)\cos\left(\frac{2\pi}{C}z\right)$, where $A$ and $B$ are maximum and minimum radii, respectively, and $C$ is the corrugation period. The total number of corrugations along the structure is denoted as $N_{\text{crg}}$. 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 $A=17.1$ mm, $B=12.9$ mm, $C=7.5$ mm, and $N_{\text{crg}}=14.5$ 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 $0.2\lambda_{0}$ where $\lambda_{0}$ is wavelength of the center frequency, the PML reflection coefficient, defined in [65], is as low as $-92$ dB. In the beam focusing system, a static axial magnetic field is applied over the SWS region. We enforce axisymmetric fields at the $z$-axis, by applying a perfect magnetic conductor (PMC) boundary condition there (note that only the $E_{z}$ component is present on axis). For the particles incident on the axis we use a perfectly reflecting boundary condition.

5.1.1 System performance

Fig. 13 illustrates a snapshot of the electron beam at $t=21.50$ 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 $v_{\text{p}}/c$. 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 $24.00$ ns. The particles decelerated from the initial velocity ($0.877c$ m/s) dominate the accelerated particles, so that a net loss of the beam kinetic energy leads to amplification of the $\text{TM}_{01}$ mode.

Fig. 17 presents a vector plot of self-fields generated by the electron beam at $t=76.00$ ns (in steady state). Coherent Cerenkov beam-wave interactions give rise to a strong $\text{TM}_{01}$ mode that may be observed within the SWS region.

The self-fields at the output port were analyzed in time and frequency domains, as shown in Fig. 18a and Fig. 18b, respectively. It can be seen in Fig. 18a 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 $35$ 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 $f_{\text{osc}}=15.575$ GHz, as shown in Fig. 18b.

In order to verify charge conservation, we plot the normalized residuals (NR) versus the nodal indices for the DCE and DGL [47] in Figs. 19a and Fig. 19b, respectively. These residuals, evaluated for each time step $n+\frac{1}{2}$ or $n$ and node $i$, are defined as

[TABLE]

Fig. 19a shows that NR levels for DCE remain fairly low at all nodes and very close to the double precision floor (below $10^{-15}$). 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 $10^{-10}$. NR levels for DGL are also distributed around the double precision floor as shown in Fig. 19b, which means that charge conservation is maintained within round-off errors.

The transverse dynamics of an electron beam in the cross-section plane are sensitive to gyroradius value $\rho_{g}$ which depends on the BFS magnetic field strength. Fig. 20 shows 3D electron velocity plots for a weaker BFS magnetic field reduced from $5$ T to $0.5$ T. RGB colors again indicate particle velocity magnitudes. The weaker axial magnetic field yields a relatively larger $\rho_{g}$ in the polar ($\rho\phi$) 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: ($i$) a coarse mesh, ($ii$) a double-refined mesh, ($iii$) a quadruple-refined mesh, ($iv$) a octuple-refined mesh, ($v$) a coarse mesh with staircased boundary, ($vi$) a double-refined mesh with staircased boundary, ($vii$) a quadruple-refined mesh with staircased boundary, and ($viii$) a octuple-refined mesh with staircased boundary. Fig. 21 depicts one period of the SCSWS boundary as rendered by the different meshes.

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.

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

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. 22a converge much faster than those of meshes with staircased boundaries in Fig. 22b 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 $\pi$-point frequency101010$\pi$-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 [70] 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 $\pi$-point of the modal fields. The underestimation of the $\pi$-point frequency in the staircased boundary is shown by Fig. 24, which depicts the dispersion diagrams for the $\text{TM}_{01}$ mode of each SWS. The analytic dispersion relation for the $\text{TM}_{01}$ mode is obtained based on Floquet’s theory by considering harmonics up to $6^{\text{th}}$ order [71]. It is clearly seen that the $\pi$-point frequency decreases as the staircasing error become significant. On the other hand, the cases without staircasing errors rapidly converge to the analytic prediction.

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 $\text{TE}^{\phi}$-polarized fields in the meridian ($\rho z$) 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 $\mathcal{E}-\mathcal{B}$ FETD scheme. A local and explicit time update is made possible by employing the SPAI approach [47]. 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 [41]. 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.

Acknowledgment

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 $p$-forms is that they admit a natural “pairing” with the $p$-cells of the mesh, where $p$ refers to the dimensionality i.e. $p=0$ refers to nodes, $p=1$ to edges, and so on [37]. Computationally, the pairing operation between an $i$-th $p$-cell of the grid $\sigma_{p}^{i}$ and the Whitney form $w^{p}_{j}$ associated with the $j$-th $p$-cell is effected by the integral below and yields [37, 40]

[TABLE]

for $p=0,1,2,3$ in the 3D space and where $\delta_{i,j}$ is the Kronecker delta. The generalized Stokes theorem [37, 40, 72] states

[TABLE]

where $\partial$ 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 $2$-cells into Faraday’s law on the primal mesh gives

[TABLE]

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

[TABLE]

Substituting (5) and (6) into (50) and using $\left(\partial\sigma^{i}_{2}\right)_{1}=\sum_{k=1}^{N_{1}}C_{ik}\sigma^{k}_{1}$ where $C_{ik}$ is an incidence matrix element taking a value in the set of $\left\{-1,0,1\right\}$ [37, 40, 62, 63], we obtain

[TABLE]

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

[TABLE]

for $i=1,...,N_{2}$. (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 [39, 46, 58, 59, 64] and considering a (doubly) artificial inhomogeneous medium at the meridian plane as done in (12) and (14), we can rewrite (21) and (22) as

[TABLE]

where $\vec{W}_{i}^{p}$ is a vector proxy for Whitney $p$-form associated with the $i$-th $p$-cell. The vector proxies for Whitney 1- and 2-forms are written as [41]

[TABLE]

where $i_{a}$, $i_{b}$, and $i_{c}$ denote indices for grid nodes belonging to the $i$-th $p$-cell, for $p=1$ or $2$, and $\lambda$ denote the barycentric coordinates of each node. Since Whitney forms have compact support, we can express the global discrete Hodge matrix in terms of sum of local matrices for the $i$-th triangle as

[TABLE]

where $\Delta_{i}$ is the area of $i$-th triangle, $t_{a}=\frac{\vec{\nabla}\lambda_{1}\cdot\vec{\nabla}\lambda_{1}}{5}$, $t_{b}=\frac{\vec{\nabla}\lambda_{1}\cdot\vec{\nabla}\lambda_{2}}{5}$, $t_{c}=\frac{\vec{\nabla}\lambda_{2}\cdot\vec{\nabla}\lambda_{2}}{5}$, and $\rho_{1}$, $\rho_{2}$, $\rho_{3}$ denote the $\rho$ coordinate in polar coordinate systems with origin at the grid nodes belonging to the $i$-th $2$-cell (triangle), and

[TABLE]

The elements of $\left[\star_{{\mu}^{-1}}\right]$ in 2D take the form of

[TABLE]

The sparsity patterns for discrete Hodge matrices of the mesh exemplified in Fig. 26 are provided in Fig. 25. Due to finite support of Whitney forms, these matrices are sparse. In the 2-D case, each row of the permittivity Hodge matrix has only 5 nonzero elements (except for the rows associated to boundary edges) and the permeability Hodge matrix is a diagonal matrix.

Appendix C Barycentric dual lattice relations

The barycentric dual lattice [40, 55] has a similar contraction identity to (47), given by

[TABLE]

where $n$ is the ambient space dimension and $p$ is the grid element dimension. Since we consider Hodge duals of electric current and charge densities in the primal mesh, we need to express them in dual formulation to be used in discrete Ampere’s law. The electric current density can be expressed by

[TABLE]

In order to obtain the discrete representation for Ampere’s law, we combine the pairing $\left<\tilde{\sigma}_{1}^{i},\left(\cdot\right)\right>$ with Ampere’s law. The electric current density term become

[TABLE]

Therefore, we can write above in matrix form as

[TABLE]

where $\left[\mathbb{J}\right]$ is a column vector with all DoFs for $\mathcal{J}$ expanded in terms of Whitney 1-forms on the dual mesh, and $\left[\mathcal{I}\right]$ is the identity matrix. A similar procedure can be done for the electric charge density in Gauss’ law:

[TABLE]

where $\left[\mathbb{Q}\right]$ is a column vector with all DoFs for $\mathcal{Q}$ expanded in terms of Whitney 2-forms on the dual mesh.

Appendix D Structure of the incidence matrices

Incidence matrices represent the spatial differential operators on the mesh (at the discrete level), distilled from the metric structure (i.e., the discrete exterior derivative) [37, 40]. Since the discrete exterior derivative can be seen as the dual of the boundary operator acting on the mesh elements by the generalized Stokes theorem, Eq. (48), incidence matrices encode the relationship between each element of the mesh and its boundary elements (say, between an edge and its boundary nodes, between a face element and its boundary edges, and so on).

In order to illustrate the structure of the incidence matrices, we consider a simple mesh with perfect magnetic conductor (or free edges) boundaries as depicted in Fig. 26. Red-colored numbers denote nodal indices, black-colored numbers edge indices, and green-colored number face indices. Intrinsic edge orientation is defined by ascending index order of the two nodes associated with any given edge. The sign the corresponding $\left[\mathcal{D}_{\text{curl}}\right]$ matrix elements is determined by comparing the intrinsic orientation of each edge with the curl and divergence operator directions indicated by the blue and beige thick arrows in Fig. 26: if they are opposite, the element is $-1$, otherwise it is $+1$. For example, if we consider $\left[\mathcal{D}_{\text{curl}}\right]$, of size $N_{2}\times N_{1}$, there are three edges wrapping face number 6: edges 8, 9, and 20. As a result, $\left[\mathcal{D}_{\text{curl}}\right]_{6,8}=-1$, $\left[\mathcal{D}_{\text{curl}}\right]_{6,9}=+1$, and $\left[\mathcal{D}_{\text{curl}}\right]_{6,9}=-1$. Furthermore, $\left[\mathcal{D}_{\text{curl}}\right]_{6,k}=0$ for all other edges $k$. This is represented in Fig. 27a, which shows the entire $\left[\mathcal{D}_{\text{curl}}\right]$ for this mesh.

Likewise, if we consider $\left[\tilde{\mathcal{D}}_{\text{div}}\right]$, of size $N_{0}\times N_{1}$, there are six edges connected to node 5: edges 7, 10, 12, 13, 14, and 15. The corresponding elements are $\left[\mathcal{D}_{\text{div}}\right]^{T}_{5,7}=-1$, $\left[\mathcal{D}_{\text{div}}\right]^{T}_{5,10}=-1$, $\left[\mathcal{D}_{\text{div}}\right]^{T}_{5,12}=+1$, $\left[\mathcal{D}_{\text{div}}\right]^{T}_{5,13}=+1$, $\left[\mathcal{D}_{\text{div}}\right]^{T}_{5,14}=+1$, and $\left[\mathcal{D}_{\text{div}}\right]^{T}_{5,15}=+1$. This is represented in Fig. 27b.