Magnetic and Combined Field Integral Equations Based on the Quasi-Helmholtz Projectors

TL;DR

This paper introduces new electromagnetic boundary integral equations based on quasi-Helmholtz projectors that address common numerical issues like ill-conditioning and spurious resonances, improving solution stability across frequencies.

Contribution

The paper develops a novel quasi-Helmholtz projector-based formulation for magnetic and combined field integral equations that are immune to multiple numerical drawbacks.

Findings

The new formulations are immune to low-frequency ill-conditioning.

They effectively eliminate spurious resonances.

Numerical results demonstrate improved stability and accuracy.

Abstract

Boundary integral equation methods for analyzing electromagnetic scattering phenomena typically suffer from several of the following problems: (i) ill-conditioning when the frequency is low; (ii) ill-conditioning when the discretization density is high; (iii) ill-conditioning when the structure contains global loops (which are computationally expensive to detect); (iv) incorrect solution at low frequencies due to current cancellations; (v) presence of spurious resonances. In this paper, quasi-Helmholtz projectors are leveraged to obtain a magnetic field integral equation (MFIE) formulation that is immune to drawbacks (i)-(iv). Moreover, when this new MFIE is combined with a regularized electric field integral equation, a new quasi-Helmholtz projector combined field integral equation is obtained that also is immune to (v). Numerical results corroborate the theory and show the practical impact of the newly proposed formulations.

Equations

$$\left({\mathcal{T}}_k\mathbf{j}\right)\left(\mathbf{r}\right)$$

$$\left({\mathcal{T}}_{s,k}\mathbf{j}\right)\left(\mathbf{r}\right)$$

$$\left({\mathcal{T}}_{h,k}\mathbf{j}\right)\left(\mathbf{r}\right)$$

$$\left({\mathcal{K}}_k\mathbf{j}\right)\left(\mathbf{r}\right)$$

$$\left({\mathcal{T}}_k\mathbf{j}\right)(\mathbf{r})=\hat{\mathbf{n}}\times {\mathbf{e}}^i(\mathbf{r})$$

$$\left(\left(\frac{\mathcal{I}}{2}+{\mathcal{K}}_k\right)\mathbf{j}\right)(\mathbf{r})=\hat{\mathbf{n}}\times {\mathbf{h}}^i(\mathbf{r})$$

$$\mathbf{j}(\mathbf{r})\approx \sum _{j=1}^{N_e}{\left[\text{mathsfbfit}j\right]}_j{\mathbf{f}}_j(\mathbf{r}),$$

$$\text{mathsfbfit}T\text{mathsfbfit}j$$

$$\left(\frac{\text{mathsfbfit}{G}^{\mathrm{T}}}{2}+\text{mathsfbfit}{K}_k\right)\text{mathsfbfit}j$$

$${\left[\text{mathsfbfit}T\right]}_{ij}$$

$${\left[\text{mathsfbfit}{T}_s\right]}_{ij}$$

$${\left[\text{mathsfbfit}{T}_h\right]}_{ij}$$

$${\left[\text{mathsfbfit}{K}_k\right]}_{ij}$$

$${\left[\text{mathsfbfit}G\right]}_{ij}$$

$${\left[\text{mathsfbfit}{v}_e\right]}_i$$

$${\left[\text{mathsfbfit}{v}_h\right]}_i$$

$$\text{mathsfbfit}j=\text{mathsfbfit}\Lambda \text{mathsfbfit}l+\text{mathsfbfit}\Sigma \text{mathsfbfit}s+\text{mathsfbfit}H\text{mathsfbfit}h$$

$$\text{mathsfbfit}{\Lambda }_{ij}$$

$$\text{mathsfbfit}{\Sigma }_{ij}$$

$$\text{mathsfbfit}\Lambda 1^{\text{mathsfbfit}\Lambda }=0.$$

$$\text{mathsfbfit}\Sigma 1^{\text{mathsfbfit}\Sigma }=0.$$

$$\text{mathsfbfit}{\Sigma }^{\mathrm{T}}\text{mathsfbfit}\Lambda =0.$$

$$\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }$$

$$\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}$$

$$\text{mathsfbfit}j=\left(\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}\text{mathsfbfit}j\right)+\left(\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }\text{mathsfbfit}j\right)$$

$$\text{mathbbmsl}{P}^{\text{mathsfbfit}\Lambda }$$

$$\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H}$$

$$\left(\frac{\mathcal{I}}{2}-{\mathcal{K}}_{-\mathrm{j}k}\right)\left(\frac{\mathcal{I}}{2}+{\mathcal{K}}_k\right)\left(\mathbf{j}\right)=\left(\frac{\mathcal{I}}{2}-{\mathcal{K}}_{-\mathrm{j}k}\right)\left({\hat{\mathbf{n}}}_r\times {\mathbf{h}}^i\right).$$

$$\text{mathbbmsl}{M}^{\mathrm{T}}\left(\frac{\text{mathsfbfit}{G}^{\mathrm{T}}}{2}-\text{mathsfbfit}{K}_{-\mathrm{j}k}\right){\left(\text{mathsfbfit}{G}^{\mathrm{T}}\right)}^{-1}\left(\frac{\text{mathsfbfit}{G}^{\mathrm{T}}}{2}+\text{mathsfbfit}{K}_k\right)\text{mathsfbfit}M\text{mathsfbfit}i$$

$$\text{mathsfbfit}M$$

$$\text{mathbbmsl}M$$

$$\begin{gathered} \text{mathsfbfit}O=(\text{mathbbmsl}{P}^{\text{mathsfbfit}\Lambda }+\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H})\text{mathsfbfit}O(\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }+\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H})= \\ \text{mathbbmsl}{P}^{\text{mathsfbfit}\Lambda }\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }+\text{mathbbmsl}{P}^{\text{mathsfbfit}\Lambda }\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}+\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H}\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }+\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H}\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}, \end{gathered}$$

$$\text{mathbbmsl}{P}^{\text{mathsfbfit}\Lambda }\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }$$

$$\text{mathbbmsl}{P}^{\text{mathsfbfit}\Lambda }\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}$$

$$\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H}\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Sigma }$$

$$\text{mathsfbfit}{K}_k=\text{mathsfbfit}{K}_0+\text{mathsfbfit}{K}_k^{\prime },$$

$$\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H}\left(\frac{\text{mathsfbfit}{G}^{\mathrm{T}}}{2}-\text{mathsfbfit}{K}_0\right){\left(\text{mathsfbfit}{G}^{\mathrm{T}}\right)}^{-1}\left(\frac{\text{mathsfbfit}{G}^{\mathrm{T}}}{2}+\text{mathsfbfit}{K}_0\right)\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}=0.$$

$${\alpha }^2\text{mathbbmsl}{P}^{\text{mathsfbfit}\Sigma \text{mathsfbfit}H}\text{mathsfbfit}O\text{mathsfbfit}{P}^{\text{mathsfbfit}\Lambda \text{mathsfbfit}H}=$$

Full text

\acresetall

I Introduction

Time-harmonic scattering by \acPEC objects oftentimes is modeled using frequency-domain boundary integral equations. Among them, electric and magnetic field integral equations (\acsEFIE and \acsMFIE) [1] are the most popular.

Although the \acEFIE is easily discretized using \acRWG basis functions [2], it suffers from ill-conditioning when the frequency is low and/or the discretization density is high. The \acMFIE, on the other hand, remains well-conditioned in both regimes, provided that a mixed discretization scheme is employed [3]. In practice, however, it is not feasible to obtain accurate results for the \acMFIE at extremely low frequencies without resorting to highly precise numerical quadrature methods. In addition to the above issues, both the \acEFIE and the \acMFIE suffer from current cancellations at low frequencies [4, 5, 6].

The \acEFIE’s conditioning and current cancellation problems can be overcome by using loop-star or loop-tree decompositions [7, 8, 9, 10, 11]. For multiply connected geometries, this requires the detection of global loops, which is computationally expensive [12]. These techniques also fail to address the dense discretization breakdown phenomena [13, 14] which causes the \acEFIE’s condition number to grow quadratically with the mesh refinement parameter. Worse still, loop-star techniques for combatting the \acEFIE’s low-frequency conditioning problems further degrade the equations dense discretization behavior [15].

Several formulations have been introduced to address these low-frequency issues without the computational burden of global loop detection [16, 17]. These solutions, however, do not address the dense discretization ill-conditioning of the \acEFIE. Both issues can be concurrently tackled by leveraging hierarchical quasi-Helmholtz decompositions [18, 19, 20, 13]. These decompositions also have been successfully coupled with other approaches such as Calderón preconditioning [21, 22, 23, 24, 25, 14, 26, 27, 28, 29] and Debye-inspired schemes [30]. The price to be paid for this dual stabilization is, once again, the need for global loop detection at very low frequencies. In addition, several of the aforementioned techniques fail to properly address low-frequency numerical cancellations occurring in the solution vector [4, 31, 5, 32]. Several of the above drawbacks have been successfully addressed by the promising scheme in [33]. Alternative remedies to current cancellations include perturbation methods [4, 31, 34] and Calderón regularization combined with loop star decompositions [25, 26]. Both families of solutions do, however, have shortcomings: the former is only applicable at low frequencies and exhibits the same spectral issues as the formulation it is applied to – high refinement breakdown for the \acEFIE or global loop detection for the \acMFIE and Calderón \acEFIE – while the latter also requires global loop detection and treatment of the high refinement instability of the loop-star decomposition. It should also be noted that some recent incarnations of augmented equations are immune to several of the above mentioned drawbacks, though they require the recovery of auxiliary quantities [35, 36].

Recently, an electric type equation based on quasi-Helmholtz projectors was proposed that is immune to all of the aforementioned issues [37]. A similar regularization has also been applied to the time domain electric field integral equation [38, 39] and both the time domain and the frequency domain PMCHWT equations [40, 41].

In this paper, quasi-Helmholtz projectors are used to obtain a new \acMFIE that no longer requires interaction integrals to be computed using extremely accurate quadrature rules. Additionally, the solenoidal and nonsolenoidal current components are scaled such that low frequency cancellations are avoided. As a result, the formulation remains accurate down to extremely low frequencies. Scattering problems involving \acPEC objects can also be solved using the \acCFIE, which is a linear combination of the \acEFIE and the \acMFIE. This equation has the added benefit that it does not support spurious resonances [42]. In this paper, the new regularization method for the \acMFIE is combined with that for the \acEFIE presented in [37]. The resulting \acCFIE is not only low-frequency stable but also immune to spurious resonances. Preliminary results of this research have previously been presented as conference contributions [43, 44].

This paper is organized as follows. To set notation, Section II defines the standard \acEFIE and \acMFIE as well as their discretizations and related quasi-Helmholtz current decompositions. In Section III, a quasi-Helmholtz decomposition is applied to a new symmetrized form of the \acMFIE. The resulting equation can be discretized accurately using standard numerical quadrature methods, and can be scaled in frequency such that no low frequency cancellations occur. In Section IV, this \acMFIE is combined with the regularized \acEFIE [37] to obtain an extremely low frequency stable \acCFIE. Section V discusses numerical results that corroborate the theory and conclusions are presented in Section VI.

II Background and Notations

The \acEFIE and \acMFIE operators $\bm{\mathcal{T}}_{k}$ and $\bm{\mathcal{K}}_{k}$ are defined as

[TABLE]

where $R=\|\bm{r}-\bm{r}^{\prime}\|$, $\Gamma$ is the boundary of a closed domain $\Omega\subset\mathbb{R}^{3}$ and $\hat{\bm{n}}$ is its exterior normal vector. Furthermore, given the angular frequency $\omega$, $k=\omega\sqrt{\mu\epsilon}$ and $\eta=\sqrt{\mu/\epsilon}$; here $\epsilon$ and $\mu$ the permittivity and permeability of vacuum, respectively. If $\Omega$ is perfectly conducting, it supports an electric current $\bm{j}(\bm{r})$ satisfying both the \acEFIE

[TABLE]

and the \acMFIE

[TABLE]

for all $\bm{r}\in\Gamma$; where $\bm{e}^{i}$ and $\bm{h}^{i}$ denote the impinging electric and magnetic fields, respectively. To numerically solve these equations via a Galerkin procedure, $\bm{j}(\bm{r})$ is expanded into RWG basis functions $\{\bm{f}_{j}(\bm{r})\}$ [2] as

[TABLE]

where $N_{e}$ is the number of edges of the mesh. Following [37], the RWG functions are normalized such that the integrated flux through their defining edges equals one.

Next, the \acEFIE (5) is tested with rotated RWG functions $\{\hat{\bm{n}}\times\bm{f}_{i}(\bm{r})\}$, while the \acMFIE (6) is tested with rotated \acBC functions [45] $\{\hat{\bm{n}}\times\bm{g}_{i}(\bm{r})\}$. The \acBC functions $\{\bm{g}_{j}\}$ are divergence-conforming functions defined on the barycentric refinement of the mesh. In addition, they are quasi curl-conforming in the sense that the mixed Gram matrix between curl-conforming rotated \acBC functions and \acRWG functions is well conditioned. For an explicit definition of these functions the reader is referred to [14, 45]. Overall, the testing procedure results in the following matrix equations:

[TABLE]

where

[TABLE]

with $\left(\bm{a},\bm{b}\right)=\int_{\Gamma}\bm{a}(\bm{r})\cdot\bm{b}(\bm{r})\mathrm{d}s$. In addition we denote by $\mathbbmsl{T}$, $\mathbbmsl{T}_{s}$ and $\mathbbmsl{T}_{h}$ the \acBC-expanded and tested counterparts of the discretized operators $\mathsfbfit{T}$, $\mathsfbfit{T}_{s}$ and $\mathsfbfit{T}_{h}$ computed with the complex wavenumber $-\mathrm{j}k$.

The solutions of (8) and (9) can be expressed as linear combinations of divergence free (loop and harmonic functions) and of non-divergence free (star functions) contributions via a quasi-Helmholtz decomposition

[TABLE]

where the first two matrices $\mathsfbfit{\Lambda}\in\mathbb{R}^{N_{e}\times N_{v}}$ and $\mathsfbfit{\Sigma}\in\mathbb{R}^{N_{e}\times N_{f}}$ represent mappings from the RWG subspace to the local loop and star subspaces, respectively. Here, $N_{v}$ and $N_{f}$ are the number of vertices and facets of the mesh, respectively [12, 46]. These two mappings can be defined using only the connectivity information of the discretized geometry as

[TABLE]

where vertices $\bm{v}_{i}^{-}$ and $\bm{v}_{i}^{+}$ define the oriented edge characterizing RWG function $i$, and $c_{i}^{-}$ and $c_{i}^{+}$ denote the corresponding cells (Figure 1). The matrix $\mathsfbfit{H}$ represents the mapping from the RWG space to the quasi-harmonic or global loop space composed of $2N_{h}$ functions, where $N_{h}$ is the number of handles in the structure. For a complete description of this mapping and the associated harmonic functions, the reader is referred to [12] and [46].

A few properties of these matrices are recalled next to facilitate further developments. For the sake of simplicity we restrict ourselves to the case of a geometry with a single closed connected component. All derivations below can be extended to arbitrary geometries using the relations in [47]. Given this assumption, $\mathsfbfit{\Lambda}$ has a null-space spanned by the all-one vector $\bm{1}^{\mathsfbfit{\Lambda}}\in\mathbb{R}^{N_{v}}$, i.e.

[TABLE]

Similarly, linear dependency of the star functions cause $\mathsfbfit{\Sigma}$ to exhibit a one-dimensional null-space spanned by the all-one vector $\bm{1}^{\mathsfbfit{\Sigma}}\in\mathbb{R}^{N_{f}}$, i.e.

[TABLE]

Finally, it is trivial to show that the loop and star subspaces are orthogonal, i.e.

[TABLE]

As $\mathsfbfit{\Lambda}$ and $\mathsfbfit{\Sigma}$ are ill-conditioned and because of the high computational cost of detecting global loops required to build $\mathsfbfit{H}$, it is convenient to leverage the quasi-Helmholtz projectors introduced in [37] to obtain a quasi-Helmholtz decomposition of the \acEFIE and \acMFIE operators. The projectors are defined as

[TABLE]

where + denotes the Moore-Penrose pseudo-inverse and $\mathsfbfit{\mathrm{I}}$ is the identity. Any RWG expansion coefficient vector can then be decomposed as

[TABLE]

where $\mathsfbfit{P}^{\mathsfbfit{\Lambda}\mathsfbfit{H}}\mathsfbfit{j}$ and $\mathsfbfit{P}^{\mathsfbfit{\Sigma}}\mathsfbfit{j}$ contain the RWG expansions of the solenoidal (loop) and non-solenoidal (star) components of the current, respectively. These operators are self-adjoint and also can be used to decompose the RWG testing space. Similarly, the dual projectors $\mathbbmsl{P}^{\mathsfbfit{\Lambda}}$ and $\mathbbmsl{P}^{\mathsfbfit{\Sigma}\mathsfbfit{H}}$, defined as

[TABLE]

decompose any linear combination of \acBC (basis or testing) functions into a non-solenoidal and solenoidal part, respectively. It should be noted that construction of these projectors does not require the detection of global loops, and that $\left(\mathsfbfit{\Sigma}^{\mathrm{T}}\mathsfbfit{\Sigma}\right)^{+}$ can be efficiently computed using multigrid preconditioners [37, 48].

III Regularizing the \acMFIE at Extremely Low Frequencies

III-A Low Frequency Behaviour of the \acMFIE

The standard RWG discretization of the \acMFIE fails to provide accurate results at low frequencies due to the unphysical scaling of the loop and star (or tree) components of the current [4]. It was shown in [49, 3] that the mixed discretization of the \acMFIE (in which \acBC or CW functions [50] are used as testing functions) improves accuracy. In particular, the loop and star components of the current obtained from this formulation scale physically [32]. This result also holds true for multiply connected geometries [51].

The mixed \acMFIE formulation still suffers from three problems. First, the physical scaling of the current can only be retrieved when interaction integrals are computed to high accuracy [32]. Second, the nonsolenoidal current component scales as $\mathcal{O}(\omega)$ whereas the solenoidal component is of $\mathcal{O}(1)$. As a result, at very low frequencies and when using finite precision, both components should be stored in different arrays to prevent the nonsolenoidal component from losing accuracy or even being cancelled out [4, 31, 5, 52, 51]. Third, the static \acMFIE (at $\omega=0$) has a null space when applied to multiply connected geometries. It follows that the discretized \acMFIE has $N_{h}$ singular values that scale as $\mathcal{O}(\omega^{2})$ [46]. Any accurate discretization of the \acMFIE operator must preserve this null-space. Standard RWG discretizations of the \acMFIE operators are not capable of correctly modelling this null space [43]. The mixed \acMFIE, on the other hand, correctly models this-null space in infinite precision. However, after discretization, the null-space associated singular values will not be more accurate than the precision of the quadrature rule.

III-B A Robust \acMFIE Formulation

To address the above described \acMFIE deficiencies we propose the following symmetrized MFIE:

[TABLE]

This equation is the magnetic field counterpart of the (localized) Calderón preconditioned electric operator in [37]. We propose to discretize (28) as

[TABLE]

where

[TABLE]

and $\mathsfbfit{M}\mathsfbfit{i}=\mathsfbfit{j}$.

The coefficient $\alpha$ allows for re-scaling of the loop and star components of the solution $\mathsfbfit{i}$ of (29) to prevent numerical cancellations. Because $\mathsfbfit{P}^{\mathsfbfit{\Sigma}}+\mathsfbfit{P}^{\mathsfbfit{\Lambda}\mathsfbfit{H}}=\mathbbmsl{P}^{\mathsfbfit{\Lambda}}+\mathbbmsl{P}^{\mathsfbfit{\Sigma}\mathsfbfit{H}}=\mathsfbfit{\mathrm{I}}$, operator $\mathsfbfit{O}$ in (29) can be decomposed as

[TABLE]

which allows for the study of the low-frequency behavior of each of the separate terms. Analysis of the frequency behavior of the first three terms is quite straightforward and yields

[TABLE]

Analysis of the last term in (32) requires special care. It is known that when decomposing $\mathsfbfit{K}_{k}$ as

[TABLE]

where $\mathsfbfit{K}_{0}$ is the static limit of $\mathsfbfit{K}_{k}$ and $\mathsfbfit{K}^{\prime}_{k}=\mathsfbfit{K}_{k}-\mathsfbfit{K}_{0}$ is the dynamic remainder, $\mathsfbfit{K}^{\prime}_{k}=\mathcal{O}(k^{2})$ as $k\to 0$ [51]. When using this decomposition in (29), it can be verified that $\mathsfbfit{K}_{0}$ satisfies

[TABLE]

The above equation holds the key to unlocking a frequency-stable \acMFIE. The proof of property (35) is provided in Appendix A. The term $\mathbbmsl{P}^{\mathsfbfit{\Sigma}\mathsfbfit{H}}\mathsfbfit{O}\mathsfbfit{P}^{\mathsfbfit{\Lambda}\mathsfbfit{H}}$ can now be studied. To this end, note that

[TABLE]

which completes the low-frequency analysis of the overall operator

[TABLE]

To choose $\alpha$, in addition to the conditioning constraint imposed by (37), we need to consider the physical scaling of the current, which for a plane wave excitation, is [5]

[TABLE]

These scaling laws reveal that for a standard formulation, a severe numerical cancellation is expected due to the fact that the non-solenoidal component of the current (which scales as $\mathcal{O}(k)$) will disappear when stored alongside the solenoidal component (which scales as $\mathcal{O}(1)$). Instead, for the regularized formulation proposed here, the equation is solved for $\mathsfbfit{i}=\mathsfbfit{M}^{-1}\mathsfbfit{j}$, which scales as

[TABLE]

It is now evident that by setting $\alpha=\sqrt{k}$, the above scaling behaviors become

[TABLE]

eliminating the low frequency cancellation and, at the same time, stabilizing the matrix at low frequencies. The latter is seen upon inserting the new scalings into (37):

[TABLE]

The deficiency of the MFIE in the static regime also is solved by the scheme proposed here. In fact, using (37) when $k=0$ we obtain

[TABLE]

which proves the existence of an exact matrix null-space in statics of dimension exactly equal to that of the harmonic subspace.

Summarizing, the proposed \acMFIE resolves the three main issues of prior standard and non-standard \acMFIE formulations and now can be linearly combined with \acpEFIE using projectors.

IV A new \acCFIE

The theoretical developments of the previous sections resulted in a magnetic field operator that can be stably discretized for arbitrarily low frequencies using standard integration rules. The electric counterpart of this operator was obtained in [37]. We will now combine these two operators, first proving the resonance-free property of their continuous combination at high frequencies, and then showing their compatibility at arbitrarily low frequencies.

Standard Calderón \acCFIE equations use a localization strategy for the \acEFIE component to obtain a resonance-free equation [23, 22]. Here, we follow the Yukawa-Calderón approach in [22]. When the Yukawa-Calderón \acEFIE is linearly combined with the new magnetic operator defined in Section III, the following symmetric Yukawa-Calderón \acCFIE is obtained:

[TABLE]

To demonstrate that this equation represents a valid Calderón CFIE, i.e. is free from internal resonances, we prove in Appendix B that the operator

[TABLE]

can be inverted for any $k$.

The discretization of the proposed Yukawa-Calderón \acCFIE follows directly from that of the new \acMFIE in Section III and that of the \acEFIE in [37]:

[TABLE]

Here $\alpha=1$ and $\alpha=\sqrt{k}$ in the high and low frequency regime, respectively. We next study the latter more in detail.

Scaling in the latter regime follows from the results of the previous section:

[TABLE]

Combining this result with the corresponding \aclRHS scalings (42) and (43) proves the overall low-frequency stability of new \acCFIE.

V Numerical Results

This section presents numerical results that validate the above properties of the proposed \acMFIE and \acCFIE.

The first set of tests involve \iacPEC sphere of radius $1\text{\,}\mathrm{m}$. Figure 2 shows the scattered far field at $f=$200\text{,}\mathrm{MHz} obtained using the new MFIE and CFIE as well as other established formulations (standard \acEFIE, \acEFIE with projectors, Calderón \acEFIE with projectors, Mixed \acMFIE, \acCFIE). For this high frequency case all formulations deliver accurate results, thus validating our implementations. A first difference in performance between our new formulations and their standard counterparts is noted when lowering the frequency. [Figure 3](#S5.F3) shows data similar to [Figure 2](#S5.F2) but for $f=$1\text{\times}{10}^{-40}\text{\,}\mathrm{Hz}. It is clear that accuracy breakdowns occur for the non-projected methods – the mixed \acMFIE, the \acEFIE, and the \acCFIE (for the latter two the lack of accuracy also is due to conditioning problems). On the other hand, all projected formulations, including the two new ones, deliver accurate results for arbitrarily low frequencies.

The low frequency stability of the new Calderón \acMFIE is further demonstrated in Figure 4, which illustrates the conditioning of the different operators for low frequencies. It is clear that the new \acMFIE remains as well-conditioned as its standard counterpart. The Calderón \acCFIE is also low-frequency stable, unlike the standard \acCFIE, which exhibits a severe ill-conditioning caused by its \acEFIE contribution.

Figure 5 shows that, despite its regularized low frequency behavior, the Calderón \acMFIE is prone to spurious resonances causing it to become periodically ill-conditioned. This issue is shared by all non-combined formulations and can be overcome by combined field strategies. It is clear from the figure that both the new Calderón \acCFIE and its standard counterpart exhibit resonance-free behaviour.

The last key property to be illustrated is the refinement stability of the proposed formulations. This property was verified by studying the dependence of the condition number of the different formulations applied to a unit radius sphere with increasing discretization density (Figure 6). These results confirm that the second kind nature of our new formulations renders them immune to the high-refinement breakdown.

In summary, the above results show that the new Calderón \acMFIE yields correct results for arbitrarily low frequencies and is well conditioned for both low frequencies and dense discretization. Additionally, when combined with the projector Calderón \acEFIE the new Calderón \acCFIE, which is low frequency stable, immune to dense discretization breakdown, and free from non-physical resonances, is obtained.

To ensure that the properties illustrated so far still persist for multiply connected structures, many of the previous analyses were repeated for a square torus. The correctness of the formulation has been verified by studying the far field scattered by the torus at high and very low frequencies, respectively (Figures 7 and 8). Since no analytic solution is readily available for the square torus, the solution of the Calderón \acEFIE was used as a reference and particular care was taken to avoid frequencies corresponding to an internal resonance. While the results are similar to those of the sphere, the reader should be aware that, because of its toroidal and poloidal null-spaces, the Calderón \acMFIE required the usage of a pseudo inversion to obtain current solutions at very low frequencies.

The low frequency stability of the Calderón \acMFIE and Calderón \acCFIE on the toroidal structure are demonstrated in Figure 9, while their resonance free behaviors are illustrated in Figure 10. Finally, the resilience of both formulations to dense discretization breakdown is illustrated in Figure 11, which presents the condition number of the integral operators with increasing discretization of the square torus.

One of the key advantages of the new Calderón \acMFIE scheme is that it does not require extremely accurate numerical integration rules because it allows explicit cancellation of near-zero terms that are challenging to obtain numerically. The slow convergence of the standard numerical integration schemes can be seen in Figure 12, in which the ratio of the norm of the term in (50) to the norm of the full operator with increasing number of integration points is presented. While this ratio does decrease with the number of Gaussian quadrature points, it does so very slowly and remains far from a machine-precision zero value. The effect of these numerical inaccuracies is evident when comparing the singular value decompositions of the Mixed \acMFIE and of the new Calderón \acMFIE in Figure 13. It is clear that the null singular values corresponding to the toroidal and poloidal subspaces of the square torus immediately reach the machine precision zero in the case of the Calderón \acMFIE, while for the Mixed \acMFIE they will require an unreasonably complex integration rules to even remotely resemble a nullspace.

Finally, to demonstrate that our schemes can be readily applied to more complex problems we studied the low frequency conditioning of our operators (Figure 15) for the complex, multiply connected geometry in Figure 14.

VI Conclusion

This paper presented a new symmetrized MFIE that can be stably and effectively discretized using quasi-Helmholtz projectors. When linearly combined with a quasi-Helmholtz projector-based Calderón \acEFIE, a new \acCFIE that is immune from all drawbacks that plague the majority of existing formulations is obtained. In fact, the proposed CFIE remains well-conditioned both at low frequencies and for high discretization densities, allows for an accurate solution at extremely low frequencies without requiring special numerical quadrature methods, does not require the detection of global loops when applied to multiply connected geometries, and is provably free from internal resonances. Numerical results confirm the theoretically predicted properties of the proposed equations.

Appendix A Proof of the fundamental matrix relationship (35)

To prove the validity of (35), i.e.

[TABLE]

we introduce $\mathsfbfit{P}^{\text{Pol}}$, $\mathsfbfit{P}^{\text{Tor}}$, $\mathbbmsl{P}^{\text{Pol}}$, $\mathbbmsl{P}^{\text{Tor}}$ the orthogonal projectors into the right and left null-spaces of the internal and external \acMFIE operators, i.e.

[TABLE]

Note that

[TABLE]

We can then define

[TABLE]

which clearly satisfies

[TABLE]

since the union of the right null spaces of the internal and external \acMFIE operators contains all the non-trivial cycles of the structure [46]. Dually,

[TABLE]

satisfies

[TABLE]

It follows that

[TABLE]

and similarly that

[TABLE]

Combining the above equations it follows that

[TABLE]

In the above expression we now insert the identity matrices $\left(\mathsfbfit{P}^{\mathsfbfit{\Lambda}\mathsfbfit{H}}+\mathsfbfit{P}^{\mathsfbfit{\Sigma}}\right)$ and $\left(\mathbbmsl{P}^{\mathsfbfit{\Sigma}\mathsfbfit{H}}+\mathbbmsl{P}^{\mathsfbfit{\Lambda}}\right)$ obtaining

[TABLE]

Given that

[TABLE]

and that

[TABLE]

and considering the property

[TABLE]

we obtain that

[TABLE]

which completes the proof.

Appendix B Resonance-free proof for the new Calderón \acCFIE operator

Since the operator $\left(\frac{\bm{\mathcal{I}}}{2}-\bm{\mathcal{K}}\right)(-\mathrm{j}k)$ always admits an inverse, the invertibility of (47) is equivalent to the invertibility of

[TABLE]

Given the anti-commutation property

[TABLE]

which follows directly from the second Calderón identity $\bm{\mathcal{T}}^{-1}\bm{\mathcal{K}}=\bm{\mathcal{T}}^{-1}\bm{\mathcal{K}}\bm{\mathcal{T}}\bm{\mathcal{T}}^{-1}=-\bm{\mathcal{T}}^{-1}\bm{\mathcal{T}}\bm{\mathcal{K}}\bm{\mathcal{T}}^{-1}=-\bm{\mathcal{K}}\bm{\mathcal{T}}^{-1}$, and defining

[TABLE]

it follows that

[TABLE]

Given this result and the fact that

[TABLE]

is a real operator, the symmetry implies it being Hermitian, so that

[TABLE]

is real and nonzero. By leveraging a straightforward extension of Theorem 3.1 in [53], it follows that

[TABLE]

is always invertible. Otherwise said, the Yukawa-Calderón \acCFIE we propose is resonance free.

Acknowledgment

This work has been funded in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC project 321, grant No. 724846).