Validating the Characteristic Modes Solvers

TL;DR

This paper introduces analytical formulas for characteristic modes of spherical shells, providing benchmarks to validate computational methods and analyze their performance across various parameters.

Contribution

It presents closed-form formulas for characteristic modes, enabling precise validation of numerical solvers and assessment of mode tracking algorithms.

Findings

Commercial and academic solvers reliably identify dominant modes.

Higher-order modes are difficult to accurately identify with current methods.

Mode tracking routines perform poorly despite recent advancements.

Abstract

Characteristic modes of a spherical shell are found analytically as spherical harmonics normalized to radiate unitary power and to fulfill specific boundary conditions. The presented closed-form formulas lead to a proposal of precise synthetic benchmarks which can be utilized to validate the method of moments matrix or performance of characteristic mode decomposition. Dependence on the mesh size, electrical size and other parameters can systematically be studied, including the performance of various mode tracking algorithms. A noticeable advantage is the independence on feeding models. Both theoretical and numerical aspects of characteristic mode decomposition are discussed and illustrated by examples. The performance of state-of-the-art commercial simulators and academic packages having been investigated, we can conclude that all contemporary implementations are capable of identifying the first dominant modes while having severe difficulties with higher-order modes. Surprisingly poor performance of the tracking routines is observed notwithstanding the recent ambitious development.

Tables

Table 1. TABLE I: Number of theoretically achievable TE and TM modes dependent on used numerical precision in Matlab and electrical size  k ​ a 𝑘 𝑎 ka . The number of TE and TM modes N 𝑁 N was calculated according to ( 21 ) with  p max TE superscript subscript 𝑝 max TE p_{\mathrm{max}}^{\mathrm{TE}} and  p max TM superscript subscript 𝑝 max TM p_{\mathrm{max}}^{\mathrm{TM}} substituted.

$ka$		Single p.	Double p.	Quadruple p.
	bits $b$	$23$	$52$	$112$
	${\log }_{10}\left(2^b\right)$	$\approx 6.924$	$\approx 15.65$	$\approx 33.72$
$0.5$	$p_{\max }^{\mathrm{TE}/\mathrm{TM}}$	$4/4$	$7/7$	$12/12$
	$N\left(p_{\max }^{\mathrm{TE}/\mathrm{TM}}\right)$	$24/24$	$63/63$	$168/168$
$1.5$	$p_{\max }^{\mathrm{TE}/\mathrm{TM}}$	$6/6$	$10/10$	$17/17$
	$N\left(p_{\max }^{\mathrm{TE}/\mathrm{TM}}\right)$	$48/48$	$120/120$	$323/323$

Equations

$$\mathcal{X}\left({\mathbf{J}}_n\right)={\lambda }_n\mathcal{R}\left({\mathbf{J}}_n\right)$$

$$\mathcal{Z}\left({\mathbf{J}}_n\right)=\mathcal{R}\left({\mathbf{J}}_n\right)+\mathrm{j}\mathcal{X}\left({\mathbf{J}}_n\right)=\hat{\mathbf{n}}\times \hat{\mathbf{n}}\times \mathbf{E}\left({\mathbf{J}}_n\right),$$

$$\begin{array}{rl}\frac{1}{2}\underset{\Omega }{\int }{\mathbf{J}}_m^{\ast }\cdot \mathcal{Z}\left({\mathbf{J}}_n\right)\mathrm{d}S=& \frac{1}{2}\underset{\Omega }{\int }{\mathbf{J}}_m^{\ast }\cdot \mathbf{E}\left({\mathbf{J}}_n\right)\mathrm{d}S\\ =& \left(1+\mathrm{j}{\lambda }_n\right){\delta }_{mn}={\kappa }_n{\delta }_{mn},\end{array}$$

$${\lambda }_n=\frac{\mathrm{Im}\left\{{\kappa }_n\right\}}{\mathrm{Re}\left\{{\kappa }_n\right\}}=\frac{\underset{\Omega }{\int }{\mathbf{J}}_n^{\ast }\cdot \mathcal{X}\left({\mathbf{J}}_n\right)\mathrm{d}S}{\underset{\Omega }{\int }{\mathbf{J}}_n^{\ast }\cdot \mathcal{R}\left({\mathbf{J}}_n\right)\mathrm{d}S},$$

$${\mathbf{E}}_{pq}^{\mathrm{TE}}\left(r\ge a\right)=-C_{pq}\left(ka\right){\gamma }_p\left(ka\right){\mathbf{M}}_{pq}\left({\mathrm{h}}_p^{(2)},r,\vartheta ,\phi \right),$$

$${\mathbf{E}}_{pq}^{\mathrm{TM}}\left(r\ge a\right)=C_{pq}\left(ka\right){\gamma }_p\left(ka\right){\mathbf{N}}_{pq}\left({\mathrm{h}}_p^{(2)},r,\vartheta ,\phi \right)$$

$${\mathbf{J}}_{pq}^{\mathrm{TE}}=C_{pq}\left(ka\right)\hat{\mathbf{r}}\times {\mathbf{N}}_{pq}\left({\mathrm{j}}_p,a,\vartheta ,\phi \right),$$

$${\mathbf{J}}_{pq}^{\mathrm{TM}}=C_{pq}\left(ka\right)\hat{\mathbf{r}}\times {\mathbf{M}}_{pq}\left({\mathrm{j}}_p,a,\vartheta ,\phi \right),$$

$${\gamma }_p\left(ka\right)=Z_{0}ka{\mathrm{j}}_p\left(ka\right)\frac{\partial \left(ka{\mathrm{j}}_p\left(ka\right)\right)}{\partial ka},$$

$$C_{pq}\left(ka\right)=\frac{k}{{\gamma }_p\left(ka\right)}\sqrt{\frac{Z_0\left(2p+1\right)\left(p-q\right)!}{\pi \left(1+{\delta }_{q0}\right)p\left(p+1\right)\left(p+q\right)!}},$$

$${\lambda }_p^{\mathrm{TE}}=-\frac{{\mathrm{y}}_p\left(ka\right)}{{\mathrm{j}}_p\left(ka\right)}$$

$${\lambda }_p^{\mathrm{TM}}=-\frac{\left(p+1\right){\mathrm{y}}_p\left(ka\right)-ka{\mathrm{y}}_{p+1}\left(ka\right)}{\left(p+1\right){\mathrm{j}}_p\left(ka\right)-ka{\mathrm{j}}_{p+1}\left(ka\right)},$$

$${\mathbf{J}}_n=\left\{{\mathbf{J}}_{pq}^{\mathrm{TE}},{\mathbf{J}}_{pq}^{\mathrm{TM}}\right\},$$

$${\delta }_n=180\left(1-\frac{1}{\pi }\arctan \left({\lambda }_n\right)\right).$$

$$\mathbf{Z}=\mathbf{R}+\mathrm{j}\mathbf{X}=\left[Z_{uv}\right]=\left[\underset{\Omega }{\int }{\mathbf{\psi }}_u^{\ast }\mathcal{Z}\left({\mathbf{\psi }}_v\right)\mathrm{d}S\right]$$

$$\mathbf{X}{\mathbf{I}}_n={\lambda }_n\mathbf{R}{\mathbf{I}}_n$$

$${\mathbf{J}}_n=\sum _u{I}_{nu}{\mathbf{\psi }}_u.$$

$${\mathbf{Z}}_{\mathrm{sym}}=\frac{1}{2}\left(\mathbf{Z}+{\mathbf{Z}}^{\mathrm{T}}\right),$$

$$\mathbf{R}{\hat{\mathbf{I}}}_n={\xi }_n{\hat{\mathbf{I}}}_n$$

$${\mathbf{R}}_{\mathrm{pos}}=\hat{\mathbf{I}}\mathbf{\xi }{\hat{\mathbf{I}}}^{\mathrm{T}},$$

$$N\left(p\right)=2p+1,$$

$${\chi }_n=\underset{p}{\max }\sqrt{\sum _{q=-p}^p{\left|\underset{\Omega }{\int }{\hat{\mathbf{J}}}_n\left(\mathbf{r}\right)\cdot {\hat{\mathbf{J}}}_{pq}^{TM/TE}\left(\mathbf{r}\right)\mathrm{d}S\right|}^2},$$

$$\hat{\mathbf{J}}=\frac{\mathbf{J}}{\sqrt{\underset{\Omega }{\int }\mathbf{J}\left(\mathbf{r}\right)\cdot \mathbf{J}\left(\mathbf{r}\right)\mathrm{d}S}}.$$

$${\zeta }_n=\frac{\underset{p}{\max }\sum _{q=-p}^p{\left|P_{pq,n}^{TE/TM}\right|}^2}{\sum _{p,q}\sum _{TE/TM}{\left|P_{pq,n}^{TE/TM}\right|}^2}$$

$$P_{pq,n}^{TE/TM}=\frac{1}{2Z_0}\underset{4\pi }{\int }{\left({\mathbf{F}}_{pq}^{TE/TM}\right)}^{\ast }\cdot {\mathbf{F}}_n\sin \vartheta \mathrm{d}\vartheta \mathrm{d}\phi .$$

Full text

I Introduction

Characteristic mode (CM) decomposition [1, 2] has become a popular tool for analyzing and designing scatterers and antennas, mainly due to the physical insight gained by modal decomposition without a particular feeding considered [3, 4]. CM decomposition yields a set of real-valued currents which form an orthonormal basis with respect to their radiation patterns and the useful properties of the CMs render this technique appealing for antenna designers [5, 6, 7]. Consequently, CM decomposition has been the subject of implementation into commercial tools, such as FEKO [8], WIPL-D [9], and CST [10], and there also exist a plethora of academic tools, employed primarily for research related to the CM [11, 12, 13, 14].

The amount of scientific data generated along with the publication activity in the field of CMs is immense. It is therefore surprising that the question of how accurate these results are is scarcely assessed. The rare exception is an early study by Mautz and Harrington, where the results of their FORTRAN implementation [15] of CM decomposition is compared to the first analytically known eigenvalues of a spherical shell [16, Table II]. This lack of detailed numerical benchmarking of the CMs was the main motivation for the developments presented in this paper.

Benchmarking activities, see [17, 18, 19, 20] and the references therein for examples, are essential for the validation and quality assessment of methods and tools of computational electromagnetics and are of particular interest for those methods with known numerical issues, as is the case of CM decomposition [21, 22, 23, 24]. As suggested in [17], four benchmarks are readily available in computational electromagnetics, namely the comparisons to a closed form solution [25, 19], to a standard problem [26, 27], to a measurement or to other modeling techniques [28, 29, 30]. The advantages and disadvantages of these possibilities are detailed in [17] and have led authors to the decision to select a comparison to the analytical model. The drawback of this choice is that analytic solutions are only available for canonical geometries such as ellipsoids and cylinders [25, 19]. The advantage of negligible error levels in the analytic model [31, 32], however, outweighs it.

In this paper we propose four independent benchmarks devoted to various aspects of CM decomposition to validate characteristic eigenvalues, their tracking and conformity between analytically and numerically calculated characteristic currents or characteristic far-field patterns. Moreover, since CMs do not take into account feeding, they can also be used for investigating the accuracy of impedance matrix assemblage which is strongly dependent on discretization [33], the selection of basis functions [34], the quadrature rules used and singularity treatment [35, 36]. As a testing object we propose a perfectly electrically conducting (PEC) sphere for which the characteristic eigenvalues and characteristic eigencurrents are known analytically [1]. The symmetry of the spherical shell also introduces eigenspace degeneration [37] which, together with the null-space of the impedance operator at internal resonances of the shell [34], introduces serious problems with modal tracking [38, 39, 40, 41, 42].

The paper is organized as follows. CM decomposition is briefly recapitulated in Section II and the analytic solution to spherical shell is provided in Section III. The matrix form of CM decomposition is defined in Section IV, including a thorough discussion of the numerical issues behind the decomposition. The benchmarks are proposed in Section V and applied on various packages in Section VI. The paper is concluded in Section VII.

II Characteristic Modes Decomposition

The CMs are introduced [2] as solutions to a generalized eigenvalue problem [43]

[TABLE]

where $\mathcal{R}$ and $\mathcal{X}$ represent the real and imaginary parts of the impedance operator [44]

[TABLE]

with $\boldsymbol{E}$ being the scattered electric field [45], $\boldsymbol{J}_{n}$ the modal current density, $\lambda_{n}$ the characteristic eigenvalue and $\boldsymbol{\hat{n}}$ the unit normal to the PEC surface $\varOmega$ [45] which, in this text, coincides with the radial direction. The time-harmonic quantities under the convention $\boldsymbol{{\mathcal{J}}}\left(\boldsymbol{r},t\right)=\mathrm{Re}\left\{\boldsymbol{J}\left(\boldsymbol{r},\omega\right)\mathrm{exp}\left(\mathrm{j}\omega t\right)\right\}$, with $\omega$ being the angular frequency, are used throughout the paper.

The CMs are commonly normalized with respect to unitary radiated power, i.e.,

[TABLE]

where $\delta_{mn}$ is the Kronecker delta [46], which allows eigenvalues $\lambda_{n}$ to be expressed as [2]

[TABLE]

in which $\mathrm{Re}\left\{\cdot\right\}$ and $\mathrm{Im}\left\{\cdot\right\}$ denote the real and imaginary parts, respectively.

Considering only the currents distributed on surfaces, the uniqueness and completeness of CM decomposition is ensured outside internal resonances [47], i.e., when all modal currents $\boldsymbol{J}_{n}$ radiate. In the light of (1) and (3), a sound definition of characteristic modes can thus be stated as follows: Characteristic modes form a basis of real-valued current densities which diagonalizes the impedance operator ${\mathcal{Z}}$ and possess orthonormal radiation patterns.

III Analytical Decomposition

The analytical solution of CMs decomposition on a spherical shell is presented in this section and the results are to be further used as a reference. We start with a short inspection of systems with potentially known analytical solutions in Section III-A and the CMs of a spherical shell are presented in Section III-B.

III-A Separable Systems

Orthonormality of far-field radiation patterns and completeness are properties shared between CMs and specific solutions to the vector Helmholtz’s equation in separable systems [48]. Particularly, the solutions to the vector Helmholtz’s equation [48] in spherical, conical, rectangular and cylindrical (circular-cylindrical, elliptical-cylindrical, paraboloidal-cylindrical) coordinate systems, orthonormalized with respect to the far-field111Known solutions to the vector Helmholtz’s equation in spheroidal coordinates do not guarantee orthonormality [49]., can be equated to characteristic modes.

The above-mentioned set of possible candidates is further restricted by a practical requirement on the finite extent of the studied structures so that the model can be discretized without using periodic boundary conditions. Consequently, we are left with two feasible coordinate systems – spherical and conical. From these two, we have chosen spherical modes othonormal with respect to spherical surfaces [50]. Such modes correspond to a surface current density distributed on a spherical shell. In comparison to conical solutions, they are formally simpler and, significantly, exhibit high order degeneracies which complicate modal tracking considerably, see Section IV.

III-B CM Decomposition of a Spherical Shell

The analytical form of characteristic currents on a spherical shell, see Fig. 1, can be found in the work of Garbacz [1] where it is a result of diagonalization of a scattering matrix, though without any derivation and with the characteristic numbers $\lambda_{n}$ from (1) presented in a slightly different form ($-1/\kappa_{n}$) which is more favorable for a scattering scenario. Here, instead, we provide a rationale to solve the problem from the perspective of (3) and then present results which can be used as a standard for numerical tests.

The orthonormal set of electric fields [50]

[TABLE]

[TABLE]

corresponding to surface current densities

[TABLE]

[TABLE]

distributed on a spherical shell of radius $r=a$ are the desired solutions to the vector Helmohltz’s equation for $\boldsymbol{r}$ denoting radial direction,

[TABLE]

functions $\boldsymbol{M}_{pq}$ and $\boldsymbol{N}_{pq}$ being defined in [50], $Z_{0}=\sqrt{\mu/\varepsilon}$ being the free-space impedance and $\mathrm{j}_{p}$ ($\mathrm{h}_{p}^{(2)}$) being the spherical Bessel (Hankel) function of $p$-th order [46]. Setting then

[TABLE]

such solutions also satisfy (3) and can be identified with the characteristic modes of a spherical shell, see Fig. 1. Substituting (5), (6), (7) and (8) into (3), the characteristic numbers are found in analytic form as

[TABLE]

and

[TABLE]

where $\mathrm{y}_{p}$ is the spherical Bessel function of the second kind and $p$-th order.

To simplify the notation, an aggregated index $n$ is adopted from [51] as

[TABLE]

and the characteristic numbers, $\lambda_{n}\in\left(-\infty,\infty\right)$, are rescaled in terms of so-called characteristic angles $\delta_{n}\in\left[90^{\circ},270^{\circ}\right]$ as [52]

[TABLE]

The characteristic numbers $\lambda_{n}$ belonging to the first six TE and TM modes (not counting degenerations) are depicted in Fig. 2 and the corresponding characteristic angles $\delta_{n}$ are depicted in Fig. 3.

IV Numerical Evaluation

A numerical solution to the characteristic modes of a spherical shell is found in this section, including a discussion of related numerical issues.

In the common treatment [34], the operator (2) is represented in a basis of piecewise functions $\left\{\,\boldsymbol{\psi}_{u}\right\}$ [34] in the form of the impedance matrix

[TABLE]

and the relation (1) is transformed into the matrix equation

[TABLE]

with modal current density from (1) calculated as

[TABLE]

Formula (16) is of general validity, and, therefore, characteristic modes of arbitrarily shaped bodies can be found in this way at the expense, however, of the occurrence of various numerical issues and artifacts which are discussed below.

IV-A Known Numerical Issues

The numerical issues connected to (16) will be demonstrated on a spherical shell discretized into a triangular mesh grid, see Fig. 4, with RWG basis functions [53] applied. The mesh grid was exported from FEKO [8] as a NASTRAN file [54] and is freely available [55]. Spherical geometry allows a near-perfect triangular mesh consisting of equiangular triangles to be generated, see Fig. 4.

When dealing with the CMs, the following issues arise:

IV-A1 Indefiniteness of a Real Part of the Impedance Matrix

The importance of algebraic properties of matrices $\mathbf{R}$ and $\mathbf{X}$ can readily be seen from (3) and (16). To obtain real characteristic numbers and vectors, these matrices have to be symmetric, i.e., Galerkin testing procedure should properly be applied [44]. If matrices are slightly non-symmetric, they can be symmetrized ex post as

[TABLE]

in which superscript $\mathrm{T}$ denotes matrix transposition. Since matrix $\mathbf{R}$ serves in (16) as a weighing operator [43] and represents radiated power, its positive definiteness, $\mathbf{R}\succ 0$, is crucial. In fact, a potential violation of this condition is one of the biggest known issue related to the CMs [16]. On the contrary, matrix $\mathbf{X}$ is usually well-posed.

In order to present this last problem numerically, matrix $\mathbf{R}$ of a spherical shell discretized into 500 triangles (left panel in Fig. 4) is decomposed as

[TABLE]

and the eigenvalues $\xi_{n}$ are depicted in Fig. 5. The results correspond to matrix $\mathbf{R}$ obtained by two different packages: the commercial package FEKO [8] and the academic tool AToM [11]. For the sake of clarity, the absolute values of the eigenvalues are shown in logarithmic scale and the negativity of the eigenvalues is indicated by the marker (star) used. Note that eigenvalues $\xi_{n}$ are proportional to the radiated power [44].

Generally, we see that only a few modes radiate well (top left part of Fig. 5) and that at least one half of all eigenvalues are negative (right half of Fig. 5). These negative eigenvalues are related to the presence of ubiquitous numerical noise and they should be removed from matrix $\mathbf{R}$ using (19) and substituting $\xi_{n}\equiv 0$ [56] in the consecutive back-composition

[TABLE]

where $\hat{\mathbf{I}}$ is a matrix containing column vectors $\hat{\mathbf{I}}_{n}$ and $\boldsymbol{\xi}$ is a diagonal matrix containing eigenvalues $\xi_{n}$. Unfortunately, formula (20) cannot cure the imminent fact that matrix $\mathbf{R}$ is ill-conditioned.

Matrix $\mathbf{R}$ originating from FEKO is non-symmetric, which causes its different appearance as compared to other curves. For the second set of eigenvalues, the matrix has been symmetrized according to (18) prior to decomposition (19). Two different orders of Gaussian quadrature rule [57] have been used in AToM [44]. In comparison to FEKO, AToM has been able to find more modes, namely those highlighted in the shaded oval denoted by number I. This fact will play an important role at a further point and is most probably caused by single precision arithmetic used in FEKO, see Appendix A for the simulation setup. It is also demonstrated by high negative eigenvalues highlighted by oval II in Fig. 5 that centroid approximation (only one quadrature point in each triangle) does not lead to a well behaved matrix $\mathbf{R}$. It is obvious that the number of properly found modes directly reflect the quality of numerics used when constructing the impedance matrix.

IV-A2 Mode Degeneracy

Degeneracy can be traced back to the geometrical symmetries of the $\varOmega$ region [58, 59, 37] and it poses complications with mode tracking [38, 42, 22]. A spherical shell has a degeneracy [50] of

[TABLE]

which means that for $p$-th order there are $N_{p}$ degenerated modes. Therefore, each solid and dashed line in Fig. 2 and Fig. 3 has $2p+1$ multiplicity. The effect of degeneracy is also seen in Fig. 5 as it appears as stairs spreading their length from left to right. Except for geometrical degeneracies (21), the numerical tracking procedure is yet more complicated for the occasional degeneracies occurring at frequencies where TM and TE modes intersect.

The possibility of reducing the number of geometrically degenerated modes lies in the utilization of the procedure from [59] which, however, relies on a particular choice of basis functions and it is therefore problem-dependent. One useful work-around is to make the mesh grid markedly non-symmetrical which helps to remove the degeneracies.

IV-A3 Internal Resonances

Internal resonances are inherent to all closed PEC surfaces [47, 60, 34] and make the matrix $\mathbf{Z}$ ill-conditioned at resonance frequencies. They occur at those frequencies at which non-radiating current sources may exist, i.e., where denominators of (11) and (12) vanish. The consequences of solving (16) at these frequencies are explained in [60, 61, 23], where it is also important to notice the comparison between resonances of the CMs and natural modes. At the internal resonance, the characteristic numbers abruptly change value between plus and minus infinity, see Fig. 2. The majority of tested packages suffer from this issue, although they can effectively be solved using Combined Field Integral Equation (CFIE) [62].

IV-A4 Dependence on Conformity and Mesh Density

Equation (16) is solved for approximate models like those depicted in Fig. 4 and the question arises of which mesh scheme describes the original smooth object best [34, 31]. Quality of mesh grid, mesh density with respect to conformity, rounding errors, and computational requirements should all be taken into account when dealing with this question. It is also important to note that, although a finer mesh better describes the original object, the increased number of potential modes is redeemed by a higher level of numerical noise resulting in more invalid modes, cf. Fig. 5, longer computational time and higher memory consumption.

IV-A5 Dependence on Integration and Singularity Treatment

As with the previous issue, this point is merely a technicality, yet it strongly influences the final results. The higher order quadrature rules [34] have a great impact on the quality of results, see Fig. 6. Higher-order basis functions can be advantageously applied as well [34]. Special care should also be taken with singularities [35, 36] occurring during the evaluation of (15).

V Benchmarks

Knowing the analytic results and common issues behind characteristic mode decomposition, the following benchmarks are proposed:

Test #1

Characteristic eigenvalues of impedance matrix $\mathbf{Z}$ for a given electrical size $ka$, see Section V-A.

Test #2

Modal tracking in a given range of $ka$, see Section V-B.

Test #3

Conformity between analytically and numerically calculated characteristic currents $\boldsymbol{J}_{n}$, see Section V-C.

Test #4

Correspondence of analytically and numerically calculated characteristic far-fields, see Section V-D.

All the aforementioned tests are extremely simple to implement and add minimal demands on post-processing. The key features to be investigated within all the tests are

• •

choice of the basis and testing functions,

• •

precision of the integration scheme used,

• •

precision and robustness of the singularity treatment used,

• •

quality of the tracking algorithm.

Notice that a fixed mesh is used and that the errors due to meshing are not considered.

V-A Test #1

The first test focuses on the quality of an impedance matrix, thus making it an efficient benchmark of the method of moments codes. The analytically known eigenvalues $\lambda_{n}$ from (11) and (12) are compared with eigenvalues obtained by CM decomposition (16) of impedance matrix $\mathbf{Z}$ from (15), its symmetrized form $\mathbf{Z}_{\mathrm{sym}}$ from (18), or of an impedance matrix with a positively definite real part (20). The test is proposed for electrical sizes $ka=0.5$ and $ka=1.5$, so that both electrically small and reasonably large objects are tested. The spectrum of eigenvalues is calculated by the generalized Schur decomposition [43] (eig with QZ algorithm in Matlab) and all eigenvalues which are infinite, complex-valued or correspond to mode with $\mathbf{I}_{n}^{\mathrm{H}}\mathbf{R}\mathbf{I}_{n}\leq 0$ are removed. The analytical results are depicted in Fig. 7. The maximum number of modes for three different numerical precisions are summarized in Table I.

V-B Test #2

The second test investigates the frequency behavior of characteristic eigenvalues. The setup is as follows: 100 modes should be calculated at 226 equidistantly spaced frequency samples between $ka=0.5$ and $ka=5$, and mode tracking should be provided. The number of modes is selected from Fig. 2, considering mode degeneracy (21) so that, theoretically, all required data to track the first six TE and TM modes are provided ($2\sum_{1}^{6}\left(2p+1\right)=96$). The frequency span is chosen to cover a couple of internal resonances. The number of frequency samples is chosen as a trade-off between sufficient $ka$ sampling ($\Delta ka=0.02$) and a computationally feasible solution. The reference solution is depicted in Fig. 2 and Fig. 3.

V-C Test #3

The conformity of numerically calculated characteristic currents $\mathbf{J}_{n}$ with the analytical results (7) and (8) is studied with the third test using

[TABLE]

and the current $\mathbf{\hat{J}}\left(\mathbf{r}\right)$ is normalized according to

[TABLE]

Ideally, the coefficient $\chi_{n}$ should be equal to unity for all $n$.

V-D Test #4

Far-fields $\boldsymbol{F}_{pq}^{\mathrm{TE/TM}}\left(\vartheta,\varphi\right)$ of analytical characteristic modes of a spherical shell can be deduced from (5), (6) and, following the definition of characteristic modes from Section II, they should form an orthogonal and complete set. Since characteristic radiation patterns are commonly used in practice [40, 63], a meaningful test is to compare numerically evaluated characteristic far-fields $\boldsymbol{F}_{n}$ with analytical ones via

[TABLE]

where

[TABLE]

The metric $\zeta_{n}$ should ideally be equal to unity for all $n$.

As for the previous test, metric $\zeta_{n}$ should mostly judge the quality of numerically evaluated current patterns which are directly reflected in the corresponding far-fields. Simultaneously, it tests the fundamental property of CMs, far-field orthogonality.

VI Results

This section presents the results of four tests from the last section performed on the numerical packages implementing method of moment solution to field integral equations, namely on:

• •

FEKO [8],

• •

CST-MWS [10],

• •

WIPL-D [9],

• •

CEM One [64]

• •

AToM [11],

• •

Makarov [14].

The most important settings of the solvers used are specified in Appendices A–E.

VI-A Test #1

The results of the first test are depicted in Fig. 8 for electric size $ka=0.5$ and in Fig. 9 for electric size $ka=1.5$, respectively. In both cases, the spherical shell was discretized into $500$ discretization elements and, unless otherwise stated, the triangular mesh grid from Fig. 4 was used. The impedance matrices of FEKO and Makarov’s code are originally slightly non-symmetric, therefore (18) was applied before the analysis. The order of Gaussian quadrature in CEM One and AToM can be controlled by the user which is why the selected integration scheme is explicitly mentioned in the parentheses. Data from CST-MWS were not analyzed222The impedance matrix is not accessible and it is thus not guaranteed that the same algorithm for eigenvalue decomposition is used for all tested cases..

The accuracy of the academic package AToM is limited by the double numerical precision used. Therefore, the six and eight lowest order non-degenerated TE and TM modes are correctly found in Fig. 8 and Fig. 9, respectively. The commercial packages are able to represent the first four or five modes, most likely for the use of single precision. Notice that, theoretically, $14$ or $20$ non-degenerated TM and TE modes could be retrieved, see Table I. In the light of this, the number of correctly found modes is relatively low. The lack of numerical precision can be partially compensated using VPA (Variable Precision Arithmetic) in Matlab or the Multiprecision Computing Toolbox [65]. However, these improvements have the potential to add only a couple of additional modes and always at the cost of extreme computational time.

The presented results suggest that the slight non-symmetry of an impedance matrix is not a major issue. However, care should be taken with the numerical precision used and the application of the high-order quadrature rule in conjunction with the precise implementation of singular terms.

VI-B Test #2

This test focuses on the mode tracking algorithm, therefore, only packages with this utility implemented could be analyzed. In all cases, the tracking analysis was performed in the frequency range from $ka=0.5$ to $ka=5$, and then, to produce lucid graphical outputs, only the problematic region of $ka\in\left[2.5,5\right]$ was depicted.

The results from the AToM package are depicted in Fig. 10. The panel (a) shows the raw data acquired from the Implicitly Restarted Arnoldi method [66] (eigs in Matlab). This method was advantageously utilized since only the first $100$ modes are required at each frequency. The gray-colored ellipses highlight two exemplary samples at which the data are missing. The panel (b) shows the modal data after the tracking procedure. A careful inspection reveals a couple of disconnected modes, one incomplete mode (missing data are highlighted by the gray-colored ellipse) and one missing inductive mode (depicted by the red dashed line). The systematic frequency shift between analytically predicted and numerically calculated internal resonances can be attributed to finite meshing is a consequence of the slightly smaller electrical size of the mesh grid.

The results from the FEKO package are depicted in Fig. 11. The raw data in panel (a) are similar to those found with the AToM package, however, the tracked modes in panel (b) are far from perfect.

The next analyzed package is WIPL-D which uses a quadrilateral mesh grid with higher-order basis functions which is why its results cannot be directly compared with other packages. However, the results seem promising, see Fig. 12, except at those places where the crossing-avoidances were incorrectly detected (see the inset).

The last package to undergo testing was CST-MWS. In comparison to the previous packages, CST-MWS uses CFIE for CMs analysis. Only the first four modes were found and, even though $226$ frequency samples were required, $1000$ interpolated values were returned, see Fig. 13. The effect of the interpolation is evident from the enlarged section. The tracking procedure seems to be relatively computationally demanding as only $15$ modes could be calculated and only four modes were found. On the other hand, these modes are well-tracked.

VI-C Test #3

The third test was performed only for the AToM and FEKO packages at $ka=1.5$ for a mesh grid composed of $500$ and $2200$ triangles, see Fig. 14 and Fig. 15. The similarity coefficient $\chi_{n}$ from (22) decreased significantly faster for the poorer mesh grid, however, the number of sufficiently represented modes, say those with $\chi_{n}>0.9$, was similar for both mesh grids. The influence of the higher-order quadrature rule is obvious and the number of well-defined currents, approximately the first $80$ TM and TE modes, including degeneracies, corresponds perfectly with the number of precisely calculated eigenvalues, see Fig. 9.

VI-D Test #4

Test metric $\zeta_{n}$ from (24) was evaluated at electrical size $ka=1.5$ for data obtained from the AToM and FEKO packages, with results depicted in Fig. 16 and Fig. 17. The test was not performed using the CST-MWS package since it internally normalizes characteristic modes via $L^{2}$ norm and not via (3). The results could, therefore, not be compared with others. The CEM One package does not have characteristic far-fields implemented.

For proper assessment, the results in both figures should be compared simultaneously since only simultaneous compliance with $\zeta_{n}\approx 1$ and $P_{\mathrm{r},n}\approx 1$ is sufficient to pass the test, with $P_{\mathrm{r},n}$ being the far-field radiated power of the $n$-th characteristic mode, i.e., the denominator of (24). The normalization of radiated power may seem to be automatic, but for higher order modes the normalization (3), performed numerically via $\mathbf{I}^{\mathrm{H}}\mathbf{R}\mathbf{I}=2$, does not imply unitary radiated power in far-field.

Only modes with real eigenvalues and positive radiated power were taken in the case of AToM. Data from FEKO were left in their original form. In all cases, the integration (25) was discretized in $N_{\vartheta}=N_{\varphi}=100$ points. This angular discretization was checked by artificially testing analytical modes which resulted in $\zeta_{n}\approx 1$ and $P_{\mathrm{r},n}\approx 1$ for all analytical modes within the depicted range.

The performance in all tested cases is quite unsatisfactory, but generally corresponds to the other tests. The results suggest that far-field decomposition into more than five TE modes and five TM modes (not including degenerations) of a spherical shell is unsafe within the used triangularization of a sphere and double precision. Higher quadrature rules and careful singularity treatment, seen in the case of AToM with an $8$-th order of quadrature, add more proper modes, but the gain is not as high as in other tests.

VII Conclusion

A set of sensitive benchmarks applicable to the majority of contemporary method of moments packages was proposed. The major advantages are formal simplicity, the low number of necessary inputs, the independence on particular feeding model and minimal required post-processing.

All tested packages showed satisfactory performance in calculating characteristic eigenvalues at a single frequency, although some of them were limited out of the single numerical precision range. Without exception, however, they performed well below theoretically achievable results. In the case of modal tracking, the differences were more severe and some packages showed unsatisfactory results. It has been demonstrated that factors, such as singularity treatment, high-order quadrature rules and used floating-point numerical precision, greatly influence the results. A persisting problem is also the numerical ill-posedness of the real part of impedance matrix caused by dominance of only couple of modes. The application of any technique increasing its dynamical range while preserving its precision is, therefore, of interest since it may significantly improve the number of correctly found characteristic modes.

Future work should be aimed toward analytical characteristic modes of more complex shapes which could further stress the precision of available method of moments codes. Promising candidates are surfaces which can be described in conical or spheroidal coordinate systems.

Appendix A FEKO Setup

FEKO (ver. 14.0-273612, [8]) was used with the following settings: a mesh structure was imported into the software using the Nastran file format [54] and the CMs and far-fields were chosen as a request for the FEKO solver. Data from FEKO were acquired using *.out, *.os, *.mat and *.ffe files. The impedance matrices were imported by an in-house wrapper [13]. Double precision was enabled for data storage in solver settings. The number of modes used for each test were as follows: Test #2 $99$ modes, Test #3 $500$ modes and Test #4 $300$ modes, respectively.

Appendix B WIPL-D Setup

WIPL-D (v13, [9]) uses higher-order basis functions with a quadrilateral mesh grid. Therefore, individual meshing was used. Surface angle tolerance was set to $15$ degrees. Additionally, a non-symmetrical mesh was created in order to obtain better stability of mode calculations. An integral solver with CM decomposition was utilized with double precision and enhanced-1 for Integral Accuracy. The matrices for Test #1 were delivered by the developer of the WIPL-D package and data for Test #2 included first $100$ modes.

Appendix C CST-MWS Setup

CST-MWS (ver. 2016.7, [10]) setting was as follows: the spherical shell was created using a sphere modeling tool and converted into a sheet. The parameter cells per wavelength was set to $5.25$ in general mesh properties which led to a mesh grid of $500$ triangular elements. The Integral equation solver was chosen and the number of requested modes was set to 16. It is worth noting that regardless of the choice of number of frequency samples in CMA accuracy settings, CST-MWS always used $1001$ frequency samples.

Appendix D CEM One Setup

CEM One (2015.2, [64]) setting was as follows: a mesh grid was imported using the Nastran file format [54]. The impedance matrix was obtained by the “Save System Matrix” in the Final Output Parameters. Quadrature order was set in a *.dat file as:

#RUMSEY simple 3 6 double 3 6

The first number represents the quadrature order for far elements and the second number represents the quadrature order for singularities. The impedance matrix was saved using a ncdump command in the E-Field command Prompt.

Appendix E AToM Setup

AToM (pre-product ver.) setting was as follows: a mesh grid was imported using the Nastran file format [54]. The number of modes was set manually to $100$ for Test #2 (the eigs routine was used) and kept as the maximum for the other tests (the eig routine). In AToM, RWG basis functions with the Galerkin procedure are adapted [53]. The Gaussian quadrature is implemented according to [57] and singularity treatment is implemented from [35]. The tracking of the CMs follows [38], including some recently added adaptive strategies [42]. Since AToM is written in Matlab, no other import procedures were needed.

Acknowledgement

The authors would like to thank Vladimir Sedenka from BUT and Michal Masek from CTU in Prague for implementing the method of moments and characteristic modes solvers within the AToM package, and would also like to thank to Gerhard Kristensson from Lund University for a fruitful discussion concerning separable systems, Kurt Schab from NC State University, Doruk Tayli from Lund University and Jasmin Music from WIPL-D for providing the impedance matrices.