# A Finite Element-Based Characteristic Mode Analysis

Konstantinos D. Paschaloudis, Constantinos L. Zekios, Stavros Georgakopoulos, George A. Kyriacou

► To cite this version:

Konstantinos D. Paschaloudis, Constantinos L. Zekios, Stavros Georgakopoulos, George A. Kyriacou. A Finite Element-Based Characteristic Mode Analysis. Ieee Open Journal of Antennas and Propagation, 2022, 3, pp.287-303. 10.1109/OJAP.2022.3150594. hal-03633546

HAL Id: hal-03633546

https://hal.science/hal-03633546

Submitted on 6 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Received 18 January 2022; revised 31 January 2022; accepted 6 February 2022. Date of publication 10 February 2022; date of current version 18 February 2022.

Digital Object Identifier 10.1109/OJAP.2022.3150594

# A Finite Element-Based Characteristic Mode Analysis

KONSTANTINOS D. PASCHALOUDIS<sup>1</sup> (Member, IEEE), CONSTANTINOS L. ZEKIOS<sup>2</sup> (Member, IEEE), STAVROS V. GEORGAKOPOULOS<sup>2</sup> (Senior Member, IEEE), AND GEORGE A. KYRIACOU<sup>3</sup> (Senior Member, IEEE)

CORRESPONDING AUTHOR: C. L. ZEKIOS (e-mail: kzekios@fiu.edu)

This work was supported by the Air Force Office of Scientific Research under Grant FA9550-18-1-0191 and Grant FA9550-19-1-0290.

**ABSTRACT** A novel Green's function-free characteristic modes formulation is introduced in this work. The desired impedance or admittance matrix is obtained utilizing and appropriately modifying the versatile finite element method. For this purpose, the generalized eigenvalue problem of the electric or magnetic field vector wave equation is formulated. In the case of the electric field wave equation, using the Schur complement, the system is reformulated and expressed only in terms of the tangential electric field over the radiating apertures, retaining the equivalent magnetic currents. Similarly, in the case of the magnetic field wave equation, the electric current density on radiating metallic surfaces is isolated using the Schur complement. In both cases, the obtained matrix is split into its real and imaginary part to yield the characteristic modes eigenvalue problem. Key advantage of the proposed formulation is that it does not require the evaluation of Green's function, thereby the study of any arbitrarily shaped, multilayered geometry loaded with anisotropic and inhomogeneous materials is feasible. To prove the validity of the proposed methodology various classical structures, with both homogeneous, and inhomogeneous and anisotropic materials, published in the bibliography are studied. Both the eigenvalues and eigenvectors compared with the published results show good agreement.

**INDEX TERMS** Characteristic modes, characteristic mode theory, finite element method, Green's function, anisotropic materials, inhomogeneous structures.

## I. INTRODUCTION

CHARACTERISTIC mode (CM) theory or the theory of characteristic modes (TCM) has become extremely popular in the last 10 years. From the late 1960s that it was initially introduced with a series of works from Garbacz *et al.* [1] and Chang and Harrington [2], until now many advances have occurred [3]. At its original form the electric field integral equation (EFIE) was used for the modal analysis of arbitrarily shaped perfect electric conductive (PEC) structures [4], [5]. For the case of dielectric and magnetic bodies a volume integral equation (VIE) formulation [6] was initially introduced. Although the VIE formulation was proven to be accurate, it requires large number of unknowns when the electric size of the structure increases. To overcome

this problem an equivalent surface integral equation (SIE) formulation was introduced in [2]. However, the formulation of [2] produces spurious modes that do not exist in the VIE-based formulation of [6]. These spurious modes have nothing to do with the complexity of the structure. They also exist for any dielectric/magnetic body. Furthermore, when the structures' complexity increases (e.g., multilayer structures with different dielectric properties), the aforementioned technique based on method of moments (MoM) suffers from inaccuracies as discussed in [3], [7]. Also, for simple and small structures, such as a sphere, some SIE-based CM formulations lead to spurious modes. To overcome these problems, significant research, with the initialization of a big number of different formulations that tries to give solution, has been

<footnote n="1">¹Institut d'Electronique et de Télécommunications de Rennes, Université de Rennes 1, 35000 Rennes, France</footnote>

<footnote n="2">²Department of Electrical and Computer Engineering, Florida International University, Miami, FL 33174, USA</footnote>

<footnote n="3">³Department of Electrical and Computer Engineering, Democritus University of Thrace, 67131 Xanthi, Greece</footnote>
employed [8], [9]. In [10], [11] post processing methods are used to remove the erroneous solutions. In [3], [12], [13] the magnetic (electric) current is eliminated from the PMCHWT (Poggio-Miller-Chang-Harrington-Wang-Tsai) equations and the eigenvalue equation is expressed only in terms of electric (magnetic) current. In [14] a modified SIE is introduced with the utilization of the generalized impedance boundary condition, a continuation of the work in [15], [16]. These formulations, although look well-defined with robust physical interpretation, are still found to produce spurious modes [7]. For example, formulations of [14]–[17] are shown to be immune to the spurious modes. Some of these formulations may, however, suffer from spurious PEC-cavity interior resonances, like the classical CM formulation of [4] does [15].

Nonetheless, all the efforts so far, with the exception of Garbacz's work [1] where the modes are evaluated from a T-matrix (transition/perturbation matrix) formulation and the very recent papers of Gustafsson et al. [18], [19], are expressed in MoM-based characteristic modes formulations, in which the existence of Green's functions requires considerably greater analytical effort [20]. Furthermore, the necessity of Green's functions in SIE-based MoM makes the study of complicated structures (e.g., structures filled with complex materials of high inhomogeneity and anisotropy) cumbersome [21]. On the contrary, VIE-based formulation is available for any penetrable material since it requires only Green's functions of the background medium. However, the determination of a volume impedance matrix remains still extremely demanding and thus computationally prohibitive [15]. For this purpose, the current work proposes a novel characteristic modes formulation based on finite element method (FEM). The FEM utilized herein is based on our previous eigen-analysis work [22]. The proposed formulation constitutes a robust numerical method appropriate for the study of any arbitrarily shaped geometry loaded with anisotropic and inhomogeneous materials. The radiation condition is taken into account by applying absorbing boundary conditions (ABCs) instead of PML (Perfect Matched Layer), utilizing our previous approach for the spurious modes suppression that arises due to the irrotational vector space spanned by the edge element basis [22]. Although, MoM is more efficient for scattering problems than FEM, the densification of the mesh near the aforementioned artificial surface provides accurate results [23]–[25].

An initial attempt with FEM-based CM formulation was made in 2015 [26], where the eigenvalues of a patch antenna were presented. A year after, a FE-Boundary Integral (BI) formulation for CM was presented [27], where also the resonant frequencies of a simple dielectric resonator were shown. In the current work the generalized eigenvalue problem of the electric (magnetic) vector wave equation is adopted. Using linear algebra operations the system is effectively expressed only in terms of the desired (electric or magnetic) tangential field, retaining the equivalent currents (magnetic or electric). This approach yields the equivalent magnetic current density

eigen-problem over radiating apertures and alternatively the equivalent electric current density on them. The important real (or physical) electric current density on metallic conductors is also formulated. The obtained matrix as in the classical characteristic mode analysis is split into its real and imaginary part, defining the well-known characteristic modes eigenvalue problem. The proposed formulation is robust and most importantly does not require the evaluation of Green's function. To validate and test the introduced formulation lossless, lossy, and anisotropic geometries are simulated. Section II presents the FEM-based CM formulation. The resulting characteristic eigenvalues and eigencurrents are cross-examined both with analytical and other numerical SIE MoM-based published results. In all the cases the validity of the results is evident.

## II. FORMULATION

The formulation of the characteristic modes based on FEM, now, is divided into two parts. The first part characterizes the electric currents ($J_s$) flowing either upon metallic parts (characterized as perfect electric conductors) or upon dielectric surfaces. The second part characterizes the magnetic currents flowing on the radiating apertures ($M_s$), which play the role of the interface between the studied geometries and the air. In both cases the field equivalence principle is utilized (e.g., [28]) to reformulate the field-based problem of the typical FEM formulation to its equivalent current-based problem of the known MoM formulation. Thus, it is a problem retrieval of impedance (admittance) matrix in which the CM is applied.

Aiming at a general formulation, an inhomogeneous three dimensional (3D) arbitrarily shaped computation domain $\Omega$ loaded with different anisotropic materials of tensor permittivity ($\overline{\epsilon}_r$) and tensor permeability ($\overline{\mu}_r$) is considered as shown in Fig. 1. This domain is source free and is enclosed by a fictitious surface ($S_{\overline{\epsilon}}$) which truncates the unbounded region. The field (electric or magnetic) vector wave equation is used to describe the electromagnetic behavior of this domain [29]:

$$ \nabla \times \overline{\epsilon}_r^k \cdot \nabla \times \mathbf{F}^k - k_0^2 \overline{\nu}_r^k \mathbf{F}^k = 0 \quad (1) $$

where $k = 1, 2$ such that $\mathbf{F}^1 = \mathbf{E}, \mathbf{F}^2 = \mathbf{H}, \overline{\epsilon}_r^1 = \overline{\mu}_r^{-1}, \overline{\epsilon}_r^2 = \overline{\epsilon}_r^{-1}, \overline{\nu}_r^1 = \overline{\epsilon}_r, \overline{\nu}_r^2 = \overline{\mu}_r$.

According to the CM formulation [4], an operator for the electric or magnetic current flowing upon a surface is defined. In the present formulation, the operator of interest is the electric or magnetic field operator applied on the wave equation:

$$ \mathcal{L}\{\mathbf{F}^k\} = \nabla \times \overline{\epsilon}_r^k \cdot \nabla \times \mathbf{F}^k - k_0^2 \overline{\nu}_r^k \mathbf{F}^k \quad (2) $$

Following the standard Galerkin procedure, (1) is written in its weak form [29], [30]:

$$ \iiint_{\Omega} (\nabla \times \mathbf{T}) \cdot \overline{\epsilon}_r^k (\nabla \times \mathbf{F}^k) d\Omega - k_0^2 \iiint_{\Omega} \mathbf{T} \cdot (\overline{\nu}_r^k \mathbf{F}^k) d\Omega \\ \mp jk_0 \tilde{Z}^k \iint_{\Gamma} \mathbf{T} \cdot (\hat{\mathbf{n}} \times \mathbf{G}^k) d\Gamma = 0 \quad (3) $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Arbitrarily shaped 3D open cavity loaded with different inhomogeneous materials denoted by distinct colours.</figcaption></figure>

where $\tilde{Z}^1 = Z_0$, $\tilde{Z}^2 = \frac{1}{Z_0}$, $Z_0$ stands for the free space intrinsic impedance, $G^1 = H$, $G^2 = E$, $k_0$ for the wavenumber, and $T$ for the test function [29]. The surface integral is defined over the surface $\Gamma$ of any element.

The infinite solution domain is truncated with the introduction of the fictitious surface $S_{fi}$. To take into account the radiation conditions, 1st kind absorbing boundary conditions (ABCs) upon the fictitious surface ($S_{fi}$) is utilized. Thus, the weak formulation (3) takes the form:

$$ \iiint_{\Omega} (\nabla \times \mathbf{T}) \cdot \bar{\bar{\tau}}_{r}^{k}(\nabla \times \mathbf{F}^{k}) d\Omega - k_{0}^{2} \iiint_{\Omega} \mathbf{T} \cdot (\bar{\bar{\nu}}_{r}^{k} \mathbf{F}^{k}) d\Omega \\ \pm jk_{0} \iint_{S_{fi}} (\hat{\mathbf{n}} \times \mathbf{T}) \cdot (\hat{\mathbf{n}} \times \mathbf{F}^{k}) dS_{fi} = 0 \quad (4) $$

As explained in our previous work [22], the last term of (4) represents the enforcement of 1st order ABCs over the fictitious surface $S_{fi}$ enclosing the structure (see Fig. 1) in order to truncate the infinite domain. This term is clearly imaginary in contrast with the classical formulation of CM [4], where the radiation leakage is expressed through the real part of the evolved matrices, while the imaginary part expresses the stored energy. Therefore, a multiplication of the whole system with the imaginary unit (j) is needed, since energy storage in Maxwell equations is related to terms including $(j\omega\epsilon)$ or $(j\omega\mu)$, thus, the ABCs' imaginary part corresponds to real power leakage.

Overall, (4) can be expressed in matrix form as:

$$ [\mathbf{S}][f^k] - k_0^2[\mathbf{M}][f^k] \pm jk_0[\mathbf{R}][f^k] = 0 \quad (5) $$

with

$$ [\mathbf{S}] = \iiint_{\Omega} (\nabla \times \mathbf{T}) \cdot \bar{\bar{\tau}}_{r}^{k} (\nabla \times \mathbf{F}^{k}) d\Omega \quad (6) $$

$$ [\mathbf{M}] = \iiint_{\Omega} \mathbf{T} \cdot (\bar{\bar{\nu}}_{r}^{k} \mathbf{F}^{k}) d\Omega \quad (7) $$

$$ [\mathbf{R}] = \iint_{S_{fi}} (\hat{\mathbf{n}} \times \mathbf{T}) \cdot (\hat{\mathbf{n}} \times \mathbf{F}^{k}) dS_{fi} \quad (8) $$

By appropriately manipulating (5) two complementary CM formulations can be derived:

a) a formulation where all the unknowns of the radiating domain are expressed in terms of equivalent magnetic currents (M-formulation), and

b) a formulation where all the unknowns of the radiating domain are expressed in terms of equivalent electric currents (J-formulation).

Notably, formulation (5) is complete including all information to yield either the accurate field distribution or the eigenvalues of the structures. Special attention has to be given to the case of spurious solutions. Spurious solutions is a term used loosely in the area of computational electromagnetics depending on the method of operation. Specifically, in the finite element method the appearance of spurious modes is a common problem in the solution of discretized Maxwell's equations. Specifically, as we showed in our previous work [22], there are certain cases where the divergence-free condition cannot be ensured only through the vector wave equation leading to erroneous solutions that although are mathematically acceptable, they are physically declined. Thus, its explicit enforcement is demanded. Therein, by appropriately modifying the discretized equations we were able to obtain a spurious-free finite element formulation. Using the same formulation here, the solution domain is expected to be also spurious-free in that sense.

In the characteristic mode analysis community, the term of spurious solutions has a different meaning. First, we can find the classification of spurious modes that can be obtained in the sense of erroneous, non-physical solutions like the ones reported in [2], [3], [12], [13]. Second, the “spurious internal resonances” can be observed due to the null-space solution of the impedance/admittance operator [21], [31], like the ones reported in [14] and [15]. Note here that the internal resonances constitute the resonances of the studied structures caused by the internal waves experiencing multiple internal reflections. Even though these modes are natural, and they have a physical meaning, they are defined as spurious ones, since the scope of the characteristic mode analysis is to mainly reveal the radiating modes.

In this work our results involve also internal resonances. The reason behind the occurrence of these solutions is due to the poor enforcement of the radiation condition at the artificial domain truncating surface when solving for open-radiating structures. Notably, these also occur and at the most established electric (EFIE) or magnetic (MFIE) field integral equations solutions. Harrington pointed out in [32] that only formulations that satisfy both tangential electric and tangential magnetic field boundary conditions (at the abject to free space interface) give unique eigensolutions at all frequencies, i.e., free from fictitious internal resonances. An elegant way to suppress these internal resonances is with the incorporation of a Dirichlet-to-Neumann boundary condition technique [33]. However, a task like that is left for future work.

Solving (5) as a generalized eigen-problem is equivalent to the estimation of complex eigenvalues. Fig. 2 demonstrates
<figure><img src="image_3.png" /><figcaption>FIGURE 2. a) Modeling of a microstrip patch antenna with its radiating metallic and dielectric apertures discretized appropriately. With red its metallic (Sj) surface is highlighted, and with green its radiating dielectric apertures (Sj) are noted. b) Evaluation of the real electric current density flowing on the outer surface of a tetrahedron that touches the metallic conductor. c) Definition of the equivalent magnetic current flowing on the outer surface of a tetrahedron that touches one of the radiating dielectric apertures.</figcaption></figure>

a typical patch antenna and its equivalent electric and magnetic current densities evaluated on its radiating metallic and dielectric apertures, respectively. Specifically, the real electric current density flowing on the metallic conductors (Fig. 2b) can be defined through the magnetic field as $J_s = \hat{n} \times H$. Alternatively, equivalent magnetic current densities $M_s = -\hat{n} \times E$ can be defined (sampled) on dielectric interfaces (on the surface of dielectric objects) through the tangential electric field as depicted in Fig. 2c. Hence, the established FEM formulation (5) can be modified to yield characteristic eigen-problems for either equivalent electric, or magnetic current densities. Below both formulations are derived.

### A. M-FORMULATION

Equation (5), after some algebraic manipulations and in the case where $k = 1$, can be expressed as:

$$[\mathbf{L}(k_0)][e] = 0 \qquad (9)$$

This system can be now written in terms of the degrees of freedom (dofs) that belong on the radiating surfaces, or on the surface of dielectric objects. Thus, taking the discretized geometry of Fig. 2 the dofs are split in those that belong on the radiating aperture $S_i$ of interest (s) and the rest ones (r). Essentially, as s are denoted the dofs of the interface between the patch antenna and the air. Likewise, as r are denoted the dofs that belong in the inner and the outer domain of the surface cover of the patch antenna, which are indeed the rest of the dofs of the solution domain.

The tangential electric field on the radiating aperture $S_i$ is directly related to the FEM dofs $e_s$ defined along the element edges lying on $S_i$ as $\mathbf{E}_{tan} = e_s \hat{l}$, where $\hat{l}$ is used to denote the tangential (colinear) unit vector along an edge. Hence, the equivalent magnetic current can be defined through the unit normal $\hat{n}$ of $S_i$ as $M_s = e_s \hat{l} \times \hat{n}$, as shown in Fig. 2c. Notably, $M_s$ lies again on the aperture surface $S_i$ being everywhere

normal to the electric field that generates it. In this manner, the magnetic current eigen-problem can be readily setup either for $E_{tan}$ or $M_s$ as follows:

$$\begin{bmatrix} \mathbf{L}_{ss}(k_0) & \mathbf{L}_{sr}(k_0) \\ \mathbf{L}_{rs}(k_0) & \mathbf{L}_{rr}(k_0) \end{bmatrix} \begin{bmatrix} e_s \\ e_r \end{bmatrix} = 0 \qquad (10)$$

To retain only the dofs of surface $S_i$, the Schur complement [34] is applied on (10):

$$\begin{align} & j[\mathbf{L}_{ss}(k_0) - \mathbf{L}_{sr}(k_0)\mathbf{L}_{rr}^{-1}(k_0)\mathbf{L}_{rs}(k_0)][e_s] = 0 \nonumber \\ & \iff [\mathbf{Y}(k_0)][e_s] = 0 \tag{11} \end{align}$$

The notation $\mathbf{Y}$ is used to denote the admittance properties of the system since the unknowns correspond to the equivalent surface electric field. The multiplication of (11) with the imaginary unit (j) is done to appropriately isolate the radiating terms and express them as real quantities (from the mathematical point of view), following the classical characteristic mode formulation [4]. Namely, in FEM the radiation is solely introduced through the ABC term [35], that is accompanied by the imaginary unit as shown in (3) and (4), and is the only purely imaginary term in a lossless formulation. Therefore, by multiplying with the imaginary unit (j) the real part of the resulting complex matrix $\mathbf{Y}(k_0) = \mathbf{G}_{rad}(k_0) + j\mathbf{B}(k_0)$ is related with the radiated power as in the classical Green function-based formulation [4]. Similarly, the impedance matrix used in the classical formulation of CM for electric currents, flowing upon conducting surfaces, is complex [4]. Furthermore, the matrix $\mathbf{Y}$ is symmetric following the properties of the FEM formulation [36], and for the case where both $\bar{\bar{\epsilon}}_r$ and $\bar{\bar{\mu}}_r$ are symmetric. Using now linear algebra, the Hermitian parts of $\mathbf{Y}$ are:

$$\mathbf{G}_{rad}(k_0) = \frac{1}{2}(\mathbf{Y}(k_0) + \mathbf{Y}^*(k_0)) \qquad (12a)$$

$$\mathbf{B}(k_0) = \frac{1}{2j}(\mathbf{Y}(k_0) - \mathbf{Y}^*(k_0)) \qquad (12b)$$

In the classical CM formulation for lossless structures, as already mentioned, the real part of impedance matrix is related to the radiated power, while the imaginary part represents the stored power. Therefore, the CM eigenvalue problem takes the form:

$$[\mathbf{G}_{rad}(k_0) + j\mathbf{B}(k_0)][e_s] = (1 + j\lambda_n)\mathbf{G}_{rad}(k_0)[e_s] \qquad (13a)$$

$$\mathbf{B}(k_0)[e_s] = \lambda_n \mathbf{G}_{rad}(k_0)[e_s] \qquad (13b)$$

At this point it is important to note the following. The theory of characteristic modes is usually defined with the use of an impedance matrix. Here, as we show with (13) an admittance matrix is used instead. This might raise the question of why an admittance-based eigenvalue problem can derive a CM eigenvalue problem. The answer to this question can be found in the seminal work of Harrington et al., in [6], where there an admittance-based eigenvalue problem is defined for the case of magnetic bodies. Specifically, the theory of characteristic modes for the magnetic bodies, and
the corresponding admittance matrix, are justified as dual to that for dielectric bodies, and the corresponding impedance matrix. For a better understanding of that the interested reader is referred to [28] where the principle of duality can be found.

As in the classical CM formulation, also here, it is important to ensure that the right-hand sides of (13a), (13b) are related only to radiation losses. In any other case, the existence of spurious solutions is inevitable [14]. Thus, lossy structures where the real part of the admittance may include dissipation (finite conductivity or dielectric hysteresis losses) is treated separately according to the recent work of Ylä-Oijala et al. [7], [14], [16], [37].

Eq. (13b) results in the desired eigenvalues ($\lambda_n$) and the corresponding eigenvectors ($e_s$). To finally observe the classical CM formulation where the eigenvectors correspond to currents, and not fields as in (13), the tangential electric field ($e_s$) is replaced by equivalent magnetic currents flowing on the radiating apertures, e.g., [28]. Thus (13b) reads as:

$$ 
\begin{aligned}
\mathbf{B}(k_0)[-\hat{\mathbf{n}} \times e_s] &= \lambda_n \mathbf{G}_{\text{rad}}(k_0)[-\hat{\mathbf{n}} \times e_s] \\
\Leftrightarrow \mathbf{B}(k_0)[m_s] &= \lambda_n \mathbf{G}_{\text{rad}}(k_0)[m_s].
\end{aligned}
\quad (14) $$

## B. J-FORMULATION

The electric current density characteristic eigenproblem can be expressed in terms of:

a) the real current density flowing on the surface of a metallic conductor, or

b) the equivalent electric current flowing on the surface of a dielectric or magnetic object.

In both cases this is defined through the tangential magnetic field $J_s = \hat{\mathbf{n}} \times \mathbf{H}$ as shown in Fig. 2b.

To derive the CM formulation in terms of the equivalent electric currents, flowing upon the metallic surface $S_j$, (5) is written for the case of $k=2$ as:

$$ [\mathbf{L}(k_0)][h] = 0 \quad (15) $$

Splitting now the dofs into those that belong on the radiating aperture $S_j$ of interest ($s$) and the rest ones ($r$), (15) takes the form:

$$ \begin{bmatrix} \mathbf{L}_{ss}(k_0) & \mathbf{L}_{sr}(k_0) \\ \mathbf{L}_{rs}(k_0) & \mathbf{L}_{rr}(k_0) \end{bmatrix} \begin{bmatrix} h_s \\ h_r \end{bmatrix} = 0 \quad (16) $$

In turn, applying Schur complement, the system is expressed in terms of the desired degrees of freedom as:

$$ 
\begin{aligned}
& j[\mathbf{L}_{ss}(k_0) - \mathbf{L}_{sr}(k_0)\mathbf{L}_{rr}^{-1}(k_0)\mathbf{L}_{rs}(k_0)][h_s] = 0 \\
\Leftrightarrow \quad & [\mathbf{Z}(k_0)][h_s] = 0
\end{aligned}
\quad (17) $$

Note that the notation $\mathbf{Z}$ is used to denote the impedance properties of the system since the unknowns correspond to the equivalent surface magnetic field. The multiplication of (17) with the imaginary unit ($j$) is needed once again to ensure the proper correlation between the impedance matrix resulted from FEM and that from classical CM formulation. Also, matrix $\mathbf{Z}(k_0) = \mathbf{R}_{\text{rad}}(k_0) + j\mathbf{X}(k_0)$ is complex

and symmetric, similar to the admittance matrix used in the M-formulation previously. Using linear algebra, the Hermitian parts of $\mathbf{Z}$ are:

$$ \mathbf{R}_{\text{rad}}(k_0) = \frac{1}{2}(\mathbf{Z}(k_0) + \mathbf{Z}^*(k_0)) \quad (18a) $$

$$ \mathbf{X}(k_0) = \frac{1}{2j}(\mathbf{Z}(k_0) - \mathbf{Z}^*(k_0)) \quad (18b) $$

Following the classical CM formulation for lossless structures, the real part of impedance matrix is related to the radiated power, while the imaginary part represents the stored power. Therefore, the related CM eigenvalue problem takes the form:

$$ [\mathbf{R}_{\text{rad}}(k_0) + j\mathbf{X}(k_0)][h_s] = (1 + j\lambda_n)\mathbf{R}_{\text{rad}}(k_0)[h_s] \quad (19a) $$

$$ \mathbf{X}(k_0)[h_s] = \lambda_n \mathbf{R}_{\text{rad}}(k_0)[h_s] \quad (19b) $$

Again, it should be ensured that the right-hand sides of (19a), (19b) involve only radiation losses in order to obtain spurious free solutions [14]. Lossy structures are again treated separately in the next section. To finally observe the classical CM formulation where the eigenvectors correspond to currents and not fields as in (19), the tangential magnetic field ($h_s$) is replaced by real/physical or equivalent electric current density, e.g., [28]. Thus (19b) reads:

$$ 
\begin{aligned}
\mathbf{X}(k_0)[\hat{\mathbf{n}} \times h_s] &= \lambda_n \mathbf{R}_{\text{rad}}(k_0)[\hat{\mathbf{n}} \times h_s] \\
\Leftrightarrow \quad \mathbf{X}(k_0)[j_s] &= \lambda_n \mathbf{R}_{\text{rad}}(k_0)[j_s]
\end{aligned}
\quad (20) $$

## C. CHARACTERISTIC EIGEN-PROBLEM OF LOSSY STRUCTURES

It is important to emphasize once again that for lossy structures the right-hand side of the characteristic problem (14) or (20) must involve only the conductance (resistance) matrix related to radiation $\mathbf{G}_{\text{rad}}(k_0)[\mathbf{R}_{\text{rad}}(k_0)]$, while the conductance (resistance) related to the dissipated power must remain in the left-hand side and represented as $\mathbf{G}_{\text{dissipated}}(k_0)[\mathbf{R}_{\text{dissipated}}(k_0)]$. This requirement is undoubtedly established by Ylä-Oijala et al., in a series of papers, e.g., [7], [14], [16], [37]. The key issue here is to ensure this in a practical situation. An efficient practical approach is to perform a simulation of the lossless structure (e.g., PEC and $\tan\delta=0$) in order to estimate the radiation conductance (resistance) as $\mathbf{G}_{\text{rad}}(k_0) = \text{Re}\{\mathbf{Y}_{\text{lossless}}(k_0)\}|[\mathbf{R}_{\text{rad}}(k_0) = \text{Re}\{\mathbf{Z}_{\text{lossless}}(k_0)\}]$. In turn, the same structure is simulated including all losses, either material losses ($\tan\delta$) and/or losses due to finite metal conductivity, in order to obtain the total conductance (resistance) matrix $\text{Re}\{\mathbf{Y}_{\text{lossy}}(k_0)\} = \mathbf{G}_{\text{rad}}(k_0) + \mathbf{G}_{\text{dissipated}}(k_0)[\text{Re}\{\mathbf{Z}_{\text{lossy}}(k_0)\} = \mathbf{R}_{\text{rad}}(k_0) + \mathbf{R}_{\text{dissipated}}(k_0)]$.

Explicitly, the radiated power $P_{\text{rad}}$ and the related admittance $\mathbf{G}_{\text{rad}}(k_0)$, or resistance $\mathbf{R}_{\text{rad}}(k_0)$ matrix result from the integration of the real part of the Pointing vector over the fictitious surface $S_{fi}$ enclosing the solution domain (Fig. 1). Also, its imaginary part yields the reactive power $P_{\text{reac}}^n$,
e.g., [7]:

$$P_n^{\text{rad}} = \frac{1}{2} \text{Re} \iint_{S_{\text{fi}}} (\mathbf{E}_n \times \mathbf{H}_n^*) \cdot \hat{\mathbf{n}} dS \quad (21a)$$

$$P_n^{\text{reac}} = \frac{1}{2} \text{Im} \iint_{S_{\text{fi}}} (\mathbf{E}_n \times \mathbf{H}_n^*) \cdot \hat{\mathbf{n}} dS \quad (21b)$$

Notably, the corresponding radiation term in the FEM formulation results from the [R] matrix defined in (8) through the absorbing boundary conditions.

Hence, in the proposed formulation, the system (5) is formulated considering a complex valued tensor for the dielectric permittivity $\bar{\epsilon}_r$. In turn, the same procedure described in M- or J-formulation section is followed in order to get the admittance $\mathbf{Y}_{\text{lossy}}$, or impedance $\mathbf{Z}_{\text{lossy}}$ matrix, respectively. Performing the same algebraic computations as before, the Hermitian parts of $\mathbf{Y}_{\text{lossy}}$ and $\mathbf{Z}_{\text{lossy}}$ are:

$$\mathbf{G}_{\text{lossy}}(k_0) = \frac{1}{2} (\mathbf{Y}_{\text{lossy}}(k_0) + \mathbf{Y}_{\text{lossy}}^*(k_0)) \quad (22a)$$

$$\mathbf{B}_{\text{lossy}}(k_0) = \frac{1}{2j} (\mathbf{Y}_{\text{lossy}}(k_0) - \mathbf{Y}_{\text{lossy}}^*(k_0)) \quad (22b)$$

$$\mathbf{R}_{\text{lossy}}(k_0) = \frac{1}{2} (\mathbf{Z}_{\text{lossy}}(k_0) + \mathbf{Z}_{\text{lossy}}^*(k_0)) \quad (22c)$$

$$\mathbf{X}_{\text{lossy}}(k_0) = \frac{1}{2j} (\mathbf{Z}_{\text{lossy}}(k_0) - \mathbf{Z}_{\text{lossy}}^*(k_0)) \quad (22d)$$

Regarding the radiation matrix in the lossy case, this is the same as in the lossless case (12a), (18a), for the M- and J-formulation, respectively, since the consideration of losses does not affect the radiation leakage. Thus, the lossless system of matrices is formulated prior to the lossy system to identify the radiation contribution.

In view of the above, the formulation for the characteristic electric currents in the lossy case is modified as:

$$\mathbf{Z}_{\text{lossy}}(k_0)[j_s] = (1 + j\lambda_n)\mathbf{R}_{\text{rad}}(k_0)[j_s] \quad (23a)$$

$$[\mathbf{X}_{\text{lossy}}(k_0) - j\mathbf{R}_{\text{dissipated}}(k_0)][j_s] = \lambda_n \mathbf{R}_{\text{rad}}(k_0)[j_s] \quad (23b)$$

where $\mathbf{Z}_{\text{lossy}}(k_0) = \mathbf{R}_{\text{lossy}}(k_0) + j\mathbf{X}_{\text{lossy}}(k_0) = \mathbf{R}_{\text{rad}}(k_0) + \mathbf{R}_{\text{dissipated}}(k_0) + j\mathbf{X}_{\text{lossy}}(k_0)$.

Likewise, the characteristic eigen-problem for the equivalent magnetic current densities over the radiating apertures becomes:

$$\mathbf{Y}_{\text{lossy}}(k_0)[m_s] = (1 + j\lambda_n)\mathbf{G}_{\text{rad}}(k_0)[m_s] \quad (24a)$$

$$[\mathbf{B}_{\text{lossy}}(k_0) - j\mathbf{G}_{\text{dissipated}}(k_0)][m_s] = \lambda_n \mathbf{G}_{\text{rad}}(k_0)[m_s] \quad (24b)$$

where $\mathbf{Y}_{\text{lossy}}(k_0) = \mathbf{G}_{\text{lossy}}(k_0) + j\mathbf{B}_{\text{lossy}}(k_0) = \mathbf{G}_{\text{rad}}(k_0) + \mathbf{G}_{\text{dissipated}}(k_0) + j\mathbf{B}_{\text{lossy}}(k_0)$.

In both of these eigen-problems the eigenvalue $\lambda_n$ becomes complex (instead of real in the lossless case), where its real part is related to reactive power while its imaginary part to dissipated power [7]:

$$\text{Re}\{\lambda_n\} = \frac{P_n^{\text{reac}}}{P_n^{\text{rad}}} \quad (25a)$$

$$\text{Im}\{\lambda_n\} = \frac{P_n^{\text{diss}}}{P_n^{\text{rad}}} \quad (25b)$$

Thus, it is clear that $\text{Re}\{\lambda_n\}$ retains the same properties as the solely real $\lambda_n$-eigenvalue of the lossless case. Although, all properties of characteristic modes require some clarification in the lossy case, all these are presented in detail by the research group of Ylä-Oijala [7], [14], [16], [37] and a series of related articles referred therein. Note that both (23a) and (24a) reduce to (20) and (14), respectively, as losses go to zero.

A very important point that needs to be addressed here is regarding the far-field orthogonality of the characteristic modes for lossy bodies. When lossy bodies are studied, the corresponding modes correspond to complex characteristic sources which their fields do not diagonalize the scattering matrix. This results to characteristic fields that are no longer orthogonal [6], [19]. Therefore, in all the cases where losses are introduced, even though we utilize expressions (12a), (12b) and (18a), (18b) that are real and symmetric, the characteristic far-fields are not anymore orthogonal.

### III. NUMERICAL RESULTS

To prove the validity of the proposed CM formulation employing the proposed finite element scheme, several examples are considered and comparisons with either previously published results or numerical results from the use of a commercial software are shown. These structures are studied using one time either the lossless (14) or lossy (24a) equivalent magnetic current eigen-problem, and other time the lossless (20) or lossy (23a) equivalent electric current eigen-problem.

#### A. PEC PLATE

The case of a rectangular PEC plate is initially taken into consideration with dimensions 150 mm × 75 mm × 1 mm, where the height is extremely small compared to length and width, resembling to a flat geometry. It should be mentioned that the chosen dimensions correspond to the size of the display in current mobile phones.

To accurately simulate this geometry, an appropriate discretization is needed on the plate's surface and on the fictitious surface (an appropriately designed parallelepiped is used), where ABCs are enforced. Therefore, the discretization step for the plate's surface is set to $\lambda_{\text{min}}/32$ ($\lambda_{\text{min}}$ is related to the maximum simulation frequency $f_{\text{max}} = 2.5$ GHz), while the fictitious surface is initially placed at a distance equal to $\lambda_{\text{max}}/4$ ($\lambda_{\text{max}}$ is related to $f_{\text{min}} = 0.5$ GHz), and is discretized with $\lambda_{\text{min}}/5$ step. Thus, the resulting mesh consists of 147, 404 tetrahedral cells, while the size of stiffness (S), mass (M) and ABCs (R) matrices is 173, 580 × 173, 580 each. At first sight, the aforementioned values might look prohibitive for numerical simulation if compared with the regular discretization used in SIE formulation. However, the dofs (edges) of interest that
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Eigenvalues of first low-order modes for the PEC plate of dimensions 150 mm × 75 mm × 1 mm. The black dots and the red crosses correspond to the proposed J-formulation for a coarse and fine mesh, respectively, while the blue dashed lines to the published results given in [38].</figcaption></figure>

lie upon the surface of the plate are only 5, 823. Therefore, the size of the impedance matrix Z comprising the CM eigen-problem is reduced to 5, 823 × 5, 823.

Fig. 3 shows the resulting eigenvalues' distributions in the frequency range of 0.5 GHz to 2.5 GHz. The black dots and red crosses stand for the results of the proposed J-formulation, while the blue dashed lines for the results of the SIE MoM-based formulation as published in [38].¹ The agreement between the proposed formulation's results and the published ones is clearly evident around the resonant frequencies, while slight deviations occur away from them. These deviations are diminished by increasing the distance between the structure and the fictitious surface from $λ_{max}/4$ to $λ_{max}/2$, and, increasing the surface's discretization step from $λ_{min}/5$ to $λ_{min}/15$, as presented by the red crosses in Fig. 3. It should be noted that the discretization step at the plate's surface is also reduced to $λ_{mid}/32$ ($λ_{mid}$ is related to $f_{mid} = 1.5$ GHz). In particular, the changes on the mesh discretization increases the number of tetrahedral elements to 234, 209, since the size of S, M, and R matrices rises to 284, 588 × 284, 588, but the final impedance matrix Z obtains a size of 2, 220 × 2, 220, since it only depends on the discretization of the plate's surface. Note that a parametric study was also performed, but omitted here for reasons of brevity, where we individually varied: (a) the distance of the fictitious surface from the PEC plate by keeping the same mesh density on its surface, and (b) the mesh density of the fictitious surface by keeping its distance constant in respect to the PEC plate. The combination of having the fictitious surface in a distance of $λ_{max}/2$ and with a surface mesh density of $λ_{min}$ was chosen as a sweet spot to achieve average errors below 1% in respect to the values obtained from the literature. Further increase of the fictitious surface's mesh density can eliminate any deviation between the results of the proposed FEM-based CM formulation and those published

1. Note that in [38] a conducting sheet of zero thickness is modeled, instead of the 1 mm thick metallic sheet that is modeled here. The introduction of such a thickness, however, is insignificant and it doesn't have any effect on the accuracy of the results.

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Modal significance curves for the modes shown in Fig. 3.</figcaption></figure>

in [38]. However, such changes would increase dramatically the mesh and consequently the time consumption, which now is estimated around 5.9 hours per frequency point for the first mesh and around 9.8 hours per frequency point for the denser one.² Nonetheless, the primary concern in antenna design is around the resonant frequencies, where the associated modes contribute more effectively to radiation [40], and where both the coarse (black dots in Fig. 3) and fine (red crosses in Fig. 3) meshes show great agreement with the reference results (dashed blue line in Fig. 3). This is also shown via the modal significance (MS), a modal parameter which is developed to investigate the modal parameters of the antenna under study, and reads [40]:

$$MS = \left| \frac{1}{1 + j\lambda_n} \right| \qquad (26)$$

The closer to unit the curve of modal significance is, the most effectively the corresponding mode radiates. Fig. 4 illustrates the MS distributions for the PEC plate in the frequency range of 0.5 GHz to 2.5 GHz. The black dots and red crosses stand for the results of the proposed J-formulation for the coarse and finer mesh, respectively, while the blue dashed line represents the published results given in [38]. As it is shown, the results of the proposed method are almost identical with the published ones near the resonance, where the values of MS reach the maximum value, for both meshes. Some deviations are observed for small MS values, which are evidently reduced with the finer mesh density. Therefore, the initial coarse discretization of the fictitious surface gives accurate results within a reasonable period of time as described above.

Fig. 5(d) illustrates the electric eigencurrents for modes 1, 3, 5 (as enumerated in Fig. 3) at 1.5 GHz. Note that the corresponding magnitudes are arbitrary as they refer to

2. Note here that the computation time is significantly higher compared to the surface integral equation-based formulation. The reason of this high computation time is the need to apply the Schur complement, and specifically the inversion of the $L_{rr}$ matrix as appears in both M- and J-formulations, [Eqs. (11) and (17), respectively]. Notably, since all these matrices are obtained from a 3D FEM formulation, they are sparse having a theoretical computational complexity for the sparse matrix factorization of $O(N^2)$ [39].
<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**FIGURE 5.** Electric current for the modes on the PEC plate at 1.5 GHz. (a) Mode 1. (a) Mode 3. (a) Mode 5. (d) Color scale.</figcaption></figure>

eigensolutions. All the eigencurrents are identical with the eigencurrents shown in [38, Fig. 3]. Explicitly, mode 1 is the first fundamental mode of a rectangular PEC plate, mode 3 has a circular weakly radiating current distribution, and 5 is a second-order mode.

<figure><img src="image_10.png" /><figcaption>**FIGURE 6.** Validation on *J*− (red dots) and *M*− (blue dots) formulations against Hu-Wang *J*-formulation (blue solid lines) [41] and Ylä-Oijala *et al.* (black dashed lines) [14]: Magnitudes of the first four eigenvalues with the smallest magnitude for a lossless non-magnetic dielectric sphere with *α* = 4.6 mm, and *ε*<sub>r</sub> = 38.</figcaption></figure>

## B. DIELECTRIC SPHERE

The next study is focused on a simple lossless dielectric sphere with radius $\alpha = 4.6$ mm, $\epsilon_r = 38$ and $\mu_r = 1$. A $\lambda_{min}/64$ discretization step is used for the sphere's surface, while the fictitious surface has a mesh density of $\lambda_{min}/10$. The resulting mesh consists of 118, 556 tetrahedral cells, while the size of stiffness (S), mass (M) and ABCs (R) matrices is 141, 440 x 141, 440 each, and the dofs (edges) of interest that lie upon the surface of the sphere are only 2, 709. Therefore, the size of the impedance (admittance) matrix Z (Y) comprising the CM eigenproblem is reduced to 2, 709 x 2, 709, which is greater than 750 edges (MoM matrix 750 x 750) of the corresponding SIE formulation, reported in [14]. However, this size is the one third of the corresponding matrix developed in the VIE-MoM method [10]. Both *J*− and *M*− formulations are utilized for the analysis of the dielectric sphere and the results are compared to related published data, whenever available.

Fig. 6 displays the magnitude of the first four eigenvalues with the smallest eigenvalue in the frequency range between 4 GHz to 8 GHz computed both with the *J*-formulation (red dots) and the *M*-formulation (blue dots) presented in Section II. Comparing with the corresponding published results (blue solid lines) computed by the corresponding Hu-Wang *J*-formulation (formulation 2 of [41]), a very good agreement is observed. A serious deviation occurs around 6.6 GHz, where the Hu-Wang formulation presents an external resonance, while the proposed formulation results at an internal resonance.

However, if we compare the same eigenvalues with the results obtained from a SIE-based TCM formulation proposed by Ylä-Oijala *et al.* [14], there is no external resonance around 6.6 GHz. This indicates that the proposed algorithm identifies correctly both the external and internal resonances of the studied geometries. Of course the existence of internal resonances is not favourable for a characteristic modes formulation, and as discussed in Section II they are considered as spurious modes. At this moment
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 7. Modal significance distributions of modes shown in Fig. 6. Red and blue dots refer to the proposed *J*− and *M*-formulations, respectively, while solid blue lines to Hu-Wang *J*-formulation [41] and black dashed lines to Ylä-Oijala et al. [14].</figcaption></figure>

of our research, however, we are not focused on eliminat-
ing these resonances, even though we understand that they
appear due to the null space of the impedance/admittance
operator [21], [31]. The elimination of the internal reso-
nances is left as a future task. Nonetheless, the efficient
distinction between external and internal resonances is of
primal importance, and is achieved through the modal sig-
nificance distributions shown in Fig. 7. There, all modes
around 6.6 GHz have MS far below the unit value, which
reveals that there is no mode that efficiently radiates at this
frequency.

Taking a closer look at the eigenvalue distribution of
Fig. 6, the existence of three almost perfectly overlapping
eigenvalues at the same frequency is revealed. Note here that
for a symmetric structure like the sphere, degenerate TE<sub>lmn</sub>
and TM<sub>lmn</sub> modes are expected [42], for m > 0. Specifically
for m ≥ n, 2n + 1 modes exist with the same resonant
frequency (the corresponding lines in all the figures overlap
to each other), but with different field distributions. This is
clearly revealed when we investigate for example the current
distributions close to the first resonant frequency of Fig. 6,
where three degenerate modes at f<sub>r</sub> = 5 GHz appear. These
modes, right before the resonant TE<sub>1</sub> mode show the same
eigenvalue (λ = −3.5), but different current distributions as
shown in Fig. 8.

The results of Fig. 6 are obtained when the fictitious sur-
face is placed at a distance of λ<sub>max</sub>/2. The same λ<sub>max</sub>/2
distance is preserved for all the forthcoming numerical stud-
ies. The mesh density on the fictitious surface plays also a
key role on the accuracy of the eigenvalues' distributions, as
discussed and shown previously with the PEC plate example
(Section III-A). Herein, the descritization of the fictitious sur-
face is chosen at λ<sub>min</sub>/10. A further increase of this mesh
density will diminish any deviations away from the reso-
nances, at the expense, though, of increased computational
time. However, such a cost is unnecessary for an antenna
design procedure, since the modes around resonance are of
primal importance [40], as already mentioned.

To investigate the effect of the discretization of the
geometry under study on the accuracy of the eigenvalue
distributions, the following study is performed; the loss-
less non-magnetic dielectric sphere with α = 4.6mm and
ε<sub>r</sub> = 38 is modeled with a surface discretization of λ<sub>min</sub>/32,
and λ<sub>min</sub>/64, respectively. Fig. 9 presents the distributions
of the first six modes for the J-formulation, showing minor
deviations and at only the higher order modes.

Apart from the finer mesh in the geometry's surface, the densification of the fictitious surface mesh has also an important impact on the resulting eigenvalues at resonance. Fig. 10 depicts the relative error of the resulting eigenvalue for the TE₁ mode in resonance at 5.2 GHz, as the fictitious surface mesh is densified. The corresponding value of Hu-Wang formulation [41], is taken as the reference value. As the mesh density increases from λ_min/5 to λ_min/10 it's obvious how the relative error decreases. A further rise of the mesh step does not offer an essential improvement of the results. Thus, the λ_min/10 discretization step is chosen in all our simulations as it offers a small relative error, and at the same time an affordable mesh size that does not consume excessive computational resources.

Similar to the *J*-formulation, the *M*-formulation is also
capable of identifying both the external and internal reso-
nances of the structure under study. However, it turns out
that the *M*-formulation is more sensitive compared to the
*J*-formulation in respect to the fictitious surface's discretiza-
tion. The sensitivity of the *M*-formulation compared to the
*J*-formulation is attributed to the bad condition number of the
corresponding **L**<sub>rr</sub> matrix in (11) which as a matter of fact
is 10³ times worse than the condition number of the corre-
sponding **L**<sub>rr</sub> matrix of (17). Namely, for a discretization step
of λ<sub>min</sub>/10 a deviation of even 10² can be observed between
the mode distributions of the *J*− and *M*-formulations, mainly
at areas far from the external resonances. These results are
not shown here for reasons of brevity. To eliminate these
differences a discretization step of λ<sub>min</sub>/20 and λ<sub>min</sub>/40 is
chosen for the *J*- and *M*-formulation, respectively.

The spherical structure studied so far is lossless as it
is made of an ideal material useful only for simulation
processes. However, real structures are not loss-free. The
suggested formulation is also able to handle lossy structures
via the incorporation of complex dielectric permittivity, the
imaginary part of which is related to the dielectric loss
tangent. Considering the same geometrical non-magnetic
sphere, but with ε<sub>r</sub> = 38 – j1.0, the resulting magnitude
of the first four eigenvalues is depicted in Fig. 11, using
both the *J*− and *M*−formulations and compared against the
results published in [14].

Examining the imaginary part of the aforementioned
eigenvalues in Fig. 12 the abrupt rise of the imaginary part
is clearly evident at the internal resonances. This increase
confirms the strong absorption of the internal resonances,
although the material losses are relatively small. Otherwise,
the imaginary part of the eigenvalues is relatively small com-
pared with the real part [43]. Thus, the classification of the
<figure><img src="image_13.png" /><figcaption>FIGURE 8. Lossless sphere eigenvectors: Real part of the eigenvectors at $f_r = 5$ GHz for the first three modes of Fig. 6 at the y-z view.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 9. Effect of mesh densification: Magnitude of the first six eigenvalues with the smallest magnitude for a lossless non-magnetic dielectric sphere with $\alpha = 4.6$ mm and $\epsilon_r = 38$ computed with the proposed J-formulation methodology when the discretization step for the surface of the sphere alters from $\lambda_{min}/32$ to $\lambda_{min}/64$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 11. Lossy sphere J- (red dots) and M- (blue dots) formulations against Ylä-Oijala *et al.* [14]: Magnitudes of the first six eigenvalues with the smallest magnitude for a lossy non-magnetic dielectric sphere with $\alpha = 4.6$ mm, and $\epsilon_r = 38 - j1.0$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 10. Relative error for the $TE_1$ resonant eigenfrequency of the lossless dielectric sphere, in respect to the corresponding value of Hu-Wang formulation [41].</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIGURE 12. Lossy sphere J- (red dots) and M- (blue dots) formulations: Imaginary parts of the first four eigenvalues with the smallest magnitude for a lossy non-magnetic dielectric sphere with $\alpha = 4.6$ mm, and $\epsilon_r = 38 - j1.0$.</figcaption></figure>

resonances between external and internal ones can be also realized by observing the imaginary parts of the eigenvalues in the case of lossy structures.

To further validate our algorithm, the real part of the first three eigenvectors of the lossy dielectric sphere for modes $TE_1$ ($f = 5.2$ GHz), $TM_1$ ($f = 7.3$ GHz), and $TE_2$ ($f = 7.45$ GHz) at their resonances is shown in Fig. 13. The

$TE$-mode eigenvectors exhibit a rotation-like distribution, while the $TM$-mode has a dipole-like behavior. These are the expected electric field patterns of spherical resonators, as also presented in bibliography [44], [45].

All the simulations throughout this article were completed in an HP workstation Z820 with Intel® Xeon CPU E5-2650 v2 @ 2.60 GHz ×17 and 251.8 GHz of DDR3 memory.
<figure><img src="image_18.png" /><figcaption>FIGURE 13. Lossy sphere eigenvectors: Real part of resonance eigenvectors for the modes of Fig. 11 at the y-z view. These are in agreement with related results from [44], [45].</figcaption></figure>

Table 1 summarizes the time and memory demands of the problem for different discretization steps of both the lossless and lossy spheres. In general, the Schur complement procedure consumes most of the memory and consequently most of the time during the solution of the system. The significant increase of computational time in the lossy case is due to the additional simulation for the determination of radiation matrix. If this matrix is known a priori, then the computational time can be reduced to half.

### C. DIELECTRIC CUBE

The proposed methodology is further tested simulating a simple non-magnetic ($\mu_r = 1$) dielectric cube with side length $l = 25.4$ mm and dielectric permittivity $\epsilon_r = 9.4$ (Cube NL in [10]). This example has been extensively studied in [10] using a VIE MoM-based formulation, and the corresponding results are used herein as reference. The structure is discretized following the discretization step criteria $\lambda_{min}/32$ for the cube, $\lambda_{min}/20$ for the fictitious surface placed $\lambda_{max}/2$ away from the structure, resulting in an impedance (admittance) matrix of the order of 2, 601 × 2, 601. Fig. 14 depicts the magnitude of the first four eigenvalues' distributions in the frequency range of 2 GHz to 4.2 GHz. The red lines represent the results of the proposed $J$-formulation, the blue lines those of the $M$-formulation, while the black dashed lines stand for the results of the VIE MoM-based formulation as published in [10]. As shown, good agreement is observed between the proposed formulation (both $J$ and $M$) and the reference results.

The total computational time for this example is 3.5 hours for each frequency point, when our formulation is used. To solve the same problem using the VIE-based formulation a slightly longer time of 3.96 hours per frequency point is needed, as reported in [10].<fnref n="3" /> Therefore, the computation time of a VIE-based formulation is comparable to the computation time of the proposed FEM-based formulation, at least at its current state. Of course in both formulations

<table><caption>TABLE 1. Time and memory consumption for the dielectric sphere solution when the fictitious surface is placed at a distance $\lambda_{max}/2$ for different mesh densities.</caption><thead><tr><th colspan="8">Lossless Dielectric Sphere</th></tr><tr><th><sup>1</sup>S<sub>f<sub>i</sub></sub><br>Discr.</th><th><sup>2</sup>FEM<br>Size</th><th><sup>3</sup>S<sub>i</sub><br>Discr.</th><th><sup>4</sup>Z(Y)<br>Size</th><th><sup>5</sup>M.A.<br>Time<br>(s)</th><th><sup>6</sup>Sol.<br>Time<br>(s)</th><th><sup>7</sup>RAM<br>(GB)</th></tr></thead><tbody><tr><td>λ<sub>min</sub>/10</td><td>41, 172<br>×</td><td>λ<sub>min</sub>/32</td><td>648<br>×</td><td>0.9</td><td>271</td><td>1.73</td></tr><tr><td>λ<sub>min</sub>/10</td><td>141, 440<br>×</td><td>λ<sub>min</sub>/64</td><td>2, 709<br>×</td><td>3.5</td><td>5, 029</td><td>8</td></tr><tr><td>λ<sub>min</sub>/20</td><td>257, 471<br>×</td><td>λ<sub>min</sub>/32</td><td>648<br>×</td><td>6.8</td><td>4, 508</td><td>22</td></tr></tbody><thead><tr><th colspan="8">Lossy Dielectric Sphere</th></tr><tr><th>λ<sub>min</sub>/10</th><th>141, 440<br>×</th><th>λ<sub>min</sub>/64</th><th>2, 709<br>×</th><th>3.7</th><th>9, 699</th><th>8</th></tr><tr><th>λ<sub>min</sub>/20</th><th>257, 471<br>×</th><th>λ<sub>min</sub>/32</th><th>648<br>×</th><th>6.9</th><th>8, 336</th><th>22</th></tr></thead></table>

appropriate techniques can be utilized (e.g., MLFMA, [46], for the MoM-based and adaptive cross approximation, [47], for the FEM-based formulation) that can significantly reduce the corresponding times. The key advantage of the FEM-based formulation compared to the VIE-based formulation

<footnote n="3">3. Explicitly, [10] reports 103 hours for 26 frequency points.</footnote>

<footnote n="1"><sup>1</sup>S<sub>f<sub>i</sub></sub> Discr.: corresponds to the fictitious' surface discretization used to apply the absorbing boundary conditions.</footnote>

<footnote n="2"><sup>2</sup>FEM Size: corresponds to the FEM's matrices size [S, M, and R are as are observed in (6), (7), (8)].</footnote>

<footnote n="3"><sup>3</sup>S<sub>i</sub> Discr.: corresponds to the sphere's surface discretization.</footnote>

<footnote n="4"><sup>4</sup>Z(Y) Size: corresponds to the CM eigenproblem matrix's size that we solve to obtain the corresponding characteristic solutions.</footnote>

<footnote n="5"><sup>5</sup>M.A. Time (s): corresponds to the matrices' [see (6), (7), (8)] assembly time.</footnote>

<footnote n="6"><sup>6</sup>Sol. Time (s): corresponds to the time needed to perform the Schur complement [see (11), (17)] and solve the corresponding eigenvalue problem [see (13b), (19b)].</footnote>

<footnote n="7"><sup>7</sup>RAM (GB): corresponds to the memory requirements for the corresponding system solution [see (13b), (19b)].</footnote>
<figure><img src="image_19.png" /><figcaption>**FIGURE 14.** Validation on *J*− (red lines) and *M*− (blue lines) formulations against VIE MoM-based formulation (black dashed lines) [10]: Magnitudes of the first four eigenvalues with the smallest magnitude for a lossless non-magnetic dielectric cube with *l* = 25.4 mm and *ϵ*<sub>r</sub> = 9.4.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**FIGURE 15.** Modal significance distributions for the modes shown in Fig. 14. Red and blue lines refer to the proposed *J*− and *M*− formulations, respectively, while dashed black lines to VIE results shown in [10].</figcaption></figure>

is not in the time consumption, but in the fact that it is a Green's function-free formulation. This explicitly means that it does not require the evaluation of Green's function, thereby the study of any arbitrarily shaped, multilayered geometry loaded with anisotropic and inhomogeneous materials is feasible without applying any modification on the formulation. A VIE-based formulation can also handle any complex geometry (e.g., multilayered geometry loaded with anisotropic and inhomogeneous materials) as long as, though, the corresponding Green's function is available. If it is not available, it can be definitely evaluated which requires this corresponding analytical burden.

## D. RECTANGULAR METALLIC PATCH ANTENNA

To further verify the proposed methodology, a rectangular metallic patch antenna of length $L_p = 60$ mm and width $W_p = 40$ mm printed on a grounded dielectric substrate of length $L_s = 100$ mm, width $W_s = 80$ mm, and height $H_s = 1.5$ mm as shown in Fig. 16 is studied. The dielectric constant of the substrate is chosen equal to $\epsilon_r = 2.3$. The *J*-formulation is used to determine the corresponding characteristic modes and is compared with the results given by the commercial software FEKO [48].

<figure><img src="image_21.png" /><figcaption>**FIGURE 16.** Rectangular metallic patch antenna printed on a grounded dielectric substrate with $\epsilon_r = 2.3$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**FIGURE 17.** Rectangular metallic patch antenna *J*-formulation: Eigenvalues distribution of the rectangular patch antenna of Fig. 16. Red dots represent the results of the proposed *J*-formulation, while dashed black lines are the results obtained from FEKO [48].</figcaption></figure>

Note here, that even though we don't know the type of integral operators FEKO uses to solve this CM problem, we use it here only for reasons of comparison, since it's the only available software to our knowledge capable of performing CM analysis on composite structures. Therefore, in what follows in this section, we only present our results without implying that FEKO or our formulation is less or more accurate.

To perform this analysis we choose to discretize all the metallic parts of the patch antenna with a step size of $\lambda_{min}/32$, and the fictitious surface with a corresponding step of $\lambda_{min}/20$. Fig. 17 compares the results of our *J*-formulation with FEKO's response showing a fair agreement between them. Specifically, regarding the number of the observed modes it looks that both our formulation and FEKO have the same number of modes. However, there is indeed some slight deviation between some of the results. FEKO has an almost vertical line that crosses zero at 1.6 GHz, while the proposed formulation produces a less vertical response (red-dots) that crosses zero at the same frequency. At higher frequencies the proposed formulation captures some modes which also seems to match with a portion of FEKO's response, even though FEKO doesn't give us results prior to some frequencies. This behavior is due to the well-known FEKO's failure of mode tracking [49], which results in jagged responses in eigenvalues' distribution. These responses are eliminated from Figs. 17, 18 through a post-processing procedure.

What is important to note here is the behavior of the observable eigenvalues. As we know, patch antennas are
<figure><img src="image_23.png" /><figcaption>FIGURE 18. Modal significance distributions for the modes shown in Fig. 17. Red lines refer to the proposed J-formulation, while dashed black lines to those obtained from FEKO.</figcaption></figure>

narrowband in principle. By observing Figs. 17 and 18 we find, however, eigenvalues that cross zero and show high modal significance for a wide band of frequencies. This can be misleading and can create the impression that patch antennas are wideband, which is not the case. Let us explain this using Figs. 18 and 19. In Fig. 18 there are two different types of eigenvalue responses. Eigenvalues that show steep response versus frequency and eigenvalues that show smooth response versus frequency. The eigenvalues with steep response have the expected narrowband behavior, and correspond essentially to the modes that are related to the patch-radiator itself. Namely, at 1.6 GHz and 2.3 GHz the traditional *TM*<sub>100</sub> and *TM*<sub>010</sub> modes are, respectively, obtained, that can be also analytically evaluated following classical microstrip patch antennas' analysis, e.g., Balanis [50]. The eigenvalues with smooth response versus frequency are attributed to the finite grounded substrate, they do not radiate efficiently, as they usually don't take any energy, and if excited, they radiate in the end-fire direction. In general, these finite substrate modes are not the desired radiation modes we obtain, or we are seeking for to have, when we excite a patch antenna. To prove that indeed the smooth responses correspond to finite ground plane substrate modes that do not radiate efficiently, we compare the current distributions for some of them with the current distributions of traditional patch radiating modes. Let us take for instance the frequency of operation of 2.3 GHz. According to Fig. 18 three red lines dominate at this frequency, implying the existence of three radiating modes. One of the lines has steep response while the other two are parts of lines with smooth responses that appear for a wide band of frequencies. By plotting the corresponding currents of these modes in Fig. 19, the nature of each mode is clear. Namely, the first from the left figure shows the traditional current distribution with high current intensity on the top metallic layer of our patch antenna, flowing perpendicular to its long edge. On the other hand, for both the middle and the right figures, the current intensity is significantly lower with minimal current flowing over the top metallic layer of the patch, indicating

that these modes are not able to efficiently radiate. Note
that all three current distributions are plotted using the same
scale.

What is interesting from the CM point of view is to
develop a mechanism that can identify the true patch-radiator
modes solely and distinguish them from the “undesired”
(under most of the circumstances) substrate modes. However,
this is a subject that is beyond the tasks of the present effort.

E. ANISOTROPIC MULTILAYERED SPHERICAL
RESONATOR

So far, the proposed formulation has been tested and val-
idated against typical SIE CM formulations. However, the
use of FEM allows the study of arbitrarily shaped mul-
tilayered structures filled with anisotropic materials. This
unique advantage is also offered by VIE but with a pro-
hibitive increase of computational demands [10]. Herein,
the anisotropic spherical resonator of Fig. 20 is analyzed,
utilizing the *M*-formulation. The resonator consists of two
concentric spheres with different dielectric properties. The
inner sphere has radius 0.7 mm, and is uniaxially anisotropic
with $\bar{\epsilon}_{r1} = (90, 90, 9)$. The external spherical shell has outer
radius 3.5 mm, and $\epsilon_{r2} = 10$. The inner sphere herein is dis-
cretized with $\lambda_{\min}/50$ while the outer with $\lambda_{\min}/32$ resulting
in a CM admittance matrix *Y* of size 4, 201 × 4, 201. As
expected, the study of multilayered structures increases the
number of dofs at the “radiating” surface $S_i$ of the outer
sphere. Although the discretization criteria for the outer
sphere are the same with the dielectric sphere studied earlier
($\lambda_{\min}/64$ for the surface of the outer sphere and $\lambda_{\min}/10$ for
the fictitious surface) and its radius is smaller, the resulting
dofs are significantly increased due to the necessary finer dis-
cretization of the inner sphere which has significant impact
on the whole mesh. This inconvenience appears mainly in
spherical topologies, where a denser mesh is needed in order
to represent the curved interfaces more accurately. However,
this rise of CM admittance matrix *Y* is not excessive to
render the solution of the system prohibitive. Notably, the
uniaxial anisotropy of the inner sphere does not cause any
increase neither in the dofs nor in the computational cost. The
distribution of the first eight eigenvalues with the smallest
magnitude is presented in Fig. 20, where two modes res-
onate, the mode *TM*<sub>110</sub> (orange line) at *f*<sub>r</sub> = 11.41 GHz,
and the *TE*<sub>110</sub> (blue line) at *f*<sub>r</sub> = 12.97 GHz. The resonance
of *TE*<sub>110</sub> complies with Wolff’s analytical solution (vertical
black dashed line) [51]. The electric field distributions of two
modes at resonance are depicted in Fig. 21 and clearly reveal
the rotation-like behavior of TE modes [Fig. 21(a)], and the
corresponding dipole-like behavior of TM ones [Fig. 21(b)],
in accordance with the electric field patterns of spherical
resonators presented in bibliography [44], [45].

As it is shown throughout this work, the proposed Green's function-free methodology for CM is able to easily handle multilayered geometries loaded with anisotropic materials without any further modification in the formulation. An initial attempt to solve multilayered structures with the use
<figure><img src="image_24.png" /><figcaption>FIGURE 19. Current distribution for the three red lines of Fig. 18 at 2.3 GHz for the patch antenna of Fig. 16.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 20. Anisotropic multilayered sphere M-formulation: First eight eigenvalues of an anisotropic spherical resonator consisting of two concentric non-magnetic dielectric spheres. The radius and dielectric permittivity of the inner sphere are $r_1 = 0.7$ mm and $\bar{\epsilon}_{r1} = (90,90,9)$, respectively, while the radius and dielectric permittivity of the outer sphere are $r_2 = 3.5$ mm and $\epsilon_{r2} = 10$, respectively. The vertical black dashed line represents the analytical resonant frequency given in [51].</figcaption></figure>

of the conventional SIE-based CM analysis is made in the latest work of Ylä-Oijala [7]. Therein, the initial formulation is modified in order to express the equivalent electric and magnetic currents at all interfaces between adjacent layers. On the contrary, the herein proposed FEM-CM formulation can be readily applied without the need of modification.

## F. CIRCULAR PATCH ANTENNA WITH FERRITE POST

The last example describes a fully inhomogeneous geometry filled with an anisotropic material. This is a circular patch antenna loaded with a ferrite post offering frequency tunability as shown in Fig. 22. To make the antenna able to resonate at the ISM band/ tunable around 2.42 GHz, for its dominant $TM_{110}$ mode, the radius of the circular patch is chosen as $b = 7.5$ mm. This geometry was proposed and extensively studied in the past by our research group [52]. Following [52], the radius of the ferrite post, $a$, is chosen to be half of the radius of the circular patch, thus $a = 0.5b$. Ceramic-PTFE $TMM10$ with $\epsilon_{rD} = 9.2$, $\tan\delta = 0.0009$, and thickness $h = 0.8$ mm is used for our substrate. The

<figure><img src="image_26.png" /><figcaption>FIGURE 21. Anisotropic multilayered sphere: Real part of resonant eigenvectors for: a) $TM_{110}$ at $f_r = 11.41$ GHz, and b) $TE_{110}$ $f_r = 12.97$ GHz modes of Fig. 20 at the y-z view.</figcaption></figure>

ferrite is a YIG-Al doped type GA-65 ferrite Domen with $4\pi M_s = 650$ Gauss [$\equiv 10^{-4}$ T], $\epsilon_{rF} = 14.2$ and $\Delta H = 45$ Oe [$\equiv (4\pi)^{-1} \times 10^3$ A/m].
<figure><img src="image_27.png" /><figcaption>FIGURE 22. Geometry of a circular patch antenna tuned by transversely magnetized ferrite post.</figcaption></figure>

The anisotropic nature of this geometry stems from the transversely magnetized-biased ferrite post, whose magnetic permeability is of the form:

$$ \bar{\bar{\mu}} = \begin{bmatrix} \mu_F & j\kappa_F & 0 \\ -j\kappa_F & \mu_F & 0 \\ 0 & 0 & \mu_0 \end{bmatrix} \qquad (27) $$

where:

$$ \mu_F = \mu_0 \left( 1 + \frac{\omega_m (\omega_0 + j\alpha_F \omega)}{(\omega_0 + j\alpha_F \omega)^2 - \omega^2} \right) = \mu_0 \mu_{rF} \quad (28) $$

and

$$ \kappa_F = \mu_0 \frac{\omega \omega_m}{(\omega_0 + j\alpha_F \omega)^2 - \omega^2} = \mu_0 \kappa_{rF} \quad (29) $$

where $\omega = 2\pi f$ refers to the operating angular frequency, $\alpha_F = \frac{\mu_0\gamma\Delta H}{2\omega_0}$ is the damping factor, $\gamma$ is the gyromagnetic ratio, $\Delta H$ is the ferrite linewidth [53], $\omega_0 = \mu_0\gamma H_0 = 2\pi f_0$, $\omega_m = \mu_0\gamma M_s = 2\pi f_m$, $H_0$ is the DC biasing magnetic field, and $M_s$ is the saturation magnetization. The subscript F is used to denote the ferrite. Since the scope, herein, is to explore the capabilities of FEM-CM formulation, tunability is beyond our purpose. Thus, some fixed magnetic biases are assumed. Indicatively, $H_0 = 63661.98$ A/m results in $f_0 = 2.8$ MHz/A/m $\cdot H_0(A/m) = 2.24$ GHz, $f_m = 2.8$ MHz/A/m $\cdot 4\pi M_s(T) = 1.82$ GHz, and $a_F = 2.8$ MHz/A/m $\cdot \Delta H(A/m)/2\omega_0 = 0.0045$ according to [53]. These values are then substituted in (28) and (29) to yield $\mu_{rF} = -3.842 - j0.3036$ and $\kappa_{rF} = -5.2326 - j0.3027$.

Fig. 23 presents the real part of the resulting eigenvalue distributions for different DC magnetizations, compared to the results provided by the analytical method in [52]. An almost perfect agreement is achieved for the resonant frequencies between the proposed algorithm and the published results. The slight deviation of 0.6% around the frequency of $f_r = 3.59$ GHz can be due to the approximations of [52] and/or the mesh density herein.

Fig. 24 presents a 1-D representation of the absolute value for the normalized $\phi$ component-magnetic current distribution around the edge of the circular patch antenna

<figure><img src="image_28.png" /><figcaption>FIGURE 23. Real part of the eigenvalues for a circular patch antenna loaded with a ferrite post. Colourful lines show the proposed method results, while the black symbols stand for the resonant frequencies shown in [52].</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>FIGURE 24. Absolute value of the normalized $\phi$ component-magnetic current distribution around the edge of the circular patch antenna for the TM110 mode. The blue solid line indicates the result from the proposed FEM-CM formulation, while the dotted red line the results from [54].</figcaption></figure>

for the TM110 mode. This is the first resonant mode at $f_r = 2.42$ GHz, when the DC biasing magnetic field is $H_0 = 63661.98$ A/m. The resulting distribution of the proposed methodology is qualitatively compared to its analytical counterpart. Namely, this normalized $\phi$ magnetic current distribution around the edge of the circular patch antenna is in accordance with the expected theoretical analysis met in [54] for the case of TM110 mode. The observed deviations are due to the fact that the theoretical analysis is made for a lossless circular patch (here we model a circular patch antenna loaded with lossy material), and due to the relatively coarse mesh discretization in this area.

## IV. CONCLUSION

This work presented thoroughly a novel method for the study of electromagnetic structures utilizing a characteristic modes analysis based on the finite element scheme. The validity of the method was proved through the study of known structures, whose resulting eigenvalues and eigenvectors were cross-examined with the corresponding results in bibliography. Thus, a wide new era arises in the field of characteristic modes, as any arbitrary shaped multilayered electromagnetic structure loaded with possibly anisotropic and/or highly inhomogeneous materials can be studied based on the robustness
and flexibility of finite element method. It should be men-
tioned that no modification of the formulation is needed in
order to study either isotropic or anisotropic material with
and without inhomogeneity.

The method presented in this work is a computational
paradigm shift for formulating CM problems, opening sev-
eral new directions of research. Amongst them we have
identified three main problems that are seeking solution:

* Investigate the reasons of internal resonances in both proposed *J*− and *M*−formulations, and eliminate them with the appropriate modification of the corresponding formulation (e.g., introduction of an equivalent formulation to the combined field integral equation approach [15], or introduction of a DtN approach [33]).

* Investigate the reasons of why the proposed *M*-formulation is more sensitive compared to the proposed *J*-formulation.

* Investigate the reasons of why the CM analysis of a patch antenna returns eigenvalues that appear to be close to zero across a very wide frequency band, implying that this structure has wideband performance, which is against the common narrowband behavior of a patch antenna.

REFERENCES

[1] R. J. Garbacz, “A generalized expansion for radiated and scat-
tered fields,” Ph.D. dissertation, Dept. Elect. Eng., Ohio State Univ.,
Columbus, OH, USA, 1968.

[2] Y. Chang and R. Harrington, "A surface formulation for characteristic modes of material bodies," *IEEE Trans. Antennas Propag.*, vol. AP-25, no. 6, pp. 789-795, Nov. 1977.

[3] Y. Chen and C.-F. Wang, *Characteristic Modes: Theory and Applications in Antenna Engineering*. Hoboken, NJ, USA: Wiley, 2015.

[4] R. F. Harrington and J. R. Mautz, “Theory of characteristic modes for conducting bodies,” *IEEE Trans. Antennas Propag.*, vol. AP-19, no. 5, pp. 622–628, Sep. 1971.

[5] R. F. Harrington and J. R. Mautz, "Computation of characteristic modes for conducting bodies," *IEEE Trans. Antennas Propag.*, vol. AP-19, no. 5, pp. 629–639, Sep. 1971.

[6] R. Harrington, J. Mautz, and Y. Chang, "Characteristic modes for dielectric and magnetic bodies," *IEEE Trans. Antennas Propag.*, vol. AP-20, no. 2, pp. 194–198, Mar. 1972.

[7] P. Ylä-Oijala, "Generalized theory of characteristic modes," *IEEE Trans. Antennas Propag.*, vol. 67, no. 6, pp. 3915-3923, Jun. 2019.

[8] A. H. Nalbantoğlu, "New computation method for characteristic modes," *Electron. Lett.*, vol. EL-18, no. 23, pp. 994-996, Nov. 1982.

[9] R. T. Maximidis, C. L. Zekios, T. N. Kaifas, E. E. Vafiadis, and G. A. Kyriacou, “Characteristic mode analysis of composite metal-dielectric structure, based on surface integral equation/moment method,” in *Proc. 8th Eur. Conf. Antennas Propag. (EuCAP)*, 2014, pp. 2822–2826.

[10] Z. T. Miers and B. K. Lau, "Computational analysis and verifications of characteristic modes in real materials," *IEEE Trans. Antennas Propag.*, vol. 64, no. 7, pp. 2595-2607, Jul. 2016.

[11] Z. Miers and B. K. Lau, "On characteristic Eigenvalues of complex media in surface integral formulations," *IEEE Antennas Wireless Propag. Lett.*, vol. 16, pp. 1820-1823, 2017.

[12] Y. Chen and C.-F. Wang, "Surface integral equation based characteristic mode formulation for dielectric resonators," in *Proc. IEEE Antennas Propag. Soc. Int. Symp. (APSURSI)*, Memphis, TN, USA, 2014, pp. 846-847.

[13] Y. Chen, “Alternative surface integral equation-based characteristic mode analysis of dielectric resonator antennas,” *IET Microw. Antennas Propag.*, vol. 10, no. 2, pp. 193–201, 2016.

[14] P. Ylä-Oijala, H. Wallén, D. C. Tzarouchis, and A. Sihvola, “Surface integral equation-based characteristic mode formulation for penetrable bodies,” *IEEE Trans. Antennas Propag.*, vol. 66, no. 7, pp. 3532–3539, Jul. 2018.

[15] Q. I. Dai, Q. S. Liu, H. U. I. Gan, and W. C. Chew, "Combined field integral equation-based theory of characteristic mode," *IEEE Trans. Antennas Propag.*, vol. 63, no. 9, pp. 3973–3981, Sep. 2015.

[16] P. Ylä-Oijala, J. Lappalainen, and S. Järvenpää, “Characteristic mode equations for impedance surfaces,” *IEEE Trans. Antennas Propag.*, vol. 66, no. 1, pp. 487–492, Jan. 2018.

[17] R. Lian, J. Pan, and S. Huang, “Alternative surface integral equation formulations for characteristic modes of dielectric and magnetic bodies,” *IEEE Trans. Antennas Propag.*, vol. 65, no. 9, pp. 4706–4716, Sep. 2017.

[18] M. Gustafsson, L. Jelinek, K. Schab, and M. Capek, "Unified theory of characteristic modes: Part I—Fundamentals," 2021, arXiv:2109.00063.

[19] M. Gustafsson, L. Jelinek, K. Schab, and M. Capek, "Unified theory of characteristic modes: Part II—Tracking, losses, and FEM evaluation," 2021, arXiv:2110.02106.

[20] R. E. Collin, *Field Theory of Guided Waves* (IEEE Press Series on Electromagnetic Wave Theory). New York. NY. USA: McGraw-Hill, 1969.

[21] W. C. Chew, M. S. Tong, and B. Hu, *Integral Equation Methods for Electromagnetic and Elastic Waves* (Synthesis Lectures on Computational Electromagnetics). San Rafael, CA, USA: Morgan Claypool Publ., 2009.

[22] C. L. Zekios, P. C. Allilomes, and G. A. Kyriacou, “DC and imaginary spurious modes suppression for both unbounded and lossy structures,” *IEEE Trans. Microw. Theory Techn.*, vol. 63, no. 7, pp. 2082–2093, Jul. 2015.

[23] F. J. C. Meyer and D. B. Davidson, “Adaptive-mesh refinement of finite-element solutions for two-dimensional electromagnetic problems,” *IEEE Antennas Propag. Mag.*, vol. 38, no. 5, pp. 77–83, Oct. 1996, doi: 10.1109/74.544405.

[24] F. J. C. Meyer and D. B. Davidson, "Error estimates and adaptive procedures for the two-dimensional finite element method," *Electron. Lett.*, vol. 30, no. 12, pp. 936-938, 1994.

[25] F. Gundes, K. Zhao, and J.-F. Lee, "Implementing higher order absorbing boundary conditions in vector finite element methods," in *Proc. IEEE Antennas Propag. Soc. Int. Symp.*, Albuquerque, NM, USA, 2006, pp. 1787–1790.

[26] C. L. Zekios, D. G. Makris, R. T. Maximidis, P. C. Allilomes, and G. A. Kyriacou, “A characteristic mode analysis of parallel resonances exploiting the finite element scheme,” in *Proc. Eur. Microw. Conf. (EuMC)*, Paris, France, 2015, pp. 179–182.

[27] F.-G. Hu and C.-F. Wang, “FE-BI formulations for characteristic modes,” *IEEE Trans. Microw. Theory Techn.*, vol. 64, no. 5, pp. 1396–1401, May 2016.

[28] R. F. Harrington, *Time-Harmonic Electromagnetic Field*. New York, NY, USA: Wiley, 2001.

[29] Y. Zhu and A. C. Cangellaris, *Multigrid Finite Element Methods for Electromagnetic Field Modeling*. Piscataway, NJ, USA: IEEE Press, 2006.

[30] J. Jin, *The Finite Element Method in Electromagnetics*, 2nd ed. New York, NY, USA: Wiley, 1993.

[31] M. Capek, V. Losenicky, L. Jelinek, and M. Gustafsson, "Validating the characteristic modes solvers," *IEEE Trans. Antennas Propag.*, vol. 65, no. 8, pp. 4134–4145, Aug. 2017.

[32] R. F. Harrington, "Boundary integral formulations for homogeneous material bodies," *J. Electromagn. Waves Appl.*, vol. 3, no. 1, pp. 1–15, 1989.

[33] P. C. Allilomes and G. A. Kyriacou, “A nonlinear finite-element Leaky-waveguide solver,” *IEEE Trans. Microw. Theory Techn.*, vol. 55, no. 7, pp. 1496–1510, Jul. 2007.

[34] F. Zhang, *The Schur Complement and Its Applications*. New York, NY, USA: Springer, Jan. 2005.

[35] B. Enquist and A. Majda, "Radiation boundary conditions for acoustic and elastic wave calculations," *Commun. Pure Appl. Math.*, vol. 32, pp. 313–357, May 1979.
<figure><img src="image_30.png" /></figure>

**CONSTANTINOS L. ZEKIOS** (Member, IEEE) received the Diploma degree (with Hons.) in electrical and computer engineering, the M.S. degree (with Hons.) in electrical and computer engineering communication and satellite telecommunication systems, and the Ph.D. degree (with Hons.) in electrical and computer engineering from the Democritus University of Thrace, Xanthi, Greece, in 2008, 2011, and 2015, respectively. From January 2016 to May 2018, he was a Postdoctoral Researcher with the Electrical and Computer Engineering Department, University of Massachusetts at Amherst, Amherst, MA, USA, and from May 2018 to August 2020, he was a Postdoctoral Research Fellow with the Transforming Antennas Center, Florida International University, Miami, FL, USA, where he is currently a Research Assistant Professor with the Department of Electrical and Computer Engineering. His main research interests include computational electromagnetics, optimization methods, antennas, antenna arrays, electromagnetic surfaces, beamforming networks, and microwave engineering.

<figure><img src="image_31.png" /></figure>

**STAVROS V. GEORGAKOPOULOS** (Senior Member, IEEE) received the Diploma degree in electrical engineering from the University of Patras, Patras, Greece, in June 1996, and the M.S. and Ph.D. degrees in electrical engineering from Arizona State University, Tempe, AZ, USA, in 1998 and 2001, respectively. From 2001 to 2007, he was the Principal Engineer with SV Microwave, Inc. Since 2007, he has been with the Department of Electrical and Computer Engineering, Florida International University (FIU), Miami, FL, USA, where he is currently a Professor, the Director of the Transforming Antennas Center (a research center on foldable/origami, physically reconfigurable, and deployable antennas), and the Director of the RF Communications, Millimeter-Waves, and Terahertz Lab. His current research interests relate to novel antennas, arrays, RFID, microwave and RF systems, novel sensors, and wireless powering of portable, wearable, and implantable devices. He received the 2015 FIU President's Council Worlds Ahead Faculty Award, which is the highest honor FIU extends to a Faculty Member for excelling in research, teaching, mentorship, and service. He served as an Associate Editor for the IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION from 2013 to 2019, and has been serving as an Associate Editor for the IEEE OPEN JOURNAL OF ANTENNAS AND PROPAGATION since 2019.

<figure><img src="image_32.png" /></figure>

**GEORGE A. KYRIACOU** (Senior Member, IEEE) was born in Famagusta, Cyprus, in March 25, 1959. He received the Electrical Engineering Diploma and Ph.D. degrees (with Hons.) from the Democritus University of Thrace, Xanthi, Greece, in 1984 and 1988, respectively. Since January 1990, he has been with the Department of Electrical and Computer Engineering, Democritus University of Thrace, where he is currently a Professor and the Director of the Microwaves Laboratory and was the Director of the Graduate Studies from 2005 to 2010. He has authored over 250 journal and conference papers and supervised eight Ph.D. and 15 M.Sc. theses and more than 110 Diploma theses. His main research interests include microwave engineering, open waveguides and antennas in anisotropic media, software defined and cognitive radio, computational electromagnetics, and biomedical engineering. He is a member of the Technical Chamber of Greece and the European Microwave Association. He has been serving as an Associate Editor for IET MICROWAVES ANTENNAS & PROPAGATION since 2015.

[36] P. Monk, *Finite Element Methods for Maxwell's Equations*. Oxford, U.K.: Oxford Scholarship Online, Sep. 2007, doi: 10.1093/acprof:oso/9780198508885.001.0001.

[37] J. Lappalainen, P. Ylä-Oijala, D. C. Tzarouchis, and A. Sihvola, "Resonances of characteristic modes for perfectly conducting objects," *IEEE Trans. Antennas Propag.*, vol. 65, no. 10, pp. 5332–5339, Oct. 2017.

[38] P. Ylä-Oijala, A. Lehtovuori, H. Wallén, and V. Viikari, "Coupling of characteristic modes on PEC and lossy dielectric structures," *IEEE Trans. Antennas Propag.*, vol. 67, no. 4, pp. 2565–2573, Apr. 2019.

[39] T. A. Davis, *Direct Methods for Sparse Linear Systems*. Philadelphia, PA, USA: Soc. Ind. Appl. Math., 2006.

[40] M. Cabedo-Fabres, E. Antonino-Daviu, A. Valero-Nogueira, and M. F. Bataller, "The theory of characteristic modes revisited: A contribution to the design of antennas for modern applications," *IEEE Antennas Propag. Mag.*, vol. 49, no. 5, pp. 52–68, Oct. 2007.

[41] F.-G. Hu and C.-F. Wang, "Integral equation formulations for characteristic modes of dielectric and magnetic bodies," *IEEE Trans. Antennas Propag.*, vol. 64, no. 11, pp. 4770–4776, Nov. 2016.

[42] M. R. Currie, "The utilization of degenerate modes in a spherical cavity," *J. Appl. Phys.*, vol. 24, no. 8, pp. 998–1003, 1953.

[43] H. Wallén, P. Ylä-Oijala, D. C. Tzarouchis, and A. Sihvola, "Mie scattering and characteristic modes of lossy dielectric objects," in *Proc. 2nd URSI Atlantic Radio Sci. Meeting (AT-RASC)*, Gran Canaria, Spain, 2018, pp. 1–4.

[44] C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles*. New York, NY, USA: Wiley, 1983.

[45] M. Pascale, G. Miano, R. Tricarico, and C. Forestiere, "Full-wave electromagnetic modes and hybridization in nanoparticle dimers," *Sci. Rep.*, vol. 9, no. 1, Oct. 2019, Art. no. 14524.

[46] C. Lu and W. C. Chew, "A multilevel algorithm for solving a boundary integral equation of wave scattering," *Microw. Opt. Technol. Lett.*, vol. 7, pp. 466–470, Jul. 1994.

[47] K. Zhao, M. N. Vouvakis, and J.-F. Lee, "The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems," *IEEE Trans. Electromagn. Compat.*, vol. 47, no. 4, pp. 763–773, Nov. 2005.

[48] "Altair® Feko®," 2020. [Online]. Available: https://www.altair.com/feko/

[49] B. D. Raines and R. G. Rojas, "Wideband characteristic mode tracking," *IEEE Trans. Antennas Propag.*, vol. 60, no. 7, pp. 3537–3541, Jul. 2012.

[50] C. A. Balanis, *Antenna Theory: Analysis and Design*. Chichester, U.K.: Wiley, 2005.

[51] I. Wolff, "A generalized description of the spherical three-layer resonator with an anisotropic dielectric material," in *Proc. IEEE MTT-S Int. Microw. Symp. Dig.*, vol. 1. Palo Alto, CA, USA, 1987, pp. 307–310.

[52] A. Mavridis, G. Kyriacou, and J. Sahalos, "On the design of patch antennas tuned by transversely magnetized lossy ferrite including a novel resonating mode," *Progr. Electromagn. Res.*, vol. 62, pp. 165–192, May 2006.

[53] D. M. Pozar, *Microwave Engineering*. Hoboken, NJ, USA: Wiley, 2012.

[54] R. Garg, P. Bhartia, I. Bahl, and A. Ittipiboon, *Microstrip Antenna Design Handbook*. Boston, MA, USA: Artech House Inc., 2001.

<figure><img src="image_33.png" /></figure>

**KONSTANTINOS D. PASCHALOUDIS** (Member, IEEE) was born in Pylaia, Thessaloniki, Greece, in January 1991. He received the Diploma degree in electrical and computer engineering and the Ph.D. degree (with Hons.) in electrical and computer engineering from the Democritus University of Thrace, Xanthi, Greece, in 2015 and 2021, respectively. Since January 2022, he has been a Postdoctoral Fellow with the Institut d'Electronique et de Télécommunications de Rennes, Université de Rennes 1, France. His

main research interests include computational electromagnetics, microwave engineering, and analysis and design of microwave/mmWave antennas.