Mixed Finite Element Simulation for Solving Eigenmodes of Cavity Resonators Filled With Both Electric and Magnetic Lossy, Anisotropic Media

This article presents a mixed finite element method to find the resonant frequencies of 3-D closed cavity resonators filled with both electric and magnetic lossy, anisotropic media. By introducing a dummy variable with zero value in the 3-D linear vector Maxwell’s eigenvalue problem for the electric field, the divergence-free condition for electric flux density is enforced in a weak sense. In addition, by introducing a dummy variable with an arbitrary constant in the 3-D linear vector Maxwell’s eigenvalue problem for the magnetic field, the divergence-free condition for magnetic flux density is enforced in a weak sense. At the same time, we prove that the eigenmodes of the second-order linear Maxwell’s eigenvalue problems with introduced the dummy variables are equivalent to the ones of the original second-order Maxwell’s eigenvalue problems. As a consequence, the novel method of introducing dummy variables can be free of all the spurious modes in solving eigenmodes of 3-D closed cavity problems. Finally, two numerical experiments are carried out to show that the numerical eigensolver supported by this article can eliminate all the spurious modes, including spurious zero modes.


I. INTRODUCTION
It is well known that the microwave resonant cavity has a wide range of applications in microwave engineering, such as microwave oven, microwave signal source, microwave filter and so on. For example, Zhang and Lin [1] designed a dielectric-loaded bandpass filter, which will have some potential applications in 5G mobile communication systems and microwave radar communication systems. The resonant frequency is a key parameter in microwave resonator, so finding several resonant eigenmodes of 3-D closed cavity problems is one of the classic problems in computational electromagnetics.
The divergence condition supported by Gauss's law is critical to numerical simulations of electromagnetic problems.
The associate editor coordinating the review of this manuscript and approving it for publication was Guillaume Parent .
Jiang et al. [2], [3] pointed out that the divergence equations in Maxwell equations were not redundant for transient and time harmonic problems and it was very important to solve Maxwell equations. Note that the computational domains considered in [2] are homeomorphic to the sphere in R 3 , but the computational domain in this article does not have this limitation. Wu and Jiang [4] utilized the least-squares finite element method to solve the electromagnetic scattering problems and directly incorporated two divergence equations in the discretization process. In [4], nodal-based functions are employed in the least-squares finite element method, but nodal basis functions are abandoned in this article and edge element basis functions of the lowest order are used instead. It is known spurious modes will appear in numerical results provided that numerical methods cannot preserve real physical properties of the electromagnetic field in electromagnetic eigenvalue problems correctly. The divergence-free condition supported by Gauss's law is critical to numerical simulations of electromagnetic eigenvalue problems. If the divergence-free condition is artificially ignored in a numerical method, a large number of spurious zero modes will be introduced in the numerical results. For instance, some early articles [5]- [7] about employing edge elements to solve 3-D closed cavity problems do not enforce the divergencefree condition, as a consequence, numerical algorithms given in these articles do not remove spurious zero modes. For numerical simulations of electromagnetic waveguide problems, Rahman and Davies [8], Winkler and Davies [9], and Kobelansky and Webb [10] applied a penalty function to enforce the divergence-free condition, however, this numerical method did not eliminate all the spurious modes in waveguide problems. Therefore, it is important to select appropriate computational methods to solve electromagnetic eigenvalue problems.
For 3-D closed resonant cavity problems filled with lossless media, Venkatarayalu and Lee [11] enforced the divergence-free condition based on tree-cotree partitioning technique and Jiang et al. [12] utilized a Lagrange multiplier to enforce the divergence-free condition in a weak sense. Consequently, the numerical methods given in [11], [12] can remove all the spurious modes, including spurious zero modes. Jiang et al. [13] applied mixed finite element method (MFEM) to solve the 3-D resonant cavity problem filled with electric and/or magnetic lossless media, excluding the case that both electric and magnetic lossy media. For all of these cavity problems, the MFEM can remove all the spurious modes, including spurious zero modes, because MFEM employs a Lagrange multiplier to enforce the divergence-free condition in a weak sense. However, the MFEM given in [12], [13] can not be applied to solve closed resonant problems with both electric and magnetic lossy media, because there are some difficulties in establishing appropriate mixed variational forms. Based upon the idea of the MFEM, Liu et al. [14] proposed a two-grid method to solve the abovementioned 3-D cavity problems. Moreover, the method given in [14] can remove all the spurious modes. For the 3-D cavity problem filled with fully conductive media in the whole cavity, Jiang et al. [15] applied edge element of the lowest order and the standard linear element to solve this problem successfully. It is worthwhile to point out that the medium considered in [15] is only conductive loss, not dielectric loss.
In 2020, Jiang and Liu [16] gave three numerical eigensolvers to solve the 3-D resonant cavity problem. In [16], the projection method can solve 3-D resonant cavity problem filled with both electric and magnetic lossy media. However, this method needs to get the null space of a given sparse matrix. The projection method utilizes a complete singular value decomposition to find a set of normalized orthogonal bases of a null space in order to maintain the numerical stability of the numerical algorithm. It is known that the computational cost of the complete singular value decomposition is very huge, this implies that the projection method given in [16] is not efficient. For the 3-D resonant cavity problem filled with both electrical and magnetic lossy, anisotropic media, the penalty method supported in [16] cannot remove all the spurious modes. Furthermore, the augmented method given in [16] cannot guarantee that the numerical modes obtained by the augmented method are all physical. Therefore, numerical methods given in [16] cannot solve 3-D closed resonant cavity problems with both electric and magnetic lossy media very well.
For these reasons above, this article continues to investigate the 3-D resonant cavity problem filled with both electrical and magnetic lossy, anisotropic media. In this article, on one hand, by introducing a dummy variable with zero value in the 3-D vector Maxwell eigenvalue problem for the electric field, the divergence-free condition for electric flux density is enforced in a weak sense. On other hand, by introducing a dummy variable with an arbitrary constant in the 3-D vector Maxwell eigenvalue problem for the magnetic field, the divergence-free condition for magnetic flux density is enforced in a weak sense. Moreover, it is theoretically proved that the novel method of introducing dummy variables can be free of all the spurious modes, including spurious zero modes. In this article, the method for the constraint of the divergencefree condition is generally consistent with the method used in literature [2].
The outline of this article is as follows. The well posed governing equations for 3-D closed cavity problems filled with both electric and magnetic lossy media are given in Section II. In Section III, we present two mixed variational forms and provide the MFEM to solve them. In Section IV, two numerical experiments are carried out to verify that the numerical method supported by the article is free of all the spurious modes, including spurious zero modes.

II. THE WELL POSED GOVERNING EQUATIONS OF 3-D RESONANT CAVITY PROBLEMS
Let be the three-dimensional geometry of a given cavity resonator, which is usually a bounded domain in R 3 . Note that the geometry may not be are homeomorphic to the sphere in R 3 . The boundary of the cavity resonator is denoted by ∂ . Letn be the unit outward normal vector on the boundary ∂ . The permittivity and permeability in vacuum are denoted by 0 and µ 0 , respectively. This article mainly focuses on establishing a numerical method to find the resonant eigenmodes of the resonant cavity filled with both electric and magnetic lossy, anisotropic media. This shows that the relative permittivity tensor r and the relative permeability µ r of an anisotropic medium are not Hermitian [17], that is to say that both † = and µ † r = µ r are valid, where the symbol † stands for complex conjugate transpose.
According to the classic electromagnetic theory [18], it is known that the governing equations of 3-D closed cavity problems are the source-free Maxwell's equations of the first order. In addition, the governing equations of the 3-D closed cavity problem include perfectly electric conductor (PEC) boundary condition. It is very difficult to solve this first order system numerically because of the constraint of divergencefree conditions. On one hand, eliminating the magnetic field H from the source-free Maxwell's equations of the first order, one can obtain a linear PDE system of the second order with respect to the electric field E, and that is where = ω 2 0 µ 0 is the square of the wavenumber in vacuum. Note that the eigenvalue consists of countable complex numbers since the material in the cavity is both electric and magnetic lossy. In finite element simulation, the PEC boundary condition (1c) is an essential boundary condition. Once the eigenmode ( , E) with = 0 is solved in (1), then the corresponding magnetic field H is given by On the other hand, eliminating the electric field E from the source-free Maxwell's equations of the first order, one can obtain a linear PDE system of the second order with respect to the magnetic field H , and that is Note that the boundary condition (2c) can be deduced from PEC boundary condition (1c). In finite element simulation, two boundary conditions (2c) and (2d) in PDEs (2) are two natural boundary conditions. Once the eigenmode ( , H ) with = 0 is solved in (2), then the corresponding electric field E is given by It is easy to know that the nonzero eigenmodes between PDEs (1) and PDEs (2) are equivalent if the electromagnetic field is at least twice differentiable. However, when the geometry of the resonant cavity has a very complex geometric topology such that both PDEs (1) and PDEs (2) have physical zero modes, these physical zero modes between PDEs (1) and PDEs (2) are not equivalent if the closed cavity is filled with lossless media. For details, please see [19].
It is clear that the divergence-free condition (1b) can be deduced from the equation (1a) when = 0 in (1a) and the divergence-free condition (2b) can be deduced from the equation (2a) when = 0 in (2a). If the divergence-free conditions (1b) and (2b) are ignored in PDEs (1) and (2), respectively, then many spurious zero eigenvalues will be introduced in these two physical models. Furthermore, the algebraic multiplicity of these zero eigenvalues is infinite and these zero eigenvalues are lack of physical significance. Therefore, PDEs (1) and PDEs (2) must include the divergence-free condition in order to remove nonphysical zero modes. From here it can be seen that the divergencefree condition is very important for the resonant cavity problem. The physical modes in resonant cavity must satisfy the divergence-free condition, however, the nonphysical modes in resonant cavity don't satisfy the divergence-free condition.
It is known that a three-dimensional vector function has three scalar function components, therefore there are respectively three unknown scalar functions in PDEs (1) and PDEs (2), however, there are respectively four single equations in (1a)-(1b) and (2a)-(2b). That is to say that PDEs (1) and PDEs (2) are two linear overdetermined systems. Based on this fact, two new scalar functions p and q are respectively introduced in (1a) and (2a), such that PDEs (1) and PDEs (2) are well posed.
The well posed eigenvalue problem associated with PDEs (1) is as follow: where ( , E, p) with E = 0 is unknown and α is an arbitrarily given nonzero constant.
The well posed eigenvalue problem associated with PDEs (2) is as follow: where β is an arbitrarily given nonzero constant. Here, our aim is to find ( , H , q) with H = 0 such that PDEs (4) are valid. From the consideration in numerical rounding error, how to select the appropriate constant α in (3a) and β in (4a) are given in the next section. Now it can be proved that the unknown scalar function p in PDEs (3) is a dummy variable with zero value. As a matter of fact, by taking the divergence at the both of (3a) and then utilize (3b) and (3d), at last we get the following PDE: It follows that p in (5a) is a harmonic function. According to the classic extremum property of the harmonic function [20], it is easy to obtain the conclusion that the function p vanishes in the whole domain . As a consequence, the solution ( , E, 0) to (3) is the solution ( , E) to (1). Obviously, the solution ( , E) to (1) is also the solution ( , E, p) to (3) as long as we take p = 0 in . Therefore, the solution between (1) and (3) is equivalent. Next it can be proved the scalar function q in PDEs (4) is a dummy variable with any constant value. As a matter of fact, by taking the divergence at the both of (4a) and then use (4b) and (4e), at last we get the following PDE: Using the first Green's theorem of scalar function, one can obtain Substituting (6a) and (6b) into (7) gives ∇q · ∇qd = 0.
The above equation (8) implies ∇q = 0 in . As a result, the scalar function q is an arbitrary constant. That is to say that the solution ( , H , C) to (4) is the solution ( , H ) to (2), where C is an arbitrary constant. Obviously, the solution ( , H ) to (2) is also the solution ( , H , q) to (4) as long as we take q = C in . Therefore, the solution between (2) and (4) is also equivalent.

III. MIXED FINITE ELEMENT METHOD DISCRETIZATION
In this section, the MFEM in computational mathematics is applied to solve PDEs (3) and PDEs (4). In order to give the mixed variational forms associated with PDEs (3) and PDEs (4), at first, we need to introduce the following Hilbert spaces: 3 : ∇ × F ∈ (L 2 ( )) 3 H 0 ( ) = F ∈ H( ) :n × F = 0 on ∂ Next, let us introduce the mixed variational forms associated with PDEs (3) and PDEs (4).

A. MIXED VARIATIONAL FORMS
Using Green formulas, the mixed variational form associated with PDEs (4) reads as: Note that * stands for taking the complex conjugate of a complex vector. The mixed variational form associated with PDEs (3) reads as: is a good approximation of H 1 0 ( ). The definition of W h 1 is given in [12]. For details, please see [12].
Restricting the mixed variational form (9) on W h 1 × S h 1 , one can get the discrete mixed variational form associated with (9), and that is Restricting the mixed variational form (10) on W h 2 × S h 2 , one can get the discrete mixed variational form associated with (10), and that is Note that ( h , H h , q h ) in (11) and ( h , E h , p h ) in (12) are an approximation of the exact solution ( , H , q) in (9) and ( , E, p) in (10) respectively. Next, the MFEM is applied to solve the discrete mixed variational forms (11) and (12).   and W h 1 is given in [12]. It's important to point out these three boundary conditions (4c)-(4e) are three natural boundary conditions in FEM. Hence, the boundary conditions (4c)-(4e) do not need a special treatment in the numerical implementation of MFEM. However, the boundary conditions (3c)-(3d) need a special treatment in the numerical implementation of MFEM, since the boundary conditions (3c)-(3d) are two essential boundary condition in FEM. Here, we only list the matrix eigenvalue problem corresponding to the mixed variational form (11), and that is

B. GENERALIZED EIGENVALUE PROBLEM
where According to actual calculations, we find that the maximum values of the absolute value of the element in matrices A and B are not in the same magnitude order. In order to decrease the rounding error in numerical calculation, we take β = (A)/ (B) in (13), where (A) and (B) are the maximum values of the absolute value of the element in matrices A and B, respectively. In (12), the constant α is also selected in this way.
After solving the generalized matrix eigenvalue problem (13) by the implicitly restarted Arnoldi methods in package ARPACK [22], we can obtain several numerical physical eigenvalues.

IV. NUMERICAL EXPERIMENTS
In this section two numerical experiments are carried out to show that the numerical method given in the article can remove all the spurious modes, including spurious zero modes.
Because it is very difficult to solve the resonant cavity problem filled with anisotropic media by an explicit analytical method, numerical methods must be applied to solve this problem. How to measure whether the numerical result of a numerical solution is correct or not? One way is in terms of the numerical examples in the existing references as benchmarks. The other way is to use commercial software such as COMSOL Multiphysics, ANSYS and then compare their numerical results with those of our numerical method. Here, COMSOL Multiphysics is utilized to simulate the 3-D closed resonant cavity problems. The computational model in COMSOL is as follows: It is clear that PDEs (14) ignores the divergence-free condition (1b). As a consequence, there are many spurious zero eigenvalues in the numerical eigenvalues from COMSOL Multiphysics.
Our MFEM codes are written in MATLAB, and they run on a computer with 8-cores Intel Core i9-11900K 3.5GHz CPU and 128GB-RAM. The mesh data of our MFEM codes are from COMSOL. In the numerical simulations of COMSOL, the FEM based on the edge element of the lowest order is applied to simulate PDEs (14).
In order to distinguish the numerical eigenvalues from (11) and (12). The numerical eigenvalues from (11) and (12) are denoted by h (cu, H ) and h (cu, E), respectively. The numerical eigenvalues obtained by the projection method in [16] to solve PDEs (2) are denoted by h (pr, H ). In addition, the numerical eigenvalues obtained by COMSOL Multiphysics to solve PDEs (14) are denoted by h (COMSOL, E).

A. CYLINDRICAL CAVITY
The radius and height of cylindrical cavity is denoted by r and b, respectively. Set r = 0.2 m and b = 0.5 m. Suppose that the relative permittivity and permeability tensor of the anisotropic medium in the whole cylindrical cavity are It is clear that both † r = r and µ † r = µ r are valid, and this implies that the anisotropic medium is both electric and magnetic lossy.
In [16], the penalty method, augmented method and projection method are used to solve the above numerical example and the solved model is just PDEs (2). Among these three computational methods, only the projection method does not introduce any spurious modes. Hence, in Table 1, we only list the numerical eigenvalue h (pr, H ) of the projection method. The time t 1 and t 2 are obtained by h (cu, E) and h (cu, H ) using the computational methods in the current article respectively. The time t 3 are obtained by the projection method in [16]. For simplicity, we only list the numerical eigenvalue corresponding to dominate mode in Table 1.
As shown in Table 1, h (cu, E) and h (cu, H ) are approximately equal to h (pr, H ). t 1 ≈ t 2 < t 3 implies that the computational method supported by this article is more efficient than the projection method in [16]. The reason for the low computational efficiency of the projection method is that this projection method uses a complete singular value decomposition algorithm to find the dense basis vector of the null space of a rectangular matrix. In contrast, the algorithm in this article does not involve the calculation of dense matrices. It is worth pointing out that the algorithm provided by this article can remove all the spurious modes when it is used to solve the resonant cavity problem filled with a both electric and magnetic lossy, anisotropic medium.

B. TORUS CAVITY
The parametric definition of a torus is as follows: where ρ 1 and ρ 2 are major and minor radiuses, respectively. Set ρ 1 = 0.8 m, ρ 2 = 0.4 m. Assume that the torus cavity is filled with a homogeneous, both electric and magnetic lossy, anisotropic medium. The medium parameter of the material in this torus cavity is given by This torus cavity problem is simulated by the numerical methods in this article and COMSOL Multiphysics. The numerical eigenvalues h (cu, E), h (cu, H ) and h (COMSOL, E) associated with the dominant mode are listed In Table 2. It can be observed In Table 2 that h (cu, E) ≈ h (cu, H ) ≈ h (COMSOL, E). This implies that the implementation of our MATLAB codes is correct. It is worthwhile to point out that there are not any numerical spurious results from our MATLAB codes, however, there are many spurious zero eigenvalues from COMSOL.
In addition, according to the above two numerical examples A and B, we do find that the entries from ζ in (13) are all the same and p h ≈ 0 in (12a) is valid. These verify the correctness of the theory in Section II.

V. CONCLUSION
By introducing dummy variables, this article enforces the divergence-free condition supported by Gauss's law in electromagnetics in a weak sense. As a consequence, this novel method can be free of all the spurious modes in solving the resonant cavity problem filled with both electric and magnetic lossy, anisotropic media.
In the current article, the MFEM based upon the edge element of the lowest order and the standard linear element is to find the eigenmodes of the resonant cavity. It is known that the convergence rate of this type of MFEM is very slow, therefore, we need to accelerate the convergence rate of the MFEM. This leads to the MFEM based on higher order basis functions. In the future, the MFEM based upon higher order basis functions is applied to solve 3-D classic resonant cavity problem. Notice that this MFEM also has an ability of removing all the spurious modes.