# Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning

Simon Adrian, Alexandre Dely, Davide Consoli, Adrien Merlini, Francesco Andriulli

► To cite this version:

Simon Adrian, Alexandre Dely, Davide Consoli, Adrien Merlini, Francesco Andriulli. Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning. IEEE Open Journal of Antennas and Propagation, 2021, 2, pp.1143-1174. 10.1109/OJAP.2021.3121097. hal-03472123

HAL Id: hal-03472123

https://hal.science/hal-03472123

Submitted on 10 Dec 2021

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>

Distributed under a Creative Commons Attribution 4.0 International License
Received 2 July 2021; revised 7 October 2021; accepted 8 October 2021. Date of publication 19 October 2021; date of current version 6 December 2021.

Digital Object Identifier 10.1109/OJAP.2021.3121097

# Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning

SIMON B. ADRIAN<sup>1</sup> (Member, IEEE), ALEXANDRE DÉLY<sup>2</sup>,

DAVIDE CONSOLI<sup>2</sup> (Graduate Student Member, IEEE), ADRIEN MERLINI<sup>4</sup> (Member, IEEE),

AND FRANCESCO P. ANDRIULLI<sup>2</sup> (Senior Member, IEEE)

CORRESPONDING AUTHOR: F. P. ANDRIULLI (e-mail: francesco.andriulli@polito.it)

This work was supported in part by the European Research Council (ERC) through the European Union's Horizon 2020 Research and Innovation Programme under Grant 724846 (Project 321); in part by the Italian Ministry of University and Research within the Program PRIN2017, EMVISIONING, under Grant 2017HZJXSZ, CUP:E64190025300; in part by the Italian Ministry of University and Research within the Program FARE, CELER, under Grant R187PMFXA4; in part by the Région Bretagne and the Conseil Départemental du Finistère through the Project "TONNERRE": in part by the ANR Labex CominLabs through the Project "CYCLE"; and in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grant SFB 1270/1-299150580.

**ABSTRACT** Integral equation formulations are a competitive strategy in computational electromagnetics but, lamentably, are often plagued by ill-conditioning and by related numerical instabilities that can jeopardize their effectiveness in several real case scenarios. Luckily, however, it is possible to leverage effective preconditioning and regularization strategies that can cure a large majority of these problems. Not surprisingly, integral equation preconditioning is currently a quite active field of research. To give the reader a propositive overview of the state of the art, this paper will review and discuss the main advancements in the field of integral equation preconditioning in electromagnetics summarizing the strengths and weaknesses of each technique. The contribution will guide the reader through the choices of the right preconditioner for a given application scenario. This will be complemented by new analyses and discussions which will provide a further and more intuitive understanding of the ill-conditioning of the electric field (EFIE), magnetic field (MFIE), and combined field integral equation (CFIE) and of the associated remedies.

**INDEX TERMS** Integral equations, boundary element method, computational electromagnetic, preconditioning, EFIE, MFIE.

## I. INTRODUCTION

INTEGRAL equation formulations, solved by the boundary element method (BEM), have become a well established tool to solve scattering and radiation problems in electromagnetics [1]–[4]. What makes these schemes so suitable for electromagnetic analyses is that, differently from approaches based on differential equations such as the finite element method (FEM) or the finite-difference time-domain method (FDTD), they naturally incorporate radiation conditions without the need for artificial absorbing boundary conditions, they only set unknowns on boundary surfaces (two-dimensional manifolds) instead of discretizing

the entire volume, and they are mostly free from numerical dispersion. On the other hand, linear system matrices arising from differential equations schemes are sparse [5], while those arising in BEM are, in general, dense. This drawback, however, can be overcome if a fast method such as the multilevel fast multipole method (MLFMM) [6], the multilevel matrix decomposition algorithm (MLMDA) [7] and later equivalents [8]–[12] are used at high frequency, or the adaptive cross approximation (ACA)/*H*-matrix methods and related schemes [13]–[18] are used at lower frequencies. These schemes are often capable of performing matrix-vector products in $O(N \log N)$ or even $O(N)$ complexity, where N

<footnote n="1">¹Fakultät für Informatik und Elektrotechnik, Universität Rostock, 18051 Rostock, Germany</footnote>

<footnote n="2">²Department of Electronics and Telecommunications (DET), Politecnico di Torino, 10129 Turin, Italy</footnote>

<footnote n="3">³Microwave Department, IMT Atlantique, 29285 Brest, France</footnote>
denotes the number of unknowns (the linear system matrix dimension). Thus the complexity to obtain the BEM solution of the electromagnetic problem is, when an iterative solver is used, $O(N_{iter}N \log N)$ (or $O(N_{iter}N)$ in the low-frequency regime), where $N_{iter}$ is the number of iterations.

The number of iterations $N_{iter}$ is generally correlated with the condition number of the linear system matrix, that is, the ratio between the largest and smallest singular values of the matrix [19]. This number is often a function of $N$ and, when the BEM formulation is set in the frequency domain, of the wavenumber $k$. This can potentially result in a solution complexity greater, and sometimes much greater, than $O(N^2)$, something that would severely jeopardize the other advantages of using BEM approaches.

For this reason it is of paramount importance to address and solve all sources of ill-conditioning for integral equations and, not surprisingly, this has been the target of substantial research in the last decade that this work will analyze, review, and summarize.

For surface integral equations (SIEs) that model scattering or radiation problems for perfect electrical conductors (PEC) geometries, we can typically distinguish the following sources of ill-conditioning: i) the low-frequency breakdown, ii) the h-refinement (dense-discretization) breakdown, iii) high-frequency issues (including internal resonances and the high-frequency breakdown), and iv) the lack of linear independence in the basis elements (including lack of orthogonality and mesh irregularities).

Some of the first methods explicitly addressing electromagnetic integral equation ill-conditioning date back to the 1980s, when the focus was on the low-frequency breakdown [20] and on the problem of interior resonances [21]. Since then, a plethora of schemes and strategies addressing one or more of the issues i)–iv) have been presented and some of these strategies are still the topic of intense research. In the past, a few review articles have appeared that dealt with aspects of stabilizing ill-conditioned electromagnetic integral equations. Most recently, Antoine and Darbas [22] presented an extensive review on operator preconditioning with a focus on high-frequency issues. A few years ago, Ylä-Oijala et al. [23] discussed issues in finding a stable and accurate integral equation formulation and they addressed certain open issues in preconditioning, and Carpentieri discussed preconditioning strategies with a focus on large-scale problems [24], [25].

Finally, although for space limitation this paper will focus on the electric field integral equation (EFIE) and magnetic field integral equation (MFIE) operators (which are the fundamental building blocks for several other formulations), the reader should note that a substantial amount of literature and quite effective preconditioned methods have been presented for modeling penetrable bodies both homogeneous and inhomogeneous [26]–[29]. The reader should also be aware that domain decomposition schemes can play a fundamental role in managing and solving electromagnetic problems containing even severely ill-conditioned operators [30], [31].

These approaches, however, are per se a discipline within Computational Electromagnetics and any brief treatment outside of a dedicated review would inevitably be insufficient and partial. Moreover, domain decomposition algorithms are not competing with the strategies discussed here but, often times, complementary [32]. For these reasons, we will not treat domain decomposition strategies in this review, but rather refer the interested reader to the excellent contributions in [33]–[35]. Similarly, discontinuous Galerkin and related methods for handling non-conformal meshes will not be treated here as extensive additional treatments would be required; the reader can refer to [36], [37], and references therein for specific discussions on this family of methods.

The purpose of this article is two-fold: on the one hand, we review and discuss the strategies that have been devised in the past to overcome the sources of ill-conditioning i)-iv) summarizing strengths and weaknesses, guiding the reader through the choices of the right preconditioner for a given application scenario. On the other hand, we complement the overview with new results that contribute to better characterizing the ill-conditioning of the EFIE and MFIE. Finally, we will complement our discussions with a spectral analysis of the formulations on the sphere, which will provide a further and more intuitive understanding of the ill-conditioning of the EFIE, MFIE, and combined field integral equation (CFIE) and of the associated potential remedies. In contrast to [22], our focus will include low-frequency effects and wideband stable formulations as well as Calderón and quasi-Helmholtz projection strategies. Moreover, whenever appropriate, we will provide implementational considerations and details that will enable the reader to dodge all practical challenges that are usually faced when engineering the most effective preconditioning schemes.

This paper is organized as follows: Section II introduces the background material and sets up the notation, Section III reviews the connection of the spectrum of matrices and the role of the condition number in the solutions of the associated linear systems. Section IV focuses on low-frequency scenarios analyzing their main challenges and solution strategies. Section V presents problems and solutions associated with highly refined meshes, while Section VI focuses on scenarios in the high-frequency regime. Section VII considers the role of mesh and basis functions quality on the overall conditioning and Section VIII presents the conclusions and final considerations.

## II. NOTATION AND BACKGROUND

We are interested in solving the electromagnetic scattering problem where a time-harmonic, electromagnetic wave ($e^i, h^i$) in a space with permittivity $\epsilon$ and permeability $\mu$ impinges on a connected domain $\Omega_-$ $\subset \mathbb{R}^3$ with PEC boundary $\Gamma := \partial\Omega_-$ resulting in the scattered wave ($e^s, h^s$). The total electric $e := e^i + e^s$ and magnetic $h := h^i + h^s$ fields satisfy Maxwell's equations

$$ \nabla \times e(\mathbf{r}) = +i k \mathbf{h}(\mathbf{r}) \quad \text{for all } \mathbf{r} \in \Omega_+, \quad (1) $$

$$ \nabla \times \mathbf{h}(\mathbf{r}) = -i k e(\mathbf{r}) \quad \text{for all } \mathbf{r} \in \Omega_+, \quad (2) $$
where $\Omega_+ := \overline{\Omega}_-^c$, $k := \omega\sqrt{\epsilon\mu}$ is the wave number, $\omega$ the angular frequency, and $e, h$ must satisfy the boundary conditions for PEC boundaries

$$ \hat{\mathbf{n}} \times \mathbf{e} = \mathbf{0} \quad \text{for all } \mathbf{r} \in \Gamma, \qquad (3) $$

$$ \hat{\mathbf{n}} \times \mathbf{h} = j_{\Gamma} \quad \text{for all } \mathbf{r} \in \Gamma, \qquad (4) $$

where $j_Γ$ is the induced electric surface current density. In addition, $e^s$ and $h^s$ must satisfy the Silver-Müller radiation condition [38], [39]

$$ \lim_{r \to \infty} |\mathbf{h}^s \times \mathbf{r} - r\mathbf{e}^s| = 0. \qquad (5) $$

We assumed (and suppressed) a time dependency of $e^{-i\omega t}$ and normalized $\mathbf{h}$ with the wave impedance $\eta := \sqrt{\mu/\epsilon}$.

To find $(e^s, h^s)$, we can solve the EFIE

$$ \mathcal{T}_k j_\Gamma = -\hat{\mathbf{n}} \times \mathbf{e}^i \qquad (6) $$

for $j_Γ$, where $\hat{\mathbf{n}}$ is the surface normal vector directed into $\Omega_+$ and

$$ \mathcal{T}_k := i k \mathcal{T}_{A,k} + 1/(i k) \mathcal{T}_{\Phi,k} \qquad (7) $$

is the EFIE operator composed of the vector potential operator

$$ (\mathcal{T}_{A,k} j_{\Gamma})(\mathbf{r}) = \hat{\mathbf{n}} \times \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}') j_{\Gamma}(\mathbf{r}') \, dS(\mathbf{r}') \qquad (8) $$

and the scalar potential operator

$$ (\mathcal{T}_{\Phi,k} j_{\Gamma})(\mathbf{r}) = -\hat{\mathbf{n}} \times \operatorname{grad}_{\Gamma} \times \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}') \operatorname{div}_{\Gamma} j_{\Gamma}(\mathbf{r}') \, dS(\mathbf{r}') \qquad (9) $$

where

$$ G_k(\mathbf{r}, \mathbf{r}') = \frac{e^{i k |\mathbf{r} - \mathbf{r}'|}}{4\pi |\mathbf{r} - \mathbf{r}'|} \qquad (10) $$

is the free-space Green's function. A definition of the surface differential operators $\operatorname{grad}_\Gamma$ and $\operatorname{div}_\Gamma$ can be found in [40, Appendix 3] or [41, Ch. 2]. Once $j_Γ$ is obtained, $e^s, h^s$ can be computed using the free-space radiation operators.

Alternatively, one can solve the MFIE for the exterior scattering problem

$$ \hat{\mathbf{n}} \times \mathbf{h}^i = \mathcal{M}_k^+ j_\Gamma := + (\mathcal{I}/2 + \mathcal{K}_k) j_\Gamma, \qquad (11) $$

where $\mathcal{I}$ is identity operator, $\mathcal{M}_k^+$ is the MFIE operator for the exterior scattering problem, and

$$ (\mathcal{K}_k j_\Gamma)(\mathbf{r}) := -\hat{\mathbf{n}} \times \int_{\Gamma} \nabla G_k(\mathbf{r}, \mathbf{r}') \times j_\Gamma \, dS(\mathbf{r}'). \qquad (12) $$

The MFIE operator for the interior scattering problem is $\mathcal{M}_k^- := -\mathcal{I}/2 + \mathcal{K}_k$ and will be used later in the construction of preconditioners.

The EFIE and the MFIE have non-unique solutions for resonance frequencies. A classical remedy is the use of the CFIE [21]

$$ \begin{aligned} -\alpha \mathcal{T}_k j_\Gamma + (1-\alpha) \hat{\mathbf{n}} \times \mathcal{M}_k^+ j_\Gamma &= \alpha \hat{\mathbf{n}} \times \mathbf{e}^i + (1-\alpha) \hat{\mathbf{n}} \times \hat{\mathbf{n}} \times \mathbf{h}^i \end{aligned} \qquad (13) $$

which is uniquely solvable for all frequencies.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. The vector field of an RWG function. The vector $e_n$ denotes the directed edge, $c_n^+$ and $c_n^-$ denote the domains of the cells, $v_n^+$ and $v_n^-$ denote vertices on the edge $e_n$, and $r_n^+$ and $r_n^-$ are the vertices opposite to the edge $e_n$.</figcaption></figure>

For the discretization of the EFIE, we employ Rao-Wilton-Glisson (RWG) basis functions $f_n \in X_f$ which are here—in contrast to their original definition in [42]—not normalized with the edge length, that is,

$$ f_n = \begin{cases} \frac{\mathbf{r} - \mathbf{r}_n^+}{2A_{c_n^+}} & \text{for } \mathbf{r} \in c_n^+, \\ \frac{\mathbf{r}_n^- - \mathbf{r}}{2A_{c_n^-}} & \text{for } \mathbf{r} \in c_n^- \end{cases} \qquad (14) $$

using the convention depicted in Figure 1.

Following a Petrov-Galerkin approach, we obtain the system of equations

$$ \mathcal{T}_k j = (i k \mathcal{T}_{A,k} + 1/(i k) \mathcal{T}_{\Phi,k}) j = -e^i \qquad (15) $$

that can be solved to obtain an approximation of the solution in the form $j_\Gamma \approx \sum_n j_n f_n$ and where

$$ [\mathcal{T}_{A,k}]_{nm} := \langle \hat{\mathbf{n}} \times f_n, \mathcal{T}_{A,k} f_m \rangle_{\Gamma}, \qquad (16) $$

$$ [\mathcal{T}_{\Phi,k}]_{nm} := \langle \hat{\mathbf{n}} \times f_n, \mathcal{T}_{\Phi,k} f_m \rangle_{\Gamma}, \qquad (17) $$

$$ \mathbf{e}_n^i := \langle \hat{\mathbf{n}} \times f_n, \hat{\mathbf{n}} \times \mathbf{e}^i \rangle_{\Gamma}, \qquad (18) $$

with

$$ \langle \mathbf{f}, \mathbf{g} \rangle_{\Gamma} := \int_{\Gamma} f(\mathbf{r}) \cdot \mathbf{g}(\mathbf{r}) \, dS(\mathbf{r}). \qquad (19) $$

Even though we are testing with $\hat{\mathbf{n}} \times f_n$, the resulting system matrix $\mathcal{T}_k$ is the one from [42] (up to the fact that the RWG functions we are using are not normalized), because our definition of the EFIE operator includes an $\hat{\mathbf{n}} \times$ term (in contrast to [42]).

For the discretization of the MFIE, functions dual to the RWGs must be used for testing [43]. Historically, the first dual basis functions for surface currents where introduced by Chen and Wilton for a discretization of the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) equation [44]. Later and independently, Buffa and Christiansen introduced the
<figure><img src="image_3.png" /><figcaption>FIGURE 2. The vector field of a BC function.</figcaption></figure>

Buffa-Christiansen (BC) functions [45], which differ from the Chen-Wilton (CW) functions in that the charge on the dual cells is not constant. Figure 2 shows a visualization of a BC function. In our implementation, we are using BC functions and denote them as $\tilde{f} \in X_f$, where the tilde indicates that the function is defined on the dual mesh. The analysis is, however, applicable to CW functions as well, and thus, we will mostly speak of “dual functions” to stress the generality of our analysis. For a definition of the BC functions as well as implementation details, we refer the reader to [46]. For the discretization of the MFIE, we obtain

$$M_k^+ j := (1/2 G_{\tilde{n} \times \tilde{f} f} + K_k) j = h^i, \quad (20)$$

where

$$[K_k]_{nm} := \langle \hat{n} \times \tilde{f}_n, K_k f_m \rangle_{\Gamma}, \quad (21)$$

$$[h^i]_n := \langle \hat{n} \times \tilde{f}_n, \hat{n} \times h^i \rangle_{\Gamma}, \quad (22)$$

and where the Gram matrix for any two function spaces $X_f$ and $X_g$ is defined as

$$[G_{f,g}]_{mn} := \langle f_m, g_n \rangle_{\Gamma}, \quad (23)$$

with $f_m \in X_f$ and $g_n \in X_g$.

For the discretization of the CFIE, we have

$$\begin{aligned} C_k j &:= (-\alpha T_k + (1-\alpha) G_{f,f} G_{\tilde{n} \times \tilde{f} f}^{-1} M_k^+) j \\ &= \alpha e^i + (1-\alpha) G_{f,f} G_{\tilde{n} \times \tilde{f} f}^{-1} h^i \end{aligned} \quad (24)$$

with the combination parameter $0 < \alpha < 1$.

### III. CONDITION NUMBERS, ITERATIVE SOLVERS, AND COMPUTATIONAL COMPLEXITY

To solve the linear system of equations arising from boundary element discretizations, such as (15), one can resort either to (fast) direct or to iterative solvers. For direct solvers, the time to obtain a solution is independent from the right-hand side, whereas for iterative solvers, the right-hand side as well as the spectral properties of the system matrix influence

the solution time. Standard direct solvers such as Gaussian elimination have a cubic complexity, which renders them unattractive for large linear systems. Recent progress in the development of fast direct solvers has improved the overall computational cost [47]–[50].

Iterative solvers, on the other hand, start from an initial guess of the solution, $x_{(0)}$, and compute a sequence of approximate solutions, where the following element of such a sequence is based on the previously computed one, until a desired accuracy is achieved. Formally, given a linear system of equations

$$Ax = b, \quad (25)$$

an iterative solver should stop when $\|Ax_{(i)} - b\|/\|b\| < \varepsilon$, where $\varepsilon > 0$ is the solver tolerance and $x_{(i)}$ the approximate solution after the $i$th iteration. Whether an iterative solver will converge or not, depends on the chosen solver and the properties of $A$, as we will discuss in the following.

To assess the overall complexity in $N$ for obtaining an approximation of $x$ within the tolerance $\varepsilon$, a relation between $N$-iter and $N$ is needed. One way to obtain such a relationship is via the condition number of the matrix, which is defined as

$$\mathrm{cond} \, A = \|A\|_2 \|A^{-1}\|_2 = \frac{\Sigma_{\max}(A)}{\Sigma_{\min}(A)}, \quad (26)$$

where $\|\cdot\|_2$ is the spectral norm, and $\Sigma_{\max/\min}$ denotes the maximal and minimal singular value.

In the case of the conjugate gradient (CG) method, which requires $A$ to be Hermitian and positive definite, there is an upper bound on the error $e_{(i)} := x_{(i)} - x$ given by [51]

$$\|e_{(i)}\|_{A} \le 2 \left( \frac{\sqrt{\mathrm{cond} \, A} - 1}{\sqrt{\mathrm{cond} \, A} + 1} \right)^{i} \|e_{(0)}\|_{A}, \quad (27)$$

where $\|\cdot\|_A$ is the energy norm defined by $\|x\|_A := (x^\dagger Ax)^{1/2}$ and $x^\dagger$ denotes the conjugate transpose of $x$. If the objective is to reduce the relative error $\|e_{(i)}\|/\|e_{(0)}\|$ below $\varepsilon$ and by considering limits for $\mathrm{cond} \, A \gg 1$, one notes [51] that

$$i \le \left[ \frac{1}{2} \sqrt{\mathrm{cond} \, A \ln\left(\frac{2}{\varepsilon}\right)} \right] \quad (28)$$

iterations are at most needed (assuming an exact arithmetic). If the condition number grows linearly in $N$, as observed for the EFIE when the mesh is uniformly refined, this implies that the complexity is at most $O(N^{1.5} \log N)$.

One could argue that this is an overly simplified picture of the situation; indeed, the CG method is not applicable to standard frequency domain integral equations as the resulting system matrices are neither Hermitian nor positive definite. One strategy to still obtain a bound on the number of iterations is to use the CG method on the normal equation

$$A^\dagger Ax = A^\dagger b. \quad (29)$$

The price for this, however, is that the condition number of the resulting system matrix is $(\mathrm{cond} \, A)^2$ and thus this
approach is, for the standard formulations, of little practical value. In addition, round-off errors due to finite precision can lead to a non-converging solver—despite the theory dictating that CG should converge in at most N steps [52], [53]. Thus, the condition number bound is relevant in practice often only in the case that cond A is small.

The problem with other popular Krylov methods such as the generalized minimal residual (GMRES) or the conjugate gradient squared (CGS) method is that, for general matrices, no bound on the number of iterations in terms of the condition number alone is available. In fact, even if two matrices have the same condition number, their convergence behavior can significantly differ: the distribution of the eigenvalues in the complex plane impacts the convergence behavior as well [22]. Typically, a better convergence can be observed if all the eigenvalues are located on either the real or the imaginary axis and if they are either strictly positive or negative (if they are on the imaginary axis, then positive or negative with respect to Im($\Lambda_f$)). We will see in the following that, under certain conditions, for low-frequency electromagnetic problems it is possible to cluster the eigenvalues on the real axis and that the condition number becomes a good indicator of the convergence behavior. Moreover, some preconditioning strategies, such as the refinement-free Calderón preconditioner which will be discussed in Section V-A2, give rise to a Hermitian, positive-definite system, and thus the CG and the associated convergence theory is applicable.

For frequency-independent problems, it is customary to call a formulation well-conditioned if cond A is asymptotically bounded by a constant C which is independent from the average edge length h of the mesh. For dynamic problems, however, we also need to study the condition number as a function of the frequency $f := \omega/(2\pi)$, and one must specify if a formulation is well-conditioned with respect to h, to f, to both, or only in a particular regime, for example, for frequencies where the corresponding wavelength is larger than the diameter of $\Gamma$.

The classical remedy to overcome ill-conditioning and thus improve the convergence behavior of iterative solvers is to use a preconditioning strategy. Such a strategy results, in the general case, in a linear system

$$P_L A P_R y = P_L b, \quad (30)$$

where $x = P_R y$ and the matrices should be chosen such that, if possible,

$$\mathrm{cond}(P_L A P_R) \le C, \quad (31)$$

where C is a constant both independent of h and f (in which case the preconditioner is optimal). Normally, the matrix-vector products in (30) are not formed explicitly and, to be an efficient preconditioner, the cost of a matrix-vector product should not jeopardize the lead complexity set by the fast method. In practice, to obtain an optimal preconditioner, the nature of the underlying operators must be taken into account. Thus, in the following sections, we will analyze the spectral properties of the (discretized) EFIE, MFIE, and

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Radar cross sections calculated, with different formulations, for the sphere of unit radius discretized with an average edge length of 0.15 m, and excited by a plane wave of unit polarization along $\hat{x}$ and propagation along $\hat{z}$ oscillating at $f = 10^{-20}$ Hz. The “EFIE” and “MFIE” labels refer to the standard formulations (15) and (20), while the “Loop-star EFIE” and “P-EFIE” refer to the EFIE stabilized with the loop-star (61) and quasi-Helmholtz projectors (72), respectively.</figcaption></figure>

CFIE operator, discuss the causes of their ill-conditioning as well as potential remedies.

## IV. LOW-FREQUENCY SCENARIOS

The low-frequency breakdown of the EFIE, that is, the growth of the condition number when the frequency f decreases, was one of the first sources of ill-conditioning of the EFIE to be studied. From a physical point of view, several problems at low-frequency are rooted in the decoupling of the electric and the magnetic field in the static limit: magnetostatic loop currents excite the magnetic field and electrostatic charges excite the electric field [20]. Both the EFIE and the MFIE suffer from computational challenges at low-frequencies. As we will see in this section, the EFIE suffers from conditioning issues when the frequencies decreases and so does, albeit for different reasons, the MFIE when applied to non-simply connected geometries (i.e., geometries containing handles like the torus illustrated in Figure 4, for example). The condition number growth is, however, only one of the possible problems: finite machine precision and inaccuracies due to numerical integration that result in catastrophic round-off errors are also plaguing the otherwise low-frequency well-conditioned integral equations such as the MFIE on simply-connected geometries. Together, these issues make the two formulations increasingly inaccurate as the frequency decreases, which is attested by the low-frequency radar cross sections illustrated in Figure 3 that show wildly inaccurate results for the standard formulations.

The low-frequency analysis of electromagnetic integral equations benefits from the use of Helmholtz and quasi-Helmholtz decompositions that we will summarize here for the sake of completeness and understanding. The well-known Helmholtz decomposition theorem states that any vector field
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Illustration of a torus and the corresponding toroidal (in blue) and poloidal (in orange) loops.</figcaption></figure>

can be decomposed into a solenoidal, irrotational, and a harmonic vector field, which in the case of a tangential surface vector field such as $j_Γ$ leads to [41, p. 251]

$$j_Γ = \mathbf{curl}_Γ \Phi + \mathrm{grad}_Γ \Psi + \mathbf{H}, \quad (32)$$

where $\Phi$ and $\Psi$ are sufficiently smooth scalar functions, $\mathbf{curl}_Γ \Phi := \mathrm{grad}_Γ \Phi \times \hat{n}$, and $\mathrm{div}_Γ \mathbf{H} = \mathbf{curl}_Γ \mathbf{H} = 0$; here, $\mathbf{curl}_Γ$ is the adjoint operator of $\mathbf{curl}_Γ$, that is, we have $\langle \mathbf{curl}_Γ f, g \rangle_Γ = \langle f, \mathbf{curl}_Γ g \rangle_Γ$ (see [41, eq. (2.5.194)]). The space of harmonic functions $H_\mathbf{H}(\Gamma)$ is finite dimensional with dim $H_\mathbf{H}(\Gamma) = 2g$ on a closed surface, where $g$ is the genus of $\Gamma$. The Helmholtz subspaces are all mutually orthogonal with respect to the $L^2(\Gamma)$-inner product.

When $j_Γ$ is a linear combination of div- but not curl-conforming functions (e.g., RWG and BC functions), only a quasi-Helmholtz decomposition is possible, where $j_Γ$ is decomposed into a solenoidal, a non-solenoidal, and a quasi-harmonic current density. It is not possible to obtain irrotational or harmonic current densities, since the curl of div-conforming (but not curl-conforming) functions such as the RWGs (or their dual counterparts) is, in general, not existing as a classical derivative; therefore, it is termed quasi-Helmholtz decomposition. Next we introduce the quasi-Helmholtz decompositions for primal (i.e., RWGs) and dual (i.e., BCs) functions that we will use for our analysis in the next section.

Just as the Helmholtz decomposition (32) decomposes the continuous solution $j_Γ$, a quasi-Helmholtz decomposition decomposes the discrete solution $j$ as

$$\sum_{n=1}^{N} [j]_n f_n = \sum_{n=1}^{N_V} [j_\Lambda]_n \Lambda_n + \sum_{n=1}^{N_C} [j_\Sigma]_n \Sigma_n + \sum_{n=1}^{2g} [j_H]_n H_n, \quad (33)$$

where $\Lambda_n \in X_\Lambda$ are solenoidal loop functions, $\Sigma_n \in X_\Sigma$ are non-solenoidal star functions, and $H_n \in X_\mathbf{H}$ are quasi-harmonic global loops [54] and where $j_\Lambda$, $j_\Sigma$, and $j_\mathbf{H}$ are the vectors containing the associated expansion coefficients; moreover, $N_V$ is the number of vertices and $N_C$ is the number of cells of the mesh.

We highlight some of the properties which we are going to use throughout this article. First, and most importantly,

the functions $\Lambda_n$, $H_n$, and $\Sigma_n$ can be represented in terms of RWG functions [54]. Thus the expansion coefficients are linked by linear transformation matrices $\Lambda$, $H$, and $\Sigma$. For the loop transformation matrix, we have

$$[\Lambda]_{ij} = \begin{cases} 1 & \text{for } v_j = v_i^-, \\ -1 & \text{for } v_j = v_i^+, \\ 0 & \text{otherwise,} \end{cases} \quad (34)$$

where $v_j$ is the $j$th vertex of the mesh (inner vertex if $\Gamma$ is open), and for the star transformation matrix

$$[\Sigma]_{ij} = \begin{cases} 1 & \text{for } c_j = c_i^+, \\ -1 & \text{for } c_j = c_i^-, \\ 0 & \text{otherwise,} \end{cases} \quad (35)$$

where $c_j$ is the $j$th cell of the mesh, following the conventions depicted in Figure 1. With the definition of these matrices, the quasi-Helmholtz decomposition in (33) can be equivalently written as

$$j = \underbrace{\Lambda j_\Lambda}_{=j_{\mathrm{sol}}} + \underbrace{H j_H}_{=j_{\mathrm{qhar}}} + \underbrace{\Sigma j_\Sigma}_{=j_{\mathrm{nsol}}} = j_{\mathrm{sol}} + j_{\mathrm{qhar}} + j_{\mathrm{nsol}}. \quad (36)$$

The linear combinations of RWGs implied by the coefficient vectors $j_{\mathrm{sol}}$, $j_{\mathrm{nsol}}$, and $j_{\mathrm{qhar}}$ are solenoidal, non-solenoidal, and quasi-harmonic current densities. These decompositions are not unique: if we were to use, for example, the loop-tree quasi-Helmholtz decomposition, we would obtain different coefficient vectors $j_{\mathrm{sol}}$, $j_{\mathrm{nsol}}$, and $j_{\mathrm{qhar}}$. The decomposition is, however, unique with respect to the loop-star space, that is, when the linear dependency of loop and of star functions (see [55] and references therein) is not resolved by arbitrarily eliminating a loop and a star function; what sets the loop-star basis apart from other quasi-Helmholtz decompositions is the symmetry with respect to dual basis functions. A symmetry that we are now going to further highlight.

First, we give to $\Lambda$ and $\Sigma$ a meaning that goes beyond merely interpreting them as basis transformation matrices. The matrices $\Lambda$ and $\Sigma$ are edge-node and edge-cell incidence matrices of the graph defined by the mesh and they are orthogonal, that is, $\Sigma^T \Lambda = 0$. It follows that $j_{\mathrm{sol}}$ and $j_{\mathrm{nsol}}$ are $l^2$-orthogonal, that is, $j_{\mathrm{nsol}}^T j_{\mathrm{sol}} = 0$. We find this noteworthy for two reasons: i) the loop $\Lambda_i$ and star functions $\Sigma_j$ are, in general, not $L^2$-orthogonal (after all, $\Sigma_j$ is not irrotational); ii) the $l^2$-orthogonality is not true for other quasi-Helmholtz decompositions such as the loop-tree basis. In light of this consideration, the matrices $\Lambda$ and $\Sigma$ could be interpreted as the graph curl ($\Lambda$) and graph gradient ($\Sigma$) of the standard mesh, an interpretation that further increases the correspondence with the continuous decomposition (32).

For global loops $H_n$, no such simple graph-based definition exists. Indeed, they are, in general, not uniquely defined and must be constructed from a search of holes and handles. For any global loop basis so obtained, we have $\Sigma^T H = 0$; however, $\Lambda^T H = 0$ is, in general, not true. This property can be enforced by constructing $H$ as the right nullspace of $[\Lambda \ \Sigma]^T$. Such a construction is possible, for example, via a
full singular value decomposition (SVD), or, via more com-
putationally efficient randomized projections [56]. However,
the computational cost is higher, in general, compared with
using a global loop-finding algorithm, in particular, since H
will be a dense matrix.

A similar decomposition can be obtained for dual
functions

$$
\sum_{n=1}^{N} [m]_n \tilde{f}_n = \sum_{n=1}^{N_v} [m_{\tilde{\Lambda}}]_n \tilde{\Lambda}_n + \sum_{n=1}^{N_c} [m_{\tilde{\Sigma}}]_n \tilde{\Sigma}_n \\
+ \sum_{n=1}^{2g} [m_{\tilde{H}}]_n \tilde{H}_n, \qquad (37)
$$

where, in contrast to the RWG case, $\tilde{\Lambda}_n$ are non-solenoidal dual star and $\tilde{\Sigma}_n$ are solenoidal dual loop functions. In matrix notation, we have

$$
j = \underbrace{\Sigma m_{\tilde{\Sigma}}}_{m_{\text{sol}}} + \underbrace{\tilde{H} m_{\tilde{H}}}_{m_{\text{qhar}}} + \underbrace{\Lambda m_{\tilde{\Lambda}}}_{m_{\text{nsol}}} = m_{\text{sol}} + m_{\text{qhar}} + m_{\text{nsol}}. \quad (38)
$$

Note that the same matrices **Σ** and **Λ** are present both in the
decomposition of RWG functions and in the one of dual func-
tions. However, while for RWGs the transformation matrix **Λ**
describes solenoidal functions and the transformation matrix
**Σ** describes non-solenoidal functions, the opposite is true
for dual functions: it is **Σ** that describes solenoidal func-
tions, while **Λ** describes non-solenoidal. Thus on the dual
mesh, **Λ** acts as graph gradient and **Σ** as a graph curl. This
is consistent with the definition of dual functions: dual basis
functions can be interpreted as a div-conforming “rotation”
by 90° of the primal functions (note that the functions **n̂** × **f**ᵢ
are a rotation by 90°, which is not div-confirming); given
that **curl**ΓΦ := ∇ΓΦ × **n̂**, it is consistent that the roles
of **Λ** and **Σ** as graph counterparts to continuous differential
surface operators are swapped on the dual mesh with respect
to the primal mesh.

Regarding the quasi-harmonic functions, it must be
emphasized that we cannot identify $\tilde{H} = H$. This equality
is only true if $H$ is the nullspace of $[\Lambda \ \Sigma]^T$, a condi-
tion, which evidently leads to the aforementioned unique
definition of $H$. Even though the construction of $H$ as the
nullspace of $[\Lambda \ \Sigma]^T$ is cumbersome—and by introducing
quasi-Helmholtz projectors in the following, we will sidestep
it—it suggests that these global loops are capturing the ana-
lytic harmonic Helmholtz subspace better than arbitrarily
chosen global loops.

A. ELECTRIC FIELD INTEGRAL EQUATION

To put into light the low-frequency challenges that plague
the EFIE, its behavior on both the solenoidal and the
non-solenoidal subspaces must be analyzed. The following
developments focus on geometries that do not contain global
loops, however the results can be immediately extended
to the general case by considering that—in the case of
the EFIE—global and local loops have similar properties.
While they have practical limitations, loop-star bases are

a convenient tool to perform this analysis. The loop-star
transformed EFIE matrix $\mathcal{T}_k^{\text{LS}} := [\Lambda \ \Sigma]^T \mathcal{T}_k [\Lambda \ \Sigma]$ can
be represented in block matrix form as

$$
\mathcal{T}_k^{\text{LS}} = \begin{bmatrix} \Lambda^\text{T} \mathcal{T}_k \Lambda & \Lambda^\text{T} \mathcal{T}_k \Sigma \\ \Sigma^\text{T} \mathcal{T}_k \Lambda & \Sigma^\text{T} \mathcal{T}_k \Sigma \end{bmatrix}, \quad (39)
$$

and the corresponding matrix equation now reads $\mathcal{T}_k^{\text{LS}} j_{LS} = [\Lambda \ \Sigma]^T e^i$, where $j = [\Lambda \ \Sigma] j_{LS}$. In these definitions, the $\Lambda$ and $\Sigma$ matrices refer to the full-rank transformation matrices in which linearly dependent columns have been removed: for each connected component of $\Gamma$ one star basis function (column of $\Sigma$) must always be removed and one loop basis function must be removed (column of $\Lambda$) if the component is closed [57].

To evidence the different low-frequency behaviors of the
EFIE matrix on the solenoidal and non-solenoidal subspaces,
the properties $\Lambda^T \mathcal{T}_{\Phi,k} = 0$ and $\mathcal{T}_{\Phi,k} \Lambda = 0$, which follow
directly from the divergence-free nature of solenoidal func-
tions, must be enforced. In addition, the behavior of the
matrix terms must be derived by performing a Taylor series
expansion of the Green's function in both $\mathcal{T}_{A,k}$ and $\mathcal{T}_{\Phi,k}$ for
$k \to 0$. For instance,

$$
\langle \hat{n} \times \mathbf{\Sigma}_m, \mathcal{T}_{A,k} \mathbf{\Lambda}_n \rangle_{\Gamma} \underset{k \to 0}{=} \int_{\Gamma} \int_{\Gamma} \frac{\mathbf{\Sigma}_m(\mathbf{r}) \cdot \mathbf{\Lambda}_n(\mathbf{r}')}{4\pi} \left( \frac{1}{R} - \frac{k^2 R}{2} - \frac{i k^3 R^2}{6} + \mathcal{O}(k^4) \right) \mathrm{d}S(\mathbf{r}') \mathrm{d}S(\mathbf{r}), \quad (40)
$$

where $R = |\mathbf{r} - \mathbf{r}'|$ and where we have used
$\int_{\Gamma} i \mathbf{k} \mathbf{\Lambda}_n(\mathbf{r}') \mathrm{d}S(\mathbf{r}') = 0$. We can deduce that, in general,

$$
\mathrm{Re} \left\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_{A,k} \mathbf{\Lambda}_n \right\rangle_{\Gamma} \underset{k \to 0}{=} \mathcal{O}(1), \quad (41)
$$

$$
\mathrm{Im} \left\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_{A,k} \mathbf{\Lambda}_n \right\rangle_{\Gamma} \underset{k \to 0}{=} \mathcal{O}(k^3). \quad (42)
$$

This process can be repeated for both $\mathcal{T}_{A,k}$ and $\mathcal{T}_{\Phi,k}$ when both expansion and testing functions are non-solenoidal and when at least one of the two is solenoidal. In summary,

$$
\mathrm{Re} \left\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_k \mathbf{\Sigma}_n \right\rangle_{\Gamma} \underset{k \to 0}{=} \mathcal{O}(k^2), \quad (43)
$$

$$
\mathrm{Im} \left\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_k \mathbf{\Sigma}_n \right\rangle_{\Gamma} \underset{k \to 0}{=} \mathcal{O}(k^{-1}), \quad (44)
$$

$$
\mathrm{Re} \left\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_k \mathbf{\Lambda}_n \right\rangle_{\Gamma} \underset{k \to 0}{=} \mathcal{O}(k^4), \quad (45)
$$

$$
\mathrm{Im} \left\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_k \mathbf{\Lambda}_n \right\rangle_{\Gamma} \underset{k \to 0}{=} \mathcal{O}(k). \quad (46)
$$

By symmetry, both $\langle \hat{\mathbf{n}} \times \mathbf{\Lambda}_m, \mathcal{T}_k \mathbf{\Sigma}_n \rangle_{\Gamma}$ and $\langle \hat{\mathbf{n}} \times \mathbf{\Lambda}_m, \mathcal{T}_k \mathbf{\Lambda}_n \rangle_{\Gamma}$ have the same low-frequency behavior as $\langle \hat{\mathbf{n}} \times \mathbf{\Sigma}_m, \mathcal{T}_k \mathbf{\Lambda}_n \rangle_{\Gamma}$. The scaling of the behavior of the block matrix is now straightforward to obtain, resulting in

$$
\mathrm{Re} \left( \mathcal{T}_k^{\mathrm{LS}} \right) \underset{k \to 0}{=} \begin{bmatrix} O(k^4) & O(k^4) \\ O(k^4) & O(k^2) \end{bmatrix}, \quad (47)
$$

$$
\mathrm{Im} \left( \mathcal{T}_k^{\mathrm{LS}} \right) \underset{k \to 0}{=} \begin{bmatrix} O(k) & O(k) \\ O(k) & O(k^{-1}) \end{bmatrix}, \quad (48)
$$
and the dominant behavior of $\mathcal{T}_k^{LS}$ is that of its imaginary part.

These results can be used to demonstrate the issues plaguing the EFIE at low frequencies, starting with its ill-conditioning. Consider the block diagonal matrix $D_k = \text{diag}[k^{-1/2} k^{1/2}]$ in which the block dimensions are consistent with that of the loop star decomposition matrix. Clearly,

$$D_k \mathcal{T}_k^{LS} D_k = \begin{bmatrix} O(1) & O(k) \\ O(k) & O(1) \end{bmatrix}, \quad (49)$$

is a well-conditioned matrix, in the sense that $\lim_{k \to 0} \text{cond} \, D_k \mathcal{T}_k^{LS} D_k := \gamma$ is finite. It then follows that

$$\begin{aligned} \text{cond} \, \mathcal{T}_k^{LS} &= \text{cond} \, D_k^{-1} D_k \mathcal{T}_k^{LS} D_k D_k^{-1} \\ &\le (\text{cond} \, D_k)^2 \text{cond} \, D_k \mathcal{T}_k^{LS} D_k \end{aligned} \quad (50)$$

and thus $\lim_{k \to 0} \text{cond} \, \mathcal{T}_k^{LS} = O(k^{-2})$. A lower bound for the condition number of interest can be obtained through the application of the Gershgorin disk theorem after diagonalization of the bottom right block of $\mathcal{T}_k^{LS}$, which proves that $\text{cond}(\mathcal{T}_k^{LS}) \ge \Sigma_{\min} k^{-2}$ where $\Sigma_{\min}$ is the smallest singular value of $(\Sigma^T \mathcal{T}_{\Phi,0} \Sigma)$. Considering these results and that the loop-star transformation matrix is invertible and frequency independent, we conclude that $\text{cond}(\mathcal{T}_k) \sim k^{-2}$ when $k \to 0$.

The second source of instability of the EFIE at low frequencies is the loss of significant digits in the right-hand side $e^i$, solution $j$, or radiated fields. To see this effect, the behavior of the right-hand side of the EFIE must be considered. Here we will restrict our developments to the plane-wave excitation, but similar results can be obtained for other problems [58]. Following the same procedure as for the matrix elements, we can determine the behavior of the loop and star right-hand side elements

$$\text{Re} \left\langle \hat{n} \times \Lambda_m, \hat{n} \times e_{\text{PW}}^i \right\rangle_{\Gamma} = O(k^2), \quad (51)$$

$$\text{Im} \left\langle \hat{n} \times \Lambda_m, \hat{n} \times e_{\text{PW}}^i \right\rangle_{\Gamma} = O(k), \quad (52)$$

$$\text{Re} \left\langle \hat{n} \times \Sigma_m, \hat{n} \times e_{\text{PW}}^i \right\rangle_{\Gamma} = O(1), \quad (53)$$

$$\text{Im} \left\langle \hat{n} \times \Sigma_m, \hat{n} \times e_{\text{PW}}^i \right\rangle_{\Gamma} = O(k), \quad (54)$$

where $e_{\text{PW}}^i$ is the electric field of the incident plane-wave. It is crucial to remember that when the standard EFIE—with no treatment—is solved numerically in finite precision floating point arithmetic, the real parts (resp. imaginary parts) of the loop and star components of the right-hand side are stored in the same floating point number. In particular, the real part of the solenoidal component that behaves as $O(k^2)$ is summed with an asymptotically much larger non-solenoidal component behaving as $O(1)$. In the context of finite precision arithmetic, the dynamic range of the floating point number will be imposed by the larger of the two components, meaning that the floating point number will become increasingly incapable of storing accurately the smaller one. This loss of

significant digits will worsen until the solenoidal component has completely vanished from the numerical value. This phenomenon is not necessarily damageable *per se*, but can lead to drastic losses in solution accuracy. In the particular case of the plane-wave excitation, we will study the effect of this loss of accuracy on the dominant parts of the solution. Using the well-known relations on block matrix inverses [59], one can show that

$$\text{Re} \left( \mathcal{T}_k^{LS} \right)^{-1} = \begin{bmatrix} O(k^2) & O(k^4) \\ O(k^4) & O(k^4) \end{bmatrix}, \quad (55)$$

$$\text{Im} \left( \mathcal{T}_k^{LS} \right)^{-1} = \begin{bmatrix} O(k^{-1}) & O(k) \\ O(k) & O(k) \end{bmatrix}, \quad (56)$$

which, in combination with the right-hand side results yields the behavior of the solution coefficients

$$\text{Re} \left( j_{\Lambda} \right) = O(1), \quad (57)$$

$$\text{Im} \left( j_{\Lambda} \right) = O(k), \quad (58)$$

$$\text{Re} \left( j_{\Sigma} \right) = O(k^2), \quad (59)$$

$$\text{Im} \left( j_{\Sigma} \right) = O(k), \quad (60)$$

which are indeed the behavior predicted by physics [60]. Note that the inaccurate right-hand side component will only have a significant contribution to the imaginary part of the solenoidal component of the solution, which is non-dominant. As such, although the error of the current could be low, the error of the charge or field could be quite high.

Finally, the reader should note that to numerically observe these results—and successfully implement the remedies that we will see later on—the vanishing of all relevant integrals must be explicitly enforced in some way, because floating point arithmetic and numerical integration are not capable of obtaining an exact zero in their computation and will saturate at machine precision, in the best case scenarios. Indeed, had they not been enforced, the solenoidal and non-solenoidal parts of the solution would have had the same behavior and, as such, would not yield a solution behaving as predicted by physics.

## 1) LOOP-STAR/TREE APPROACHES

Historically, the loop-star and loop-tree decompositions have been used to cure the low-frequency breakdown of the EFIE [20], [54] and as such are well-known and studied [55]. The fundamental curing mechanism of these approaches is to decompose the EFIE system using a RWG-to-loop-star or RWG-to-loop-tree mapping and isolate the solenoidal and non-solenoidal parts of the system. This separation allows for a diagonal preconditioning of the decomposed matrix to cure its ill-conditioning (as was done in Section IV-A). In addition, this separation makes it possible to enforce that the required integrals and matrix products vanish and cures the loss of significant digits that plagues the EFIE, since the loop and star contributions of each entity are stored in separate floating point numbers. In the case of the loop-star
approach, the stabilized matrix system is

$$D_k T_k^{LS} D_k j_{DLS} = D_k [\Lambda \ \Sigma]^T e^i, \quad (61)$$

where $j = [\Lambda \ \Sigma]D_k j_{DLS}$, following the notations of Section IV-A. Once the intermediate solution $j_{DLS}$ has been obtained, it must be handled with particular care. If, for instance, the quantity of interest is the field radiated by the solution, the radiation operators must be applied separately on the solenoidal and non-solenoidal parts of the solution that can be retrieved as $j_{sol} = [\Lambda \ 0]D_k j_{DLS}$ and $j_{sol} = [0 \ \Sigma]D_k j_{DLS}$, because additional vanishing integrals must be enforced in the scattering operators when applied to solenoidal functions. In addition, any explicit computation of $j$ would be subject to a numerical loss of significance and would further compromise the accuracy of the fields.

The key difference between loop-tree and loop-star techniques is that, in the former, the quasi-Helmholtz decomposition leverages a tree basis in place of the star basis, as indicated by their names. To define this tree basis consider the connectivity graph joining the centroids of all adjacent triangle cells of the mesh. To each edge of this graph corresponds a unique RWG function. Then, given a spanning tree of this graph, a tree basis can be defined as the subset $\{\theta_j\}$ of the RWG functions whose corresponding edge in the connectivity graph is included the spanning tree [54], [61]. The rationale behind the technique is that, by construction, such a basis will not be capable of representing any loop function. Clearly, the construction of this basis is not unique, since it depends on the choice of spanning tree. In practice, the loop-tree approach results in a matrix system similar to (61), in which the RWG-to-loop-star mapping $[\Lambda \ \Sigma]$ is replaced by an RWG-to-loop-tree mapping $[\Lambda \ \Theta]$ and $T_k^{LS}$ becomes $T_k^{LT} := [\Lambda \ \Theta]^T T_k [\Lambda \ \Theta]$ where

$$[\Theta]_{ij} = \begin{cases} 1 & \text{if } f_i = \theta_j, \\ 0 & \text{otherwise,} \end{cases} \quad (62)$$

is the general term of the RWG-to-tree transformation matrix. The resulting preconditioned equation is

$$D_k T_k^{LT} D_k j_{DLT} = D_k [\Lambda \ \Theta]^T e^i, \quad (63)$$

where $j = [\Lambda \ \Theta]D_k j_{DLT}$.

At first glance, the computational overhead of the two methods seems low, since $\Lambda$, $\Sigma$, $\Theta$, and $D_k$ are sparse matrices. However, while both methods adequately address the low-frequency breakdown of the EFIE, in the sense that they yield the correct solution (Figure 3) and prevent the conditioning of the system to grow unbounded as the frequency decreases (Figure 5), they cause the conditioning of the system matrix to artificially worsen because the loop-star and loop-tree bases are ill-conditioned [62]. This has led to the development of a permutated loop-star and loop-tree bases to reduce the number of iterations required to solve the preconditioned system using iterative solvers [61]. In general, the loop-tree preconditioned EFIE was observed to converge

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Comparison of the conditioning of the system matrices for several formulations on a sphere of radius 1 m discretized with an average edge length of 0.15 m, for varying frequency.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Comparison of the conditioning of the loop-star, loop-tree, and projector-based preconditioned EFIE matrices on a spheres of radius 1 m discretized with an average edge length of 0.3 m (solid lines) and 0.2 m (dotted lines) as a function of the frequency. The labels “Loop-tree EFIE” and “P-CMP-EFIE” refer to the EFIE stabilized with the loop-tree approach (63) and the Calderón EFIE stabilized with quasi-Helmholtz projectors (128).</figcaption></figure>

faster than the loop-star preconditioned EFIE [63], which can be explained by the fact that $\Lambda$ and $\Sigma$ can be interpreted as the discretizations of the graph curl and graph gradient [55], [62], that are ill-conditioned derivative operators. While a rigorous proof of the effect of this ill-conditioning on the preconditioned EFIE matrix is out of the scope of this review, pseudo-differential operator theory can be used to show that the differential strength of the loop-star transformation operators is sufficiently high not to be compensated by that of the vector potential. To illustrate this adverse effect, the conditioning of the system matrices has been obtained numerically and is presented in Figure 6. Clearly, the standard EFIE matrix shows a condition number growing as the frequency
decreases. However, at moderate frequencies, the condition-
ing of the loop-star and loop-tree preconditioned matrices is
significantly higher than that of the original matrix.

2) QUASI-HELMHOLTZ PROJECTORS

From the previous sections, it is clear that although the loop-
star/tree decompositions are helpful in analyzing the reasons
behind the low-frequency breakdown and that historically
provided a cure for it, they still give rise to high condition
numbers since they introduce an ill-conditioning related to
the mesh discretization. Moreover, for non-simply connected
geometries, loop-star decompositions require a search for the
mesh global cycles, an operation that can be computationally
cumbersome.

A family of strategies to overcome the drawbacks of loop-star/tree decompositions while still curing the low-frequency breakdown is the one based on quasi-Helmholtz projectors [62], [64]. Quasi-Helmholtz projectors can decompose the current and the operators into solenoidal and non-solenoidal components (just like a loop-star/tree decomposition does) but, being projectors, have a flat spectrum that, differently from loop-star/tree decompositions, do not alter the spectral slopes of the original operators and thus do not introduce further ill-conditioning.

Starting from the quasi-Helmholtz decomposition (36)

$$
j = \Sigma j_{\Sigma} + \Lambda j_{\Lambda} + H j_{H}, \quad (64)
$$

the quasi-Helmholtz projector for the non-solenoidal part is
the operator that maps j into Σj<sub>Σ</sub>. Since

$$
\Sigma^T j = \Sigma^T \Sigma j_{\Sigma}, \quad (65)
$$

the looked for projector is

$$
P_{\Sigma} := \Sigma (\Sigma^{\mathrm{T}} \Sigma)^{+} \Sigma^{\mathrm{T}}, \quad (66)
$$

where + denotes the Moore–Penrose pseudoinverse. The pro-
jector for the solenoidal plus harmonic components can be
obtained out of complementarity as

$$
P_{\Lambda H} := \mathbf{I} - P_{\Sigma}. \qquad (67)
$$

The same reasoning for dual functions leads to the dual
definitions of the projector

$$
P_{\Lambda} := \Lambda(\Lambda^{\mathrm{T}}\Lambda)^{+} + \Lambda^{\mathrm{T}}, \quad (68)
$$

which is the non-solenoidal projector for dual functions.
The solenoidal plus harmonic projector for dual functions
is, again, obtained by complementarity as

$$
P_{\Sigma H} := \mathbf{I} - P_{\Lambda}. \qquad (69)
$$

It is important to note that, even though the projectors
presented so far include a pseudo-inverse in their defini-
tion, they can be applied to arbitrary vectors in quasi-linear
complexity by leveraging algebraic multigrid precondition-
ing [62], [65], [66] and, as such, are fully compatible with
standard fast solvers.

Quasi-Helmholtz projectors can be used to cure the dif-
ferent deleterious effects of the low-frequency breakdown
by isolating the solenoidal and non-solenoidal parts of the
system matrix, unknowns, and right-hand side and rescaling
them appropriately. Thus they are an alternative to loop-
star/tree decompositions that presents several advantages
when compared to these schemes. Quasi-Helmholtz projec-
tors have been used to cure the low-frequency breakdowns
of several formulations, however for the sake of readabil-
ity and conciseness, we will only detail their application
to the standard EFIE where it is more straightforward, but
will point to relevant papers describing their applications to
other well-known formulations. Preconditioning the original
system (15) with matrices of the form

$$
P := \alpha P_{\Lambda H} + \beta P_{\Sigma}, \quad (70)
$$

where, following a frequency analysis similar to the one
used for loop-star/tree decompositions, an optimal coefficient
choice can be found to be α ∝ k⁻¹⁄₂ and β ∝ k¹⁄₂, that is,

$$
P_k := \sqrt{C/k} P_{\Lambda H} + i \sqrt{k/C} P_{\Sigma}, \quad (71)
$$

resulting in a new system of equations

$$
P_k T_k P_k y = P_k e^i, \quad (72)
$$

where $P_k y = j$. The constant $C$ can be obtained by maximizing the components of the solution current that are recovered [60] and by enforcing an equal contribution of the vector and scalar potential components; this results in

$$
C := \sqrt{\frac{\| T_{\Phi,k} \|}{\| P_{\Lambda H} T_{A,k} P_{\Lambda H} \|}} . \quad (73)
$$

The analysis of the conditioning effect of the projector can
mimic the strategy used for the loop-star decomposition. In
particular, the EFIE preconditioned with the projectors has
a frequency-independent limit

$$
\lim_{k \to 0} P_k T_k P_k \to i CP_{\Lambda H} T_{A,0} P_{\Lambda H} + i /CT_{\Phi,0}, \quad (74)
$$

where we used that $P_{\Sigma} T_{\Phi,k} P_{\Sigma} = T_{\Phi,k}$ and $P_{\Lambda} T_{\Phi,k} = T_{\Phi,k} P_{\Lambda} = 0$ and thus

$$
\lim_{k \to 0} \mathrm{cond}(P_k T_k P_k) = \gamma, \qquad (75)
$$

where γ is a frequency independent constant. This approach
can be proved to simultaneously solve the problem of catas-
trophic round-off errors in both the current and the right-hand
side of the EFIE [60]. Finally, the use of the projectors has
clear advantages in terms of conditioning with respect to
the use of loop-star or related decompositions that can be
seen in Figure 6. The impact on current and right-hand side
cancellation effects can be observed in Figures 7 and 8.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 7. Comparison of the solenoidal part of the surface current density induced on a sphere of radius 1 m discretized with an average edge length of 0.3 m at $f = 10^{-20}$ Hz, computed with different formulations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8. Comparison of the non-solenoidal part of the surface current density induced on a sphere of radius 1 m discretized with an average edge length of 0.3 m at $f = 10^{-20}$ Hz, computed with different formulations.</figcaption></figure>

### 3) OTHER STRATEGIES FOR THE EFIE LOW-FREQUENCY REGULARIZATION

From previous sections, it is clear that the main drawbacks of loop-star/tree decompositions reside in their constant-in-frequency, but still high, condition number and also in the need to be enriched with global loop functions [67], [68]. Both of these drawbacks can be overcome by the use of quasi-Helmholtz projectors, as explained above, but other schemes can alternatively be used as effective cures for one or both of the drawbacks above. By using a rearranged non-solenoidal basis, for example, the conditioning of a loop-star or a loop-tree preconditioned EFIE could be further improved [61]. Moreover, to avoid the construction of global loops on multiply-connected geometries, formulations have been presented that consider the saddlepoint formulation of the EFIE [69], [70], where the charge is introduced as unknown, in addition to the current in the RWG basis. The most notable are the current-charge formulation [71]

and the augmented EFIE [72]. However, these formulations are, in general, not free from round-off errors in the current or the right-hand side so that, for example, perturbation methods need to be used [58] for further stabilization. An alternative to the perturbation method is the augmented EFIE with normally constrained magnetic field and static charge extraction, which includes a boundary integral equation for the normal component of magnetic field [73]. A disadvantage of current-charge formulations is the introduction of an additional unknown, the charge; hence, methods have been presented to save memory by leveraging nodal functions [74]. An entirely different approach is used in [75], where a closed-form expression of the inverse of the EFIE system matrix is derived based on eigenvectors and eigenvalues of the generalized eigenvalue problem.

Another class of strategies forfeits the EFIE approaches; instead, they are based on potential formulations [76], [77]. These formulations are low-frequency stable on simply- and multiply-connected without the need for searching global loops. The potential-based approaches [76], [77] are also dense-discretization stable. This property is shared with hierarchical basis and Calderón-type preconditioners. With a few exceptions [78], hierarchical basis preconditioners are based on explicit quasi-Helmholtz decomposition [79]–[84], since it then suffices to find a hierarchical basis for scalar-valued functions. While they yield an overall improved condition number with respect to classical loop-star and loop-tree approaches, they require the search for global loops on multiply-connected geometries; a suitable combination with quasi-Helmholtz projectors has been shown to alleviate the need for this search [85]. Calderón-type preconditioners will be discussed in the next section in greater detail. At this point, we are content to say that standard Calderón preconditioned EFIEs have a spectral behavior similar to that of the MFIE operator: thus certain low-frequency issues that plague the MFIE (and which we discuss in the next subsection) persist in the Calderón preconditioned EFIE. Initially remedies relied on combining Calderón preconditioners with loop-star preconditioners. However, the Gram matrix becomes ill-conditioned and global loops must be explicitly recovered [86], [87]. As we will show in Section V, this can be avoided by using quasi-Helmholtz projectors [64].

### 4) HANDLING OF THE RIGHT-HAND SIDE AND FIELD COMPUTATION

As we have already mentioned, a well-conditioned discretization alone is not sufficient to accurately compute $j$: the right-hand side suffers typically from numerical inaccuracies due to finite integration precision and from round-off errors. The main reason for this is that the quasi-Helmholtz components scale differently in frequency. As an example, for the case of the plane-wave excitation, the asymptotic behavior is noted in (51)–(54).

Strategies have been presented in the past to yield stable discretizations of the right-hand side [20], [61], which work with arbitrary right-hand side excitations. For the plane-wave
excitation, a simple solution is to not only compute $e$ as in (18), but also an $e_{\text{extracted}}$, where the static contribution is extracted. We obtain this by replacing $e^{i\hat{k}\cdot r}$ with $e^{i\hat{k}\cdot r} - 1$, where $\hat{k}$ denotes the direction of propagation. Then

$$ \sqrt{C/k} P_{\Lambda H} e_{\text{extracted}} + i \sqrt{Ck} P_{\Sigma} e \quad (76) $$

is a stable discretization of the preconditioned right-hand side. To obtain a stable discretization for small arguments of the exponential function, the subtraction in $e^{i\hat{k}\cdot r} - 1$ should be replaced by a Taylor series, where the static part is omitted.

Similarly, the far-field cannot be computed by simply evaluating

$$ \int_{\Gamma} \sum_{n=1}^{N} [\hat{j}]_n f_n(\mathbf{r}') e^{-i \hat{k} \cdot \mathbf{r}'} dS(\mathbf{r}'), \quad (77) $$

where $j = P_k y$ from (72), $\hat{r} = r/|r|$. On the one hand, by computing the unknown vector of the unpreconditioned formulation $j = P_k y$, the different asymptotic behavior in $k$ of the quasi-Helmholtz components of $j$ as denoted in (57)–(60) would lead to a loss of the solenoidal/quasi-harmonic components in the static limit due to finite machine precision. Thus for the field computation, one should keep the unpreconditioned components of $j$ separately, that is, $j_{\text{sol-qhar}} = \sqrt{C/k} P_{\Lambda H} y$ and $j_{\text{nsol}} = i \sqrt{Ck} P_{\Sigma} y$. On the other hand, it has been pointed out that also the far-field computation suffers from round-off errors [88]. To avoid these, we compute the far-field in two steps: we compute the contribution of $j_{\text{nsol}}$ to the far-field by evaluating

$$ E_{\text{nsol}}^{\text{far}}(\mathbf{r}) = \sum_{n=1}^{N} [j_{\text{nsol}}]_n \int_{\Gamma} f_n(\mathbf{r}') e^{-i \hat{k} \cdot \mathbf{r}'} dS(\mathbf{r}') \quad (78) $$

and the contribution of $j_{\text{sol-qhar}}$ by

$$ E_{\text{sol-qhar}}^{\text{far}}(\mathbf{r}) = \sum_{n=1}^{N} [j_{\text{sol-qhar}}]_n \int_{\Gamma} f_n(\mathbf{r}') (e^{-i \hat{k} \cdot \mathbf{r}'} - 1) dS(\mathbf{r}'), \quad (79) $$

where a Taylor-series expansion should be used for small arguments of the exponential; then

$$ E^{\text{far}}(\mathbf{r}) = E_{\text{nsol}}^{\text{far}}(\mathbf{r}) + E_{\text{sol-qhar}}^{\text{far}}(\mathbf{r}). \quad (80) $$

Also for the near-field computation, the separation in $j_{\text{sol-qhar}}$ and in $j_{\text{nsol}}$ must be maintained, the static contribution removed from the Green's function, and, in addition, the divergence of the scalar potential explicitly enforced by omitting it.

## B. MAGNETIC FIELD INTEGRAL EQUATION

The MFIE has, other than in the Green's function kernel, no explicit dependency on $k$ and should thus be expected to remain well-conditioned in frequency for $k \to 0$. Indeed, for simply-connected geometries $\Gamma$, we have $\text{cond}(M_k) = O(1)$ when $k \to 0$. In the case of multiply-connected geometries,

the MFIE operator exhibits a nullspace associated with the toroidal (for the exterior MFIE) or poloidal loops (for the interior MFIE) in the static limit [89]–[91]. This leads to an ill-conditioned system matrix [91]; in the following, we are going to show that $\text{cond}(M_k) \ge C/k^2$ for some constant $C \in \mathbb{R}^+$.

To prove this result on the condition number, we will consider the low-frequency behavior of block matrices that result from a discretization of the MFIE with a loop-star basis. For the analysis, we must however distinguish two types of harmonic functions, the poloidal and the toroidal loops. If $\Omega_-$ has genus $g$, then the space $H_H(\Omega_-)$ defined by the harmonic functions in $\Omega_-$ and the space $H_H(\Omega_+)$ defined by the harmonic functions in $\Omega_+$ have both dimension $g$. The space defined by $H_{\hat{H}_P}(\Gamma) := \hat{n} \times H_H(\Omega_-) |_{\Gamma}$ are the poloidal loops and the space defined by $H_{\hat{H}_T}(\Gamma) := \hat{n} \times H_H(\Omega_+) |_{\Gamma}$ are the toroidal loops [92], [93]. $H_{\hat{H}_P}(\Gamma)$ has been shown [89] to be the nullspace of $M_0^-$ and $H_{\hat{H}_T}(\Gamma)$ the nullspace of $M_0^+$ and that $\hat{n} \times \hat{H}_{Tn} \in H_{\hat{H}_P}(\Gamma)$ and $\hat{n} \times \hat{H}_{Pn} \in H_{\hat{H}_T}(\Gamma)$.

We need to address how quasi-harmonic functions formed from primal (RWG) or dual (CW/BC) functions are related to harmonic functions. On the one hand, neither with RWG nor with BC functions, we can find linear combinations that are in $H_{\hat{H}_P}(\Gamma)$ or in $H_{\hat{H}_T}(\Gamma)$, a consequence of the fact that these functions are not curl-conforming, as mentioned in Section IV. On the other hand, quasi-harmonic functions $H_n$ (and dual quasi-harmonic functions $\tilde{H}_n$, respectively) are associated with the holes and handles of the geometry. They are not, unlike the locally defined loop functions $\Lambda_n$, derived from a continuous scalar potential on $\Gamma$. Together with the fact that $H_n$ and $\tilde{H}_n$ are solenoidal but not irrotational (since the RWG/BC functions are not curl-conforming), (32) implies that quasi-harmonic loops are linear combinations of solenoidal and harmonic functions (i.e., quasi-harmonic functions are harmonic functions with solenoidal perturbation). Clearly, any quasi-harmonic basis $\{H_n\}_{n=1}^{2g}$ can be rearranged into two bases, where one basis $\{H_{Tn}\}_{n=1}^{g}$ is orthogonal to poloidal loops and the other basis $\{H_{Pn}\}_{n=1}^{g}$ is orthogonal to toroidal loops. In the following, $\hat{H}_{Tn} \in H_{\hat{H}_T}(\Gamma)$ denote the harmonic toroidal and $\hat{H}_{Pn} \in H_{\hat{H}_P}(\Gamma)$ denote the harmonic poloidal basis functions on $\Gamma$, while $H_{Tn}, H_{Pn} \in X_f$ and $\tilde{H}_{Tn}, \tilde{H}_{Pn} \in X_f$ are their quasi-harmonic counterparts.

In [91], scalings were reported of the blocks of the system matrix of $M_k^+$ in terms of a quasi-Helmholtz decomposition

$$ \hat{n} \times \tilde{\Sigma} \begin{bmatrix} \Lambda & \Sigma & H_T & H_P \\ \hat{n} \times \tilde{\Lambda} & O(k^2) & O(k^2) & O(k^2) \\ \hat{n} \times \tilde{H}_P & O(1) & O(1) & O(1) \\ \hat{n} \times \tilde{H}_T & O(k^2) & O(1) & O(k^2) \end{bmatrix}. \quad (81) $$

An analogous result can be obtained for $M_k^-$ with exchanged roles for poloidal and toroidal loops. To observe such a frequency behavior, it is necessary that the testing functions are curl-conforming. Indeed, for the historical MFIE tested with RWG functions, the scalings are not observed [94].
To derive the asymptotic behavior for $k \to 0$ of the block matrices, we start by considering the Taylor series of the Green's function kernel of the MFIE

$$ 
\begin{aligned}
\nabla G(\mathbf{r}, \mathbf{r}') &= \frac{\mathbf{R}}{4\pi R^3} (i k R - 1) e^{i k R} && (82) \\
&= \frac{\mathbf{R}}{4\pi R^3} \left( -1 + \frac{1}{2}(i k)^2 R^2 + \frac{1}{3}(i k)^3 R^3 + \dots \right), && (83)
\end{aligned}
$$

where $R = |\mathbf{r} - \mathbf{r}'|$ and $\mathbf{R} = \mathbf{r} - \mathbf{r}'$. An $O(k^2)$-scaling is observed for a block matrix if the contribution due to the static term in the Taylor series vanishes. For $\mathbf{\Sigma}_n$ as expansion function or $\hat{\mathbf{n}} \times \tilde{\mathbf{\Lambda}}$ as testing function, the static contribution does not vanish and thus we conclude that the scalings of the blocks in the second column and the second row of (81) are constant in $k$.

We now consider the static MFIE, which models the magnetostatic problem, where $j_Γ$ is either solenoidal or harmonic. In this case, we have

$$ \langle \hat{\mathbf{n}} \times \tilde{\mathbf{\Sigma}}_m, \mathcal{M}_0^{\pm} \mathbf{\Lambda}_n \rangle_{\Gamma} = 0, \quad (84) $$

$$ \langle \hat{\mathbf{n}} \times \tilde{\mathbf{\Sigma}}_m, \mathcal{M}_0^{\pm} \mathbf{H}_{Tn} \rangle_{\Gamma} = 0, \quad (85) $$

$$ \langle \hat{\mathbf{n}} \times \tilde{\mathbf{\Sigma}}_m, \mathcal{M}_0^{\pm} \mathbf{H}_{Pn} \rangle_{\Gamma} = 0. \quad (86) $$

This can be seen by considering that a solenoidal function $\mathbf{\Lambda}$ has a corresponding scalar potential

$$ \mathbf{\Lambda} = \mathrm{curl}_{\Gamma} \Phi. \quad (87) $$

Furthermore, we have the equality

$$ \langle \hat{\mathbf{n}} \times \mathbf{\Lambda}, \mathcal{M}_0^{\pm} j_\Gamma \rangle_\Gamma = - \langle \mathbf{\Lambda}, \hat{\mathbf{n}} \times \mathcal{M}_0^{\pm} j_\Gamma \rangle_\Gamma. \quad (88) $$

Inserting (87) for the testing function in the right-hand side of (88) and using the fact that $\mathrm{curl}_\Gamma$ is the adjoint operator of $\mathrm{curl}_\Gamma$, we obtain

$$ 
\begin{aligned}
\langle \hat{\mathbf{n}} \times \mathbf{\Lambda}, \mathcal{M}_0^{\pm} j_\Gamma \rangle_\Gamma &= \langle \Phi, \mathrm{curl}_\Gamma \hat{\mathbf{n}} \times (\mathcal{M}_0^{\pm} j_\Gamma) \rangle_\Gamma \\
&= -\langle \Phi, \mathrm{curl}_\Gamma \hat{\mathbf{n}}_{\perp}^{\pm} \rangle_\Gamma,
\end{aligned}
\quad (89) $$

where $\hat{\mathbf{n}}_{\perp}^{\pm} := \hat{\mathbf{n}} \times (\mathcal{M}_0^{\pm} j_\Gamma)$ is the (rotated) tangential component of the magnetic field

$$ \mathbf{h}(\mathbf{r}) = \mathrm{curl} \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}') j_{\Gamma}(\mathbf{r}') \, dS(\mathbf{r}') \quad \text{for } \mathbf{r} \in \Omega_{\pm}, \quad (90) $$

that is, $\hat{\mathbf{n}}_{\perp}^-$ is the (rotated) tangential component of $\mathbf{h}$ when $\Gamma$ is approached from within $\Omega_-$, and $\hat{\mathbf{n}}_{\perp}^+$ when $\Gamma$ is approach from within $\Omega_+$. We recall that (90) is obtained by finding a vector potential $\mathbf{a}$ such that $\mathbf{h} = \mathrm{curl}\mathbf{a}$ and noting that $\mathrm{curl}\mathrm{curl}\mathbf{a} = j_\Gamma$ under the assumption that $\mathrm{div}j_\Gamma = 0$ [40, Ch. 6.1].

From

$$ \mathrm{curl}_\Gamma (\mathbf{v}_\Gamma^{\pm}(\mathbf{r})) = \lim_{\Omega_{\pm} \ni \mathbf{r}' \to \mathbf{r}} (\mathrm{curl}(\mathbf{v}(\mathbf{r}'))) \cdot \hat{\mathbf{n}}(\mathbf{r}') \quad (91) $$

and from [40, eq. (6.17)]

$$ \mathrm{curl} \mathbf{v} = \{\mathrm{curl} \mathbf{v}\} + \hat{\mathbf{n}} \times (\mathbf{v}^{+} - \mathbf{v}^{-}) \delta_{\Gamma} \quad (92) $$

where curly braces {} mean that this part is evaluated only in $\Omega^{\pm}$ and $\delta_\Gamma$ is the surface Dirac delta function, together with (90), we have

$$ \mathrm{curl}_\Gamma \hat{\mathbf{n}}_{\perp}^{\pm} = \mathrm{curl} \hat{\mathbf{n}} \cdot \hat{\mathbf{n}} = 0. \quad (93) $$

Next we will establish that $\langle \hat{\mathbf{n}} \times \mathbf{H}_{Pm}, \mathcal{M}_0^{\pm} \mathbf{\Lambda}_n \rangle_\Gamma = 0$ and $\langle \hat{\mathbf{n}} \times \mathbf{H}_{Tm}, \mathcal{M}_0^{\pm} \mathbf{\Lambda}_n \rangle_\Gamma = 0$. First, note that the exterior MFIE operator has the mapping properties [89]

$$ \mathcal{M}_0^{\pm} \hat{\mathbf{H}}_{Tn} = \mathbf{0} \quad (94) $$

and

$$ \mathcal{M}_0^{\pm} \hat{\mathbf{H}}_{Pn} = \hat{\mathbf{H}}_{Pn}, \quad (95) $$

while for the interior MFIE operator, we have the mappings properties

$$ \mathcal{M}_0^{-} \hat{\mathbf{H}}_{Pn} = \mathbf{0} \quad (96) $$

and

$$ \mathcal{M}_0^{-} \hat{\mathbf{H}}_{Tn} = \hat{\mathbf{H}}_{Tn}. \quad (97) $$

Furthermore, for any two surface functions $f, g$, we have [89, Sec. 5]

$$ 
\begin{aligned}
\langle \mathbf{g}, \mathcal{K}f \rangle_{\Gamma} &= \langle -\hat{\mathbf{n}} \times (\mathcal{K}(\mathbf{g} \times \hat{\mathbf{n}})), f \rangle_{\Gamma} \\
&= \langle -\mathcal{K}(\mathbf{g} \times \hat{\mathbf{n}}), f \times \hat{\mathbf{n}} \rangle_{\Gamma}.
\end{aligned}
\quad (98) $$

Then we find

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Tn}, \mathcal{M}_0^{+} \mathbf{\Lambda} \rangle_{\Gamma} = \langle \mathcal{M}_0^{-} \hat{\mathbf{H}}_{Tn}, \mathbf{\Lambda} \times \hat{\mathbf{n}} \rangle_{\Gamma} \stackrel{(97)}{=} 0, \quad (99) $$

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Pn}, \mathcal{M}_0^{+} \mathbf{\Lambda} \rangle_{\Gamma} = \langle \mathcal{M}_0^{-} \hat{\mathbf{H}}_{Pn}, \mathbf{\Lambda} \times \hat{\mathbf{n}} \rangle_{\Gamma} \stackrel{(96)}{=} 0, \quad (100) $$

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Tn}, \mathcal{M}_0^{-} \mathbf{\Lambda} \rangle_{\Gamma} = \langle \mathcal{M}_0^{+} \hat{\mathbf{H}}_{Tn}, \mathbf{\Lambda} \times \hat{\mathbf{n}} \rangle_{\Gamma} \stackrel{(94)}{=} 0, \quad (101) $$

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Pn}, \mathcal{M}_0^{-} \mathbf{\Lambda} \rangle_{\Gamma} = \langle \mathcal{M}_0^{+} \hat{\mathbf{H}}_{Pn}, \mathbf{\Lambda} \times \hat{\mathbf{n}} \rangle_{\Gamma} \stackrel{(95)}{=} 0, \quad (102) $$

where we used (94)-(97) and the orthogonality of harmonic and irrotational functions $\hat{\mathbf{n}} \times \mathbf{\Lambda}$. Now consider that $\hat{\mathbf{H}}_{Pn} = \hat{\mathbf{H}}_{Pn} + \tilde{\mathbf{\Sigma}}_{Pn}$ and $\hat{\mathbf{H}}_{Tn} = \hat{\mathbf{H}}_{Tn} + \tilde{\mathbf{\Sigma}}_{Tn}$, where $\tilde{\mathbf{\Sigma}}_{P/Tn}$ is the respective perturbation. Thus by taking into account (84), we have

$$ \langle \hat{\mathbf{n}} \times \tilde{\mathbf{H}}_{Pm}, \mathcal{M}_0^{\pm} \mathbf{\Lambda}_n \rangle_{\Gamma} = 0, \quad (103) $$

$$ \langle \hat{\mathbf{n}} \times \tilde{\mathbf{H}}_{Tm}, \mathcal{M}_0^{\pm} \mathbf{\Lambda}_n \rangle_{\Gamma} = 0. \quad (104) $$

From now on, we will only consider the harmonic functions $\hat{\mathbf{H}}_{Tn}$ and $\hat{\mathbf{H}}_{Pn}$ instead of their quasi-harmonic counterpart since, as we have seen, the solenoidal perturbation will always vanish. Then for $\mathcal{M}_0^+$, we have

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Tm}, \mathcal{M}_0^{+} \hat{\mathbf{H}}_{Tn} \rangle_{\Gamma} = \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Pm}, \mathcal{M}_0^{+} \hat{\mathbf{H}}_{Tn} \rangle_{\Gamma} \stackrel{(94)}{=} 0, \quad (105) $$

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Tm}, \mathcal{M}_0^{-} \hat{\mathbf{H}}_{Pn} \rangle_{\Gamma} = \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Pm}, \mathcal{M}_0^{-} \hat{\mathbf{H}}_{Pn} \rangle_{\Gamma} \stackrel{(96)}{=} 0 \quad (106) $$

due to the nullspace. Finally, we have

$$ \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Pm}, \mathcal{M}_0^{+} \hat{\mathbf{H}}_{Pn} \rangle_{\Gamma} \stackrel{(95)}{=} \langle \hat{\mathbf{n}} \times \hat{\mathbf{H}}_{Pm}, \hat{\mathbf{H}}_{Pn} \rangle_{\Gamma} = 0, \quad (107) $$
as $\hat{n} \times \hat{H}_{Pm} \in H_{\hat{H}_T}(\Gamma)$ and this space is orthogonal to $H_{\hat{H}_P}(\Gamma)$. Likewise, we can conclude that there is at least one $m$ such that

$$ \langle \hat{n} \times \hat{H}_{Tm}, M_0^+ \hat{H}_{Pn} \rangle_{\Gamma} \neq 0. \quad (108) $$

Analogously, we obtain

$$ \langle \hat{n} \times \hat{H}_{Pm}, M_0^+ \hat{H}_{Tn} \rangle_{\Gamma} \neq 0. \quad (109) $$

This has established the scalings of the different block matrices composing $M_k^\pm$.

To prove the ill-conditioning in $k$, we consider the static limit $k \to 0$, where $j_\Sigma = 0$ and

$$ \|j_\Lambda\| = \|j_{\hat{H}_T}\| = \|j_{\hat{H}_P}\| = O(1), $$

following the argumentation in [95]. Forming the matrix-vector product, we find

$$ M_k^+ \begin{bmatrix} j_\Lambda \\ \mathbf{0} \\ j_{\hat{H}_T} \\ j_{\hat{H}_P} \end{bmatrix} = \begin{bmatrix} O(k^2) \\ O(1) \\ O(k^2) \\ O(1) \end{bmatrix}. \quad (110) $$

Thus for $k \to 0$, we have the reduced system

$$ \underbrace{\begin{bmatrix} M_{0,\hat{n}\times\tilde{\Lambda},\Lambda}^{+} & M_{0,\hat{n}\times\tilde{\Lambda},H_T}^{+} & M_{0,\hat{n}\times\tilde{\Lambda},H_P}^{+} \\ \mathbf{0} & \mathbf{0} & M_{0,\hat{n}\times\tilde{H}_T,H_P}^{+} \end{bmatrix}}_{(N_{\tilde{\Lambda}}+N_{\tilde{H}_T}) \times (N_{\Lambda}+N_{H_T}+N_{H_P})} \times \underbrace{\begin{bmatrix} j_\Lambda \\ j_{\hat{H}_T} \\ j_{\hat{H}_P} \end{bmatrix}}_{(N_{\tilde{\Lambda}}+N_{\hat{H}_P}) \times 1} = \underbrace{\begin{bmatrix} O(1) \\ O(1) \end{bmatrix}}_{(N_{\tilde{\Lambda}}+N_{\hat{H}_P}) \times 1}, \quad (111) $$

where $N_X$ is the dimensionality of the respective finite element space of the basis functions $X_n$. Since $N_{\tilde{\Lambda}} = N_\Lambda$, this implies that $M_0^+$ has an $N_{\hat{H}_P} = g$-dimensional nullspace.

To obtain a lower bound on the condition number dependency of $M_k^\pm$ on $k$, we consider that for all $j_0^\pm \in \text{null}(M_0^\pm)$, we have $\|M_k^\pm j_0^\pm\| = O(k^2)$ for $k \to 0$. Thus $\text{cond}(M_k^\pm) \gtrsim k^2$.

## 1) QUASI-HELMHOLTZ PROJECTORS

It is not possible to cure the low-frequency breakdown of the MFIE on multiply-connected geometries with projectors. The reason is that only a part of the harmonic Helmholtz subspace is affected. What, however, is more critical for the MFIE in low-frequency scenarios than the ill-conditioning, is that round-off errors and finite precision of integration rules lead to non-vanishing static components. Equation (81) shows that several block matrix entries should vanish when $k \to 0$, for example, we should observe $\|P_\Sigma M_0^\pm P_\Lambda\| = 0$; this is not observed in practice due to errors, for example, in numerical integration. To obtain accurate solutions even for $k \to 0$, the contribution of the static kernel $K_0$ must be removed. On simply-connected geometries, it is sufficient to explicitly set the term $P_\Sigma M_0^\pm P_\Lambda$ to zero. On multiply-connected geometries, this is not possible due to the different scaling of the block matrices in (81) associated with poloidal

and toroidal loops (i.e., it could be cured with projectors on the poloidal and toroidal quasi-Helmholtz subspace. These are, however, numerically expensive to construct).

Quasi-Helmholtz projectors can be used to ensure the vanishing of the static components that ought to be, according to (81), zero, but are not due to numerical integration errors. This is obtained by a combination of the interior and the exterior MFIE operator together with quasi-Helmholtz projectors as outlined in [96], that is,

$$ \tilde{P}_k \left( M_k^- G_{\hat{n}\times\tilde{f},f}^{-1} M_k^+ - \underbrace{P_{\Sigma H} M_0^- G_{\hat{n}\times\tilde{f},f}^{-1} M_0^+ P_{\Lambda H}}_{\text{remove static kernel contribution}} \right) P_k Y = \tilde{P}_k M_k^- G_{\hat{n}\times\tilde{f},f}^{-1} h^i, \quad (112) $$

where

$$ \tilde{P}_k = 1/\sqrt{k} P_{\Sigma H} + i \sqrt{k} P_{\Lambda}, \quad (113) $$

is the counterpart of $P_k$ on the dual mesh. The formulation in (112) symmetries the behavior on the quasi-harmonic Helmholtz subspace and ensures that the MFIE can be solved accurately for $k \to 0$ since no catastrophic round-off errors occur in the right-hand side and $j$, and since the numerical integration error that leads to a nonphysical contribution of the static kernel is removed. Moreover, as we will see in the next section, it can be directly combined with an EFIE resulting in a stable CFIE formulation.

It should be noted that for a stable implementation of (112), one should, instead of removing the static kernel contribution by a subtraction, build the system matrix without this contribution in the first place. To this end, we define the operator

$$ (\mathcal{K}'_k j_\Gamma)(r) := -\hat{n} \times \int_\Gamma \nabla \frac{e^{i k |r - r'|} - 1}{4\pi |r - r'|} \times j_\Gamma \, dS(r') \quad (114) $$

and its discretization as

$$ [\mathcal{K}'_k]_{nm} := \langle \hat{n} \times \tilde{f}_n, \mathcal{K}'_k f_m \rangle_{\Gamma}, \quad (115) $$

where we note that for small $k|\mathbf{r} - \mathbf{r}'|$ a Taylor-series should be used to remove the static contribution. Then the overall system matrix is

$$ \begin{aligned} & 1/\sqrt{k} P_{\Sigma H} \left( (-\mathcal{K}'_k) G_{\hat{n}\times\tilde{f},f}^{-1} M_k^+ + M_k^- G_{\hat{n}\times\tilde{f},f}^{-1} \mathcal{K}'_k \right) \left( 1/\sqrt{k} P_{\Lambda H} \right) \\ & + 1/\sqrt{k} P_{\Sigma H} M_k^- G_{\hat{n}\times\tilde{f},f}^{-1} M_k^+ \left( i \sqrt{k} P_\Sigma \right) \\ & + i \sqrt{k} P_\Lambda M_k^- G_{\hat{n}\times\tilde{f},f}^{-1} M_k^+ \left( 1/\sqrt{k} P_{\Lambda H} \right) \\ & + i \sqrt{k} P_\Lambda M_k^- G_{\hat{n}\times\tilde{f},f}^{-1} M_k^+ \left( i \sqrt{k} P_\Sigma \right), \end{aligned} \quad (116) $$

which evidently does not contain the static kernel contribution $P_{\Sigma H} M_0^- G_{\hat{n}\times\tilde{f},f}^{-1} M_0^+ P_{\Lambda H}$.
## 2) HANDLING OF THE RIGHT-HAND SIDE

Similarly to the case of EFIE, the right-hand side requires special treatment. First, we note

$$ \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} h^i = \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Sigma H} h_{\text{extracted}}^i + \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i, \quad (117) $$

where $h_{\text{extracted}}^i$ is computed analogously to $e_{\text{extracted}}^i$. Again, the equality is only true in exact arithmetic; the first term on the right-hand side needs no further treatment as the asymptotically constant blocks of $\mathcal{K}_k^-$ will dominate. The second term, however, needs to be further stabilized. We note

$$ \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i = \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} (P_{\Lambda} + P_{\Sigma H}) \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \\ = \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} P_{\Lambda} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i, \quad (118) $$

where we used $P_{\Sigma H} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} = \mathbf{0}$ from Lemma 1. Using (81), we have

$$ \begin{aligned} \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} P_{\Lambda} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i & = 1/\sqrt{k} P_{\Sigma H} (-K'_k) P_{\Lambda} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \\ & \quad + i \sqrt{k} P_{\Lambda} \mathcal{M}_k^{-} P_{\Lambda} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \\ & = 1/\sqrt{k} P_{\Sigma H} (-K'_k) \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \\ & \quad + i \sqrt{k} P_{\Lambda} \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i. \end{aligned} \quad (119) $$

Summarizing, we obtain for the right-hand side

$$ \begin{aligned} \tilde{\mathcal{P}}_k \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Sigma H} h_{\text{extracted}}^i & \\ & + 1/\sqrt{k} P_{\Sigma H} (-K'_k) \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \\ & + i \sqrt{k} P_{\Lambda} \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i, \end{aligned} \quad (120) $$

which can be rearranged into

$$ \begin{aligned} & 1/\sqrt{k} P_{\Sigma H} \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Sigma H} h_{\text{extracted}}^i \\ & + 1/\sqrt{k} P_{\Sigma H} (-K'_k) \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \\ & + i \sqrt{k} P_{\Lambda} \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} h^i \end{aligned} \quad (121) $$

using

$$ \begin{aligned} & \left\| i \sqrt{k} P_{\Lambda} \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Sigma H} h_{\text{extracted}}^i \right\| \\ & \ll \left\| i \sqrt{k} P_{\Lambda} \mathcal{M}_k^{-} \mathcal{G}_{\hat{n} \times \tilde{\mathcal{f}}_f}^{-1} P_{\Lambda} h^i \right\|. \end{aligned} \quad (122) $$

## 3) OTHER LOW-FREQUENCY REGULARIZATIONS OF THE MFIE

When it comes to the low-frequency regularization of the MFIE, the body of literature falls shorter than that of the EFIE. In light of the discussion of the previous section, it is evident that a correct low-frequency behavior can only be expected when the MFIE operator is tested with dual basis functions [43], even though improvements can be expected

for the standard MFIE by using, for example, a perturbation method [94]. Instead of quasi-Helmholtz projectors, explicit quasi-Helmholtz decompositions can be used to ensure the vanishing of the static components of the discretized MFIE. Standard bases such as the loop-star worsen, however, the conditioning in $h$ [62] since the identity operator is turned into a discretized Laplace-Beltrami operator. This can be avoided with hierarchical bases as denoted in [97] in the context of the CFIE, or by using quasi-Helmholtz projectors as outlined in Section IV-C1.

A remaining issue is the ill-conditioning due to the quasi-static nullspace of either the toroidal or the poloidal loops of the exterior or the interior MFIE operator, respectively. Strategies that have been presented in the past typically introduce an extra condition based on vector potential considerations and on an explicit detection of global loops to remove the ill-conditioning [98], [99].

To obtain a wideband stable solver, an extension of the preconditioning schemes to the CFIE are necessary. For Calderón preconditioning, we will discuss this in greater detail in the next section, but a few words should be said on the other techniques. As mentioned earlier, explicit quasi-Helmholtz decompositions can be readily extended to the CFIE, however, some such as the loop-star basis will worsen the conditioning of the CFIE. Hierarchical basis preconditioners can be extended without the worsening [97]. An augmented CFIE has been presented [100], however, it does not appear to resolve numerical round-off losses.

In concluding this section on the EFIE and the MFIE low-frequency behavior, a word should be said regarding fast methods. Typically, a fast method is employed to accelerate the matrix-vector products. One of the most popular choices, the MLFMM, is not low-frequency stable [6]. Even if a low-frequency stable method is chosen, such as the ACA, catastrophic round-off errors can appear in the EFIE due to the different scaling of the vector and scalar potential part in frequency. This can be resolved by simply storing those contributions separately [101]. More memory efficient approaches have been presented [102] at the price of a more complicated implementation.

## V. DENSE-DISCRETIZATION SCENARIOS

The dense-discretization breakdown is sometimes confused with the low-frequency breakdown. They are, however, two distinct phenomena: the dense-discretization breakdown implies that the condition number grows when $h \to 0$.

The condition number of the EFIE is known to grow as

$$ \mathrm{cond}(\mathcal{T}) = O(1/h^2). \quad (123) $$

Unlike the low-frequency breakdown, it is more difficult to prove this statement rigorously. One possible line of argumentation is based on pseudo-differential operator theory [103].

Another line of argumentation follows from a functional analytic point of view, where the link between the coefficient
<figure><img src="image_11.png" /><figcaption>FIGURE 9. The figure shows $\Sigma_{\max}(T_k)$ and $\Sigma_{\min}(T_k)$ as a function of the maximum spectral index $1/h$ for a sphere of unit radius.</figcaption></figure>

space $C^N$ and the domain and range of the respective operator is considered [104]. Based on inverse inequalities, the norm of a Sobolev space $H^m$ is bounded by the norm of $H^l$, where $l < m$. In the case of the EFIE operator, this approach needs a non-trivial extension as the $H_{\text{div}}^{-1/2}$-space necessitates a Helmholtz decomposition.

To avoid the introduction of either Sobolev space or pseudo-differential operator theory [105], one can also study the singular value spectrum. This is typically only possible for a few canonical geometries and thus the results are not a rigorous proof. They do provide, however, a more intuitive understanding of the underlying physics. In the Appendix, we have included a discussion of the singular value spectrum of the EFIE on a sphere and its link to the breakdowns.

From Appendix A-B, it is evident that $\mathcal{T}_{A,k}$ has singular values that accumulate at zero with a rate of $h$, while $\mathcal{T}_{\Phi,k}$ has singular values that grow to infinity with a rate of $1/h$ resulting in an overall condition number scaling proportional to $1/h^2$. The discretized operator inherits the spectral properties of the analytic EFIE operator, though the effect of the basis function should be removed by considering the spectrum of $G_{f,f}^{-1}T$. Figure 9 shows the maximum and minimum singular values of the EFIE operator discretized on a sphere, confirming that indeed the condition number scales as $1/h^2$.

For the MFIE, we observe that the singular values are bounded from above and below. Indeed, second-kind integral equations such as the MFIE are known to result in well-conditioned system matrices if the expansion and testing functions are $L^2$-stable (i.e., that the Gram matrix is well-conditioned) [106].

For the CFIE, the conditioning improves to

$$ \mathrm{cond}(C_k) = O(1/h), \qquad (124) $$

due to the MFIE part of the CFIE: the identity of the MFIE introduces a lower bound on the singular values, which

annihilates the ill-conditioning due to $\mathcal{T}_{A,k}$. Hence, only $\mathcal{T}_{\Phi}$ contributes to the growing condition number.

## A. CALDERÓN PRECONDITIONING

The Calderón identity

$$ \mathcal{T}_k \circ \mathcal{T}_k = \mathcal{M}_k^{-} \circ \mathcal{M}_k^{+} = -\frac{\mathcal{I}}{4} + \mathcal{K}_k^2 \qquad (125) $$

suggests that the EFIE operator can be turned into a second-kind integral operator, which, as other second-kind integral equations like the MFIE, is well-conditioned in $h$. The stable discretization of (125) is not trivial, however, (see, for example, the pioneering works [103], [107]–[109]) and research for a stable discretization lead to discovery of an alternative to CW functions in the form of the BC functions [45] that leveraged the seminal work in operator preconditioning of Steinbach and Wendland [110].

### 1) CALDERÓN PRECONDITIONING WITH BC FUNCTIONS

The Calderón multiplicative preconditioner (CMP) is an effective and well-conditioned scheme for intricate geometries [46]. With the notation of this paper, the CMP-EFIE can be expressed as

$$ \tilde{\mathcal{T}}_k G_{\hat{n} \times f \hat{f}}^{-1} \mathcal{T}_k j = - \tilde{\mathcal{T}}_k G_{\hat{n} \times f \hat{f}}^{-1} e^i, \qquad (126) $$

where

$$ [\tilde{\mathcal{T}}_k]_{nm} := \langle \hat{n} \times \tilde{f}_n, \mathcal{T}_k \tilde{f}_m \rangle_I. \qquad (127) $$

The preconditioning strategy has been used extensively and successfully in several application scenarios with effective improvements [111]–[113].

However, when the frequency is extremely low (when $k^2$ normalized by the norm ratio of the scalar and vector potentials reaches $10^{-16}$ in double precision), the standard CMP will break down: in fact, from Section IV, we can expect that (126) inherits the low-frequency issues that the MFIE is suffering from on multiply-connected geometries due to the Calderón identity $-\mathcal{I}/4+\mathcal{K}_k^2 = (\mathcal{I}/2+\mathcal{K}_k)(-\mathcal{I}/2+\mathcal{K}_k)$. To overcome the low-frequency breakdown due to the harmonic Helmholtz subspace and due to catastrophic round-off errors, a quasi-Helmholtz decomposition must be introduced.

In [64], quasi-Helmholtz projectors are leveraged yielding

$$ \tilde{\mathcal{P}}_k \tilde{\mathcal{T}}_k \tilde{\mathcal{P}}_k G_{\hat{n} \times f \hat{f}}^{-1} \mathcal{P}_k \mathcal{T}_k \mathcal{P}_k y = - \tilde{\mathcal{P}}_k \tilde{\mathcal{T}}_k \tilde{\mathcal{P}}_k G_{\hat{n} \times f \hat{f}}^{-1} \mathcal{P}_k e^i. \qquad (128) $$

The projectors ensure that the overall system matrix is well-conditioned also on multiply-connected geometries. Figure 10 shows that the condition number of the preconditioned system matrix in (128) is asymptotically bounded. A similar result holds even on the more challenging NASA almond benchmark [114] (Figure 11).

### 2) REFINEMENT-FREE CALDERÓN PRECONDITIONING

The use of BC functions increases, unfortunately, the numerical costs due to their larger support compared with RWG functions (though specific implementations can help offset
<figure><img src="image_12.png" /><figcaption>FIGURE 10. Condition number of the system matrices as a function of 1/h, which is proportional to the maximum spectral index, for a cube with side length 1 m and a frequency of 10⁷ Hz. The labels “CMP-EFIE”, “CFIE”, “RF-CMP-EFIE”, “P-CMP-MFIE”, and “P-CMP-CFIE” refer to the standard Calderón EFIE (126), the conforming CFIE, the refinement-free Calderón preconditioned EFIE (129), and the Calderón preconditioned MFIE and CFIE with projectors (112) and (142).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 11. Conditioning of the different formulations on the NASA almond — scaled to be of 1.1 m in diameter — and a simulating frequency of 5 × 10⁷ Hz.</figcaption></figure>

the overhead [115], [116]). Moreover, no extension to geometries with junctions is available. Fortunately, a Calderón preconditioned EFIE has been obtained without the use of BC functions by indirectly using the Laplace-Beltrami operator as preconditioner. The formulation reads [117]

$$P_o^\dagger T_k^\dagger P_m T P_o i = -P_o^\dagger T_k^\dagger P_m e^i, \quad (129)$$

where the outer matrix $P_o$ is

$$P_o := P_{\Lambda H}/\alpha + i P_g \Sigma / \beta \quad (130)$$

with

$$P_g \Sigma := \Sigma (\Sigma^T \Sigma)^+ G_{\tilde{\Lambda},p}^{-1} \Sigma^T, \quad (131)$$

where $G_{\tilde{\Lambda},p}$ is the Gram matrix of piecewise linear dual basis functions as defined in [45] and patch basis functions with

<figure><img src="image_14.png" /><figcaption>FIGURE 12. Condition number of the system matrices of several formulations, for a sphere of radius 1 m and average edge length of 0.15 m as a function of the frequency f. The label “Yu-CMP-CFIE” refers to the standard Calderón CFIE with modified Yukawa wavenumber (141).</figcaption></figure>

patch height $1/A_i$, where $A_i$ is the area of cell i of the mesh; furthermore, the middle matrix is

$$P_m := P_m \Lambda + P_m \Sigma \quad (132)$$

with the definitions

$$P_m \Lambda := \Lambda G_{\Lambda\Lambda}^{-1} \Lambda^T / \alpha^2 + P_{\Lambda H} / \gamma, \quad (133)$$

$$P_m \Sigma := \Sigma (\Sigma^T \Sigma)^+ G_{pp}^{-1} (\Sigma^T \Sigma)^+ \Sigma / \beta^2, \quad (134)$$

where $G_{\Lambda\Lambda}$ and $G_{pp}$ are the Gram matrices of piecewise linear and piecewise constant basis functions on the primal mesh (for a detailed definition, we refer the reader to [117]), and the scaling coefficients are

$$\alpha = \sqrt[4]{\| P_{\Lambda H} T_{A,k}^\dagger \Lambda G_{\Lambda\Lambda}^{-1} \Lambda^T T_{A,k} P_{\Lambda H} \|_2}, \quad (135)$$

$$\beta = \sqrt[4]{\| P_g^T \Sigma T_{\Phi,k}^\dagger \Sigma (\Sigma^T \Sigma)^+ G_{pp}^{-1} (\Sigma^T \Sigma)^+ \Sigma^T T_{\Phi,k} P_g \Sigma \|_2}, \quad (136)$$

$$\gamma = \|(P_{\Lambda H}/\alpha) T_{A,k}^\dagger P_{\Lambda H} T_{A,k} (P_{\Lambda H}/\alpha)\|_2, \quad (137)$$

where the norms can be estimated with the power iteration method for computing the largest singular value. We note that for $k \to 0$, we have

$$\alpha \propto \sqrt{k}, \quad (138)$$

$$\beta \propto 1/\sqrt{k}, \quad (139)$$

$$\gamma \propto k, \quad (140)$$

which are similar scalings as encountered for the quasi-Helmholtz projectors. A key property of this formulation is that the resulting system matrix is Hermitian, positive definite. This allows the use of CG, which in exact arithmetic guarantees convergence, and has a lower computational cost with respect to CGS or GMRES [53].
<figure><img src="image_15.png" /><figcaption>FIGURE 13. Number of iterations needed by GMRES(50) to converge for tolerance $10^{-4}$ when solving the linear systems arising from different formulations discretized on a cube, for a frequency of $10^7$ Hz. For the RF-CMP from (129), the iterations needed by CG are shown (since for this formulation the system matrix is Hermitian, positive definite).</figcaption></figure>

Figure 13 shows the number of iterations for the refinement-free preconditioned compared with (128), a standard EFIE and a loop-tree preconditioned EFIE. It displays that the number of iterations become bounded independently from $h$, thus confirming the dense-discretization stability of the formulation.

### 3) ALTERNATIVE CALDERÓN-TYPE STRATEGIES

Many alternative Calderón-identity based preconditioning strategies have been presented, most prominently the pioneering works by Adams and Brown [118] and by Christiansen and Nédélec [103], [119]. Often the focus was not only on obtaining a well-conditioned EFIE, but rather a well-conditioned CFIE or combined source integral equation (CSIE) [120], [121] with the pioneering works by Contopaganos et al. [108] and Adams [122]. These schemes do not use dual basis functions, but instead require ad-hoc implementations of modified operators or use Nystrom discretizations.

Another class of preconditioners is based on Dirichlet-to-Neumann operators and on-surface radiation conditions [109], [123]–[125]. As will be discussed in the next section, some of these methods have been extended to handle the high-frequency breakdown. Naturally, these extension are all capable of curing the dense-discretization breakdown. In these schemes, low-frequency issues have not always been addressed in the original papers, but we believe they could be extended to handle those as well.

To stabilize the standard Calderón preconditioner in the static limit, it has been proposed to discretize the EFIE with a loop-star basis first [86], [87], which implies ill-conditioned Gram matrices and the search for global loops. A perturbation method based approach was also suggested [126].

Another strategy employs a discretization with the $H_{\text{div}}$-inner product [127]; combined with a Calderón preconditioner, that cures both the low-frequency and the dense-discretization breakdown; however, the scheme has not yet been extended to multiply-connected geometries.

For open problems, the Calderón preconditioner in (128) yields only a logarithmic bound in $h$. Recently, a formulation has been proposed that yields an $h$-independent bound by leveraging an analytic transformation for the open disc and Lipschitz transformations of the geometry [128].

### B. QUASI-OPTIMAL HIERARCHICAL PRECONDITIONERS

One alternative to Calderón preconditioning for the EFIE are hierarchical basis approaches. Hierarchical bases (also referred to as multiresolution or wavelet methods) have been pioneered in the context of scalar integral equations; they can not only be used as preconditioner, but also to compress the resulting system matrix [129]–[140]. Both the compressibility of the system matrix [78], [141], [142], as well as the preconditioning effect has been demonstrated for the EFIE system matrix [79], [80], [83], [142], [143]. For unstructured meshes, a hierarchically preconditioned EFIE could initially only be obtained for the $T_{\Phi,k}$-part of $T_k$ [81]. Recently, a scheme based on generalized primal and dual Haar prewavelets has been presented, where the condition grows only logarithmically in $N$, and which—in contrast to other hierarchical bases—extends naturally to unstructured meshes [84].

The logarithmic growth renders hierarchical basis preconditioners quasi-optimal compared with Calderón schemes, where the condition number can be asymptotically bounded by a constant. This disadvantage can, in practice, often be compensated by the fact that unlike Calderón schemes, no second multiplication with the EFIE matrix is required. As mentioned in Section IV-A3, hierarchical basis preconditioners can be applied to the CFIE [97] and the search for global loops can be avoided by a suitable combination with quasi-Helmholtz projectors [85].

### C. ALTERNATIVE STRATEGIES

From the alternative strategies that remedy the low-frequency breakdown mentioned in Section IV-A3, only a few are dense-discretization stable without further changes. The most prominent example is the decoupled potential integral equation [76]. The augmented EFIE, on the other hand, can be combined with a Calderón preconditioner [144].

As an alternative to operator preconditioning based strategies, algebraic preconditioners have been proposed. In general, this class of preconditioners will not provide asymptotic bounds on the condition number in $h$. Very often those methods have been designed with electrically large problems in mind and they typically employ the MLFMM: wideband stability is typically not a construction criterion.

The idea of algebraic preconditioners is to obtain an approximation $B$ of the system matrix $A$, where $B^{-1}$ can be
computed explicitly (a direct inverse is obtained) or implicitly (a linear system has to be solved at each iteration step) fast.<fnref n="1" /> A detailed review of algebraic preconditioners for electromagnetic integral equations is provided in [25]. Here, we review only the most well-known techniques.

A typical example of explicit preconditioners are sparse approximations [145]–[148]. A challenge is that if a fast method such as the MLFMM is used, the full system matrix is not available. It is customary to use, for example, the near-field interactions of the MLFMM to construct a preconditioner [148]. While this matrix is indeed sparse, the goal is to obtain a sparse inverse and, in general, the inverse of a sparse matrix is not sparse. This sparse matrix is typically obtained as solution to the minimization problem $||I - A_N B||_F$, where $||\cdot||_F$ is the Frobenius norm and $A_N$ refers to the near-interaction matrices [149]. A filtering of $A_N$ is often necessary to reduce the computational cost [150] and the use of macro basis functions has also been proposed to reduce the problem size [151]. Since the methods are based on the near-interactions, it is local in nature; instead, to consider global interactions, [152] proposes to use inter-polative decompositions, while [149] employs an embedded iterative scheme. Recently, hybrid parallelizations for clusters have been presented [153]. Also, the sparse approximate preconditioners can be combined with other methods such as a two-grid spectral preconditioner [154].

Some of the most prominent among the implicit methods are based on the incomplete LU decomposition [155]–[158]. The costs for obtaining the incomplete LU decomposition are often lower than the costs of constructing a sparse approximate inverse. However, highly indefinite systems pose a significant challenge as they can give rise to ill-conditioned triangular factors [25]. As an alternative, GMRES-based near-zone preconditioning has been proposed [159].

## VI. HIGH-FREQUENCY ISSUES

For resonant frequencies, both the EFIE and the MFIE have nullspaces associated with the interior resonances and are thus ill-conditioned. To overcome the problem of interior resonances, a CFIE or CSIE formulation is typically employed [21], [160].

### A. CALDERÓN-YUKAWA COMBINED FIELD INTEGRAL EQUATION

Calderón preconditioning can also be extended to a combined field framework. A trivial combination of the Calderón preconditioned EFIE (126) and the MFIE from (20) would lead to an overall operator suffering from interior resonances: the identity (125) implies that the Calderón preconditioned EFIE and the MFIE share part of their nullspaces. To obtain a wideband stable CFIE, the EFIE from (126) should be combined with the MFIE from (20) with one important adjustment: the wavenumber $k$ in $\tilde{T}_k$ needs to be made a

complex number. One possibility for a complex $k$ is to use the Yukawa potential resulting in

$$ \begin{aligned} & \left(-\alpha \tilde{T}_{ik} G_{\hat{n} \times \hat{f} \hat{f}}^{-1} T_k + (1-\alpha) M_k^{+} \right) j \\ &= \alpha \tilde{T}_{ik} G_{\hat{n} \times \hat{f} \hat{f}}^{-1} e^i + (1-\alpha) h^i. \end{aligned} \quad (141) $$

Originally, this has been proposed in [108] and later employed in the Calderón multiplicative preconditioner for the CFIE [161]; Section V-A3 includes a discussion of Calderón-type preconditioners for the CFIE and CSIE.

This standard approach is, however, not stable until arbitrarily low frequencies, because of the numerical issues detailed in Section IV. A stable formulation can be obtained by combining the Calderón preconditioned EFIE (128) and the stabilized MFIE from (112). Similar to the preconditioning EFIE matrix, a complex wavenumber should be used for preconditioning MFIE. This results in the overall formulation

$$ \begin{aligned} & \tilde{P}_k (\tilde{T}_{ik} \tilde{P}_k G_{\hat{n} \times \hat{f} \hat{f}}^{-1} P_k T_k P_k \\ & \qquad + \xi (M_{ik}^{-} G_{\hat{n} \times \hat{f} \hat{f}}^{-1} M_k^{+} - P_{\Sigma H} M_0^{-} G_{\hat{n} \times \hat{f} \hat{f}}^{-1} M_0^{+} P_{\Lambda H})) P_k y \\ &= -\tilde{P}_k (\tilde{T}_{ik} \tilde{P}_k G_{\hat{n} \times \hat{f} \hat{f}}^{-1} P_k e^i - \xi M_{ik}^{-} G_{\hat{n} \times \hat{f} \hat{f}}^{-1} h^i), \end{aligned} \quad (142) $$

with a suitably chosen $\xi$ and in high-frequency scenarios, the scaling of the projectors must be set to a unitary value [96]. Figure 12 shows that the CFIE formulation is free from interior resonances and Figure 13 the number of iterations as a function of the spectral index $1/h$. For the presented results, we have weighted the EFIE and MFIE part equally. Finally, this combined formulation has the advantage of not exhibiting a static nullspace. This is because, unlike the standard the Calderón EFIE (126), the Calderón EFIE stabilized with projectors (128) has no static nullspace, which can additively cure that of the MFIE and results in equation (142) that is free from spurious resonances, as illustrated in Figure 14.

While both Calderón-Yukawa CFIEs (141) and (142) are free from the dense-discretization breakdown (i.e., when all parameters are kept constant expect from $h$, the condition number can be asymptotically bounded for $h \to 0$), the condition number cannot be bounded independently from $k$. This effect has been designated as the high-frequency breakdown.

### B. HIGH-FREQUENCY BREAKDOWN OF THE CALDERÓN-YUKAWA CFIE

The high-frequency breakdown corresponds to the regime in which the wavenumber $k$ increases and the mesh parameter $h$ decreases proportionally to the inverse of the wavenumber. Concretely, it corresponds to increasing the frequency while keeping the number of basis functions per wavelength constant (e.g., ten subdivisions per wavelength is a common choice). In other words, we are interested in the high-frequency regime which is characterized by $hk = \text{const}$ while $k \to \infty$.

It can be shown that in general both the standard CFIE and the Calderón-Yukawa CFIE have singular values that

<footnote n="1">1. Operator preconditioning, in this logic, strives to find $B$ as an approximation of $A^{-1}$.</footnote>
<figure><img src="image_16.png" /><figcaption>FIGURE 14. Low frequency (10⁻¹² Hz) spectra of several formulations putting into evidence their magnetostatic nullspace, or lack thereof. The matrices have been obtained by discretizing a torus with a square the cross section of 0.2 m sides and with a minimum inner diameter of 2 m. The average edge length of the mesh is of 0.2 m and the spectra have been normalized to make the figure more readable. The combination parameter ξ = i was used here.</figcaption></figure>

are unbounded in $k$. As, however, the problem is that the singular values of the analytic operator are unbounded in $k$, the condition number of the Calderón-Yukawa CFIE grows when $k \to \infty$. On spherical geometries, it has been proved that the condition number grows as $O(k^{2/3})$ [162]. As was seen in the previous sections, the spectrum of the electromagnetic boundary integral operators can be separated into parts: one for the solenoidal functions and the other for the non-solenoidal functions. The $O(k^{2/3})$ growth of the condition number is due to a maximum singular value in the solenoidal part of the spectrum that grows as $O(k^{1/3})$ and a minimum singular value in the non-solenoidal part of the spectrum that decreases as $O(k^{-1/3})$. This effect is related to the high-frequency breakdown in the scalar Helmholtz equation in acoustics, where both the Dirichlet problem solved with a standard CFIE and the Neumann problem solved with a Calderón-Yukawa CFIE have condition numbers that grow as $O(k^{1/3})$ on a sphere [163], [164]. A more detailed analysis of the high frequency breakdown on the sphere is given in the Appendix.

### C. REMEDIES

The preconditioners to cure the high-frequency breakdown that are found in the literature are based on the use of a modified wavenumber $k_m$. This modified wavenumber $k_m$ is the original wavenumber $k$ plus a shift in the complex plane. This modified wavenumber has the form [165]

$$k_m = k + i c k^{1/3} a^{-2/3}, \quad (143)$$

where $c$ is a dimensionless constant and $a$ is a length that depends on the geometry. For example, on spherical geometries, $a$ is the radius of the sphere and $c$ is often chosen to be $c \approx 0.4$, which is actually based on optimizations for the conditioning of the acoustic operators [165], but also give good

results in the electromagnetic case. A heuristic from [162] is to choose as the maximum of the absolute value of the mean curvature of the object boundary $H = 1/a$.

Among the solutions found in [162] proposes a Calderón-like preconditioning with a modified wavenumber $k_m$ in the preconditioner of the EFIE operator, then an unchanged MFIE is added in a CFIE fashion. The work [166] also uses a Calderón-like preconditioning of the MFIE with the same $k_m$. The work in [123] uses the inverse square root of the vector Helmholtz operator with the modified wavenumber as a preconditioner for the non-solenoidal part of the equation. In [167], two different symmetric formulations for the EFIE are proposed that are high-frequency stable (in addition to being low-frequency and dense-discretization stable): one based on a loop-star decomposition with a scalar Helmholtz operator (with $k_m$) and its inverse to precondition separately the loop and star components, and the other based on quasi-Helmholtz projectors to precondition independently the solenoidal and the non-solenoidal part with the vector Helmholtz operator and its inverse. Essentially, in the two cases, the EFIE is squared with the Helmholtz operator in the middle, which avoids the discretization of the EFIE with dual basis functions similar to the refinement-free Calderón preconditioner from Section V-A2.

## VII. A NOTE ON THE ROLE OF THE BASIS FUNCTION GRAM MATRIX ON CONDITIONING

Given that several of the techniques described above aim to obtain a preconditioned equation which is spectrally equivalent to an identity, it follows that the final condition number is related to the one of the (potentially mixed) Gram matrix between the rightmost expansion and leftmost testing basis functions used. In the special case in which these two bases coincide, the Gram matrix of the basis will be the one dictating ultimately the conditioning behavior. There are two main causes of ill-conditioning for a Gram matrix associated with a specific basis function set: (i) the mesh quality is low or (ii) the basis functions linear dependency increases with one of the relevant parameters (mesh size and polynomial order, for example). Below we will briefly discuss these two cases together with related solution approaches.

### A. ILL-SHAPED MESH ELEMENTS

Simply speaking, a mesh is ill-shaped if it has narrow triangles (see Figure 15 for a mesh with narrow triangles; for a more nuanced discussion on measures for the quality of meshes, see [168]). Ill-shaped meshes lead to higher condition numbers of the discretized operator compared with a discretization based on a well-shaped mesh, that is, where the angles of each triangle have roughly the same size. To quantify the impact of ill-conditioning due to the mesh, one can consider the condition number of the discretized $I$, that is, $G_{ff}$: the condition number will be a baseline for any other integral operator discretized with RWG basis functions. Clearly, this problem is not specific to surface integral
<figure><img src="image_17.png" /><figcaption>FIGURE 15. Example of ill-shaped triangles.</figcaption></figure>

equation methods, volume integral equation methods or finite element methods also suffer from this issue.

This kind of Gram matrix ill-conditioning can be cured by the optimization of the geometry or the application of mesh smoothers, a strategy to avoid this cumbersome (often manual) work has been presented by Stephenanson and Lee [169]. Based on a study of the SVD of the EFIE system matrix, degrees of freedom associated with the ill-shaped triangles are eliminated. An alternative strategy can be the use of discontinuous Galerkin (DG) strategies. In fact, one source of ill-shaped meshes stems from geometries that necessitate a finer discretization of certain parts of the geometry to either capture important details or to avoid a work intensive cleanup of the geometry. With DG individual parts of the surface are meshed independently from each other. Weak conditions are used to enforce current continuity [37]. Lastly, one can entirely avoid this issue by discretizing directly on the spline-based geometry, that is, leveraging methods from isogeometric analysis. For the electric field integral equation such approaches have recently been presented [170]–[173].

Improvements of the Calderón preconditioning in (128) can be expected by placing an additional Gram matrix resulting in

$$ 
\begin{aligned}
G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-T} \tilde{P}_{\tilde{k}} \tilde{T}_{\tilde{k}} \tilde{P}_{\tilde{k}} G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-1} P_{\tilde{k}} T_{\tilde{k}} P_{\tilde{k}} Y \\
= -G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-T} \tilde{P}_{\tilde{k}} \tilde{T}_{\tilde{k}} \tilde{P}_{\tilde{k}} G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-1} \tilde{P}_{\tilde{k}} e^{i},
\end{aligned}
\quad (144) $$

where $G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-T} = (G_{\tilde{n} \times \tilde{f} \tilde{f}}^{T})^{-1}$. It ensures that the condition number can be bounded independently from $h$ and the bases used for discretization. For a well-conditioned bases such as RWG and BC functions, this is not necessary. However, it typically improves the condition number, in particular, if the mesh contains ill-shaped triangles or if the discretization is non-uniform which occurs, for instance, in the case of meshes resulting from adaptive refinement. The cost for inverting $G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-T}$ is often outweighed by a saving in the number of iterations required to solve the preconditioned system since the cost of a matrix-vector product with $T$ or $\tilde{T}$ is often more expensive than inverting $G_{\tilde{n} \times \tilde{f} \tilde{f}}^{-T}$ iteratively.

## B. BASIS FUNCTIONS LINEAR DEPENDENCY

The choice of basis functions affects the conditioning of the Gram matrix and thus of the discretized operator (especially when the latter is of the second kind, as it often happens for preconditioned formulations). A particularly

important case of Gram matrix ill-conditioning though is the one found in higher-order basis functions. The most popular types of higher-order bases for electromagnetic integral equations are interpolatory and hierarchical higher-order bases [174]–[185].

Hierarchical higher-order basis functions should not be confused with hierarchical bases used for curing the dense-discretization breakdown, which we discussed in Section V-B. Typical hierarchical bases for curing the dense-discretization breakdown of the EFIE are of the same polynomial order as first order RWG bases. While the use of such hierarchical basis functions leads with an increasing order $p$ to a faster convergence with respect to the analytic solution, it also increases the condition number. For interpolatory bases, an exponential growth of the condition number in $p$ has been reported, in contrast to a polynomial growth for hierarchical bases [180, p. 187].

Strategies to overcome the $p$-breakdown rely on constructing orthogonal higher-order basis functions. In fact, hierarchical higher-order basis functions enforce orthogonality between functions on the same level [177], [178]. Further improvements have been presented for wire-, quadrilateral-, and brick-type elements [187]–[189].

For some of these hierarchical higher-order basis functions, Calderón-type preconditioners have been presented [189]. As a consequence from operator preconditioning theory [110], [190], it is possible to obtain a higher-order discretization of the EFIE, where the condition number is bounded independently from $h$ and $p$. It must be noted though that while this is true for the overall condition number, the (inverse) Gram matrices that appear in the Calderón preconditioner remain ill-conditioned in $p$.

## VIII. CONCLUSION

This paper has presented an overview of the state of the art in electromagnetic integral equation preconditioning. Different spectral regularization strategies have been discussed, further analyzed, and framed in the overall context of solving challenging large computational electromagnetics problems. Bibliographic reviews have been alternated with theoretical treatments and corroborating numerical results to show the practical impacts of all the strategies.

## APPENDIX A
CONTINUOUS INTERPRETATIONS

In this section, $\Gamma$ is a sphere of radius $a$ centered at the origin. The scalar spherical harmonics (SH) are noted $Y_{lm}$ with $l \ge 0$ and $|m| < l$. The vector spherical harmonics (VSH) are defined as

$$ Y_{lm}(\hat{r}) = \hat{r} Y_{lm}(\hat{r}), \quad (145) $$

$$ X_{lm}(\hat{r}) = \frac{a}{i\sqrt{l(l+1)}} \hat{r} \times \nabla Y_{lm}(\hat{r}), \quad (146) $$

$$ U_{lm}(\hat{r}) = -\frac{a}{i\sqrt{l(l+1)}} \nabla Y_{lm}(\hat{r}), \quad (147) $$
They form an orthonormal basis for vector fields on $\Gamma$ with respect to the weighted inner product

$$ \langle \mathbf{u}, \mathbf{v} \rangle = \int_{\Gamma} \bar{\mathbf{u}} \cdot \mathbf{v} \frac{dS(\mathbf{r})}{a^2} = \int_{\varphi=0}^{2\pi} \int_{\theta=0}^{\pi} \bar{\mathbf{u}} \cdot \mathbf{v} \sin\theta \, d\theta \, d\varphi. \quad (148) $$

Note that $X_{lm}$ and $U_{lm}$ are tangential, while $Y_{lm}$ are radial vector fields. Therefore, only $X_{lm}$ and $U_{lm}$ are used in the analysis of the surface integral operators. By construction, this basis is Helmholtz decomposed: $X_{lm}$ are solenoidal and $U_{lm}$ are irrotational.

In the following, $J_l$ is the Riccati-Bessel function and $H_l$ is the Riccati-Hankel function of first kind. They are related to the Bessel function $J_l$ and the Hankel function of first kind $H_l$ by $J_l(z) = \sqrt{\pi z/2} J_{l+1/2}(z)$ and $H_l(z) = \sqrt{\pi z/2} H_{l+1/2}(z)$. Also, $\ell$ denotes their derivative with respect to their argument. We thus obtain the analytic expressions of the EFIE operator [108], [191], [192]

$$ \mathcal{T}_k X_{lm} = -J_l(ka) H_l(ka) U_{lm}, \quad (149) $$

$$ \mathcal{T}_k U_{lm} = J'_l(ka) H'_l(ka) X_{lm}, \quad (150) $$

the exterior MFIE operator

$$ (\mathcal{I}/2 + \mathcal{K}_k)X_{lm} = i J'_l(ka) H_l(ka) X_{lm}, \quad (151) $$

$$ (\mathcal{I}/2 + \mathcal{K}_k)U_{lm} = -i J_l(ka) H'_l(ka) U_{lm}, \quad (152) $$

and the interior MFIE operator

$$ (-\mathcal{I}/2 + \mathcal{K}_k)X_{lm} = i J_l(ka) H'_l(ka) X_{lm}, \quad (153) $$

$$ (-\mathcal{I}/2 + \mathcal{K}_k)U_{lm} = -i J'_l(ka) H_l(ka) U_{lm}. \quad (154) $$

The tangential vector fields $X_{lm}$ and $U_{lm}$ are eigenvectors of the MFIE operator, and they are singular vectors of the EFIE operators. Assume we discretize the operator $\mathcal{T}$ with the VSH, then it is clear that the singular values associated with the right singular vectors $X_{lm}$ are $\Sigma_{l,X} = |J_l(ka) H_l(ka)|$, and the singular values associated with the right singular vectors $U_{lm}$ are $\Sigma_{l,U} = |J'_l(ka) H'_l(ka)|$. By studying the behavior of these singular values in different regimes, we can obtain the asymptotic behavior of the condition number of the discretized EFIE. Since $\nabla \cdot X_{lm} = 0$, the study of the operators applied to $X_{lm}$ gives the behavior of the operators applied to solenoidal functions. Whereas the study of the operators applied to $U_{lm}$ gives the behavior of the operators applied to non-solenoidal functions. The following asymptotic expressions can be obtained easily from the asymptotic behavior of the special functions in different regimes [193].

# A. LOW-FREQUENCY BREAKDOWN AND OTHER ISSUES

## 1) LOW-FREQUENCY BREAKDOWN

We study the condition number of the discretized EFIE for $k \to 0$, where we keep the number of basis functions constant. In this case, it corresponds to considering only the orders $l$ less than some constant $l_{max}$. Thus in the asymptotic expansion $l$ is fixed and independent of $k$. In the following,

<figure><img src="image_18.png" /><figcaption>FIGURE 16. Singular values of the EFIE associated to the first three harmonic orders ($l = 1, 2, 3$) as functions of the frequency. The low frequency behavior is given by a growth of the maximum singular value as $(ka)^{-1}$, and the decrease of minimum singular value as $(ka)$, which results in a condition number that grows proportionally to $(ka)^{-2}$ as $ka \to 0$. In higher frequencies there are the first resonant frequencies (at $ka \approx 2.74, 3.87$, etc.) around which the condition number is unbounded.</figcaption></figure>

we will use the asymptotic behavior of the Riccati-Bessel, Riccati-Hankel, and their derivatives for small arguments given by

$$ J_l(ka) \underset{k \to 0}{\sim} \frac{(ka)^{l+1}}{(2l+1)!!}, \quad (155) $$

$$ J'_l(ka) \underset{k \to 0}{\sim} \frac{(l+1)(ka)^l}{(2l+1)!!}, \quad (156) $$

$$ H_l(ka) \underset{k \to 0}{\sim} \frac{i(2l-1)!!}{(ka)^l}, \quad (157) $$

$$ H'_l(ka) \underset{k \to 0}{\sim} \frac{i l (2l-1)!!}{(ka)^{l+1}}, \quad (158) $$

where !! is the double factorial (or semifactorial) that is defined by the recurrence relation $n!! = (n-2)!!n$ and $1!! = 0!! = 1$. Therefore, as $k \to 0$, the singular values scale as

$$ \sigma_{l,X} = |J_l(ka) H_l(ka)| = O(k), \quad (159) $$

$$ \sigma_{l,U} = |J'_l(ka) H'_l(ka)| = O(k^{-1}). \quad (160) $$

As a consequence, the condition number grows as $O(k^2)$. The singular values of the EFIE operator on a sphere have been represented on Figure 16.

A way to solve this problem is by scaling the solenoidal part by $k^{-1}$ and by scaling the non-solenoidal part by $k$. In the VSH basis, it is simply a diagonal preconditioner, but in practical applications one has to use some techniques such as a loop-star decomposition or quasi-Helmholtz projectors to decompose the two subspaces before rescaling them. Another way to solve the problem is by Calderón preconditioning. In this case, the main operator is $(\mathcal{T}_k)^2$, whose singular values
are

$$|\mathcal{J}_l(ka) \mathbb{H}_l(ka) \mathcal{J}'_l(ka) \mathbb{H}'_l(ka)| = O(1) \quad (161)$$

as $k \to 0$. As a result, the condition number of the Calderón EFIE remains bounded as $k \to 0$.

## 2) LOW-FREQUENCY ROUND-OFF ERRORS IN THE EXCITATION

We consider the plane wave $f_{pw}(\mathbf{r}) = e_0 e^{i\hat{k}\cdot\hat{r}} \hat{\mathbf{x}}$. Its tangential component expands as [194]

$$\langle \mathbf{X}_{l, \pm 1}, f_{\text{pw}} \rangle = 2\pi i^l e_0 \sqrt{\frac{2l+1}{4\pi}} \frac{\mathcal{J}_l(ka)}{ka}, \quad (162)$$

$$\langle \mathbf{U}_{l, \pm 1}, f_{\text{pw}} \rangle = \pm 2\pi i^l e_0 \sqrt{\frac{2l+1}{4\pi}} \frac{\mathcal{J}'_l(ka)}{ka}. \quad (163)$$

For $k \to 0$, the dominant terms of the solenoidal and non-solenoidal parts are

$$\langle \mathbf{X}_{1, \pm 1}, f_{\text{pw}} \rangle = O(ka), \quad (164)$$

$$\langle \mathbf{U}_{1, \pm 1}, f_{\text{pw}} \rangle = O(1). \quad (165)$$

When the RWG basis functions are used for the discretization, the two components are summed together; due to finite machine precision and the associated numerical round-off errors, the solenoidal part is lost.

## 3) LOW-FREQUENCY ROUND-OFF ERRORS IN THE SOLUTION

Consider the EFIE $\mathcal{T}_k j_\Gamma = -\hat{\mathbf{n}} \times \mathbf{e}^i$ when the excitation is the plane wave $\mathbf{e}^i(\mathbf{r}) = e_0 e^{i\hat{k}\cdot\hat{r}} \hat{\mathbf{x}}$. From (149) and (162), the solution $j_\Gamma$ has the following expansion

$$\langle \mathbf{X}_{l, \pm 1}, j_\Gamma \rangle = \frac{2\pi i^l}{\mathbb{H}_l(ka)ka} \sqrt{\frac{2l+1}{4\pi}} e_0, \quad (166)$$

$$\langle \mathbf{U}_{l, \pm 1}, j_\Gamma \rangle = \pm \frac{2\pi i^l}{\mathbb{H}'_l(ka)ka} \sqrt{\frac{2l+1}{4\pi}} e_0. \quad (167)$$

For $k \to 0$, the dominant terms of the solenoidal and non-solenoidal parts are respectively

$$\langle \mathbf{X}_{1, \pm 1}, j_\Gamma \rangle = O(1), \quad (168)$$

$$\langle \mathbf{U}_{1, \pm 1}, j_\Gamma \rangle = O(ka). \quad (169)$$

Again, similarly to the round-off errors in the excitation, the solution has two Helmholtz components that scale differently with the frequency. Therefore, in a numerical solver, unless some care is taken by treating the two components separately, the non-solenoidal part will get canceled.

## B. DENSE-DISCRETIZATION BREAKDOWN

We study the condition number of the discretized EFIE when $h$ decreases and the frequency remains constant. A decrease in $h$ is equivalent to an increase of the maximum spacial frequency of the functions that are discretized. In the case of the VSH, their spacial frequency increases with their order $l$. The number of VSH whose order is

less or equal to $l$ grows quadratically with $l$, due to their multiplicity on the order $m$ ($|m| \le l$). When we use local basis functions such as RWGs to discretize the operators, the number of unknowns also grows quadratically with $a/h$. Therefore, to study the dense-discretization breakdown, we can use the correspondence $l = O(a/h)$. In the following, we will use the asymptotic behavior of the Riccati-Bessel, Riccati-Hankel, and their derivatives for large order given by

$$\mathcal{J}_l(ka) \underset{l \to +\infty}{\sim} \frac{1}{\sqrt{2e}} \left( \frac{e ka}{2l+1} \right)^{l+1}, \quad (170)$$

$$\mathcal{J}'_l(ka) \underset{l \to +\infty}{\sim} \sqrt{\frac{e l+1}{2l+1}} \left( \frac{e ka}{2l+1} \right)^l, \quad (171)$$

$$\mathbb{H}_l(ka) \underset{l \to +\infty}{\sim} -i \sqrt{\frac{2}{e}} \left( \frac{2l+1}{e ka} \right)^l, \quad (172)$$

$$\mathbb{H}'_l(ka) \underset{l \to +\infty}{\sim} \frac{i l \sqrt{2e}}{2l+1} \left( \frac{2l+1}{e ka} \right)^{l+1}. \quad (173)$$

Therefore, as $l \to +\infty$, the singular values scale as

$$\sigma_{l,X} = |\mathcal{J}_l(ka) \mathbb{H}_l(ka)| = O(l^{-1}) = O(h), \quad (174)$$

$$\sigma_{l,U} = |\mathcal{J}'_l(ka) \mathbb{H}'_l(ka)| = O(l) = O(h^{-1}). \quad (175)$$

As a consequence, the condition number grows as $O(h^{-2})$. The spectrum of the EFIE have been represented on Figure 17 for different frequencies to illustrate this growth of the condition number as the discretization density increases.

This problem can be solved using Calderón preconditioning: as $l \to +\infty$ (or equivalently as $h \to 0$), the singular values of $(\mathcal{T}_k)^2$ are

$$|\mathcal{J}_l(ka) \mathbb{H}_l(ka) \mathcal{J}'_l(ka) \mathbb{H}'_l(ka)| = \sigma_{l,X} \sigma_{l,U} = O(1). \quad (176)$$

As a result the condition number of the Calderón EFIE remains bounded as $h \to 0$.

## C. HIGH-FREQUENCY RESONANCES

We study the condition number for $ka \gg 1$. Here the point to notice is that $l$ is kept constant as $ka \to +\infty$, and

$$\mathcal{J}_l(ka) = \sin(ka - l\pi/2) + O((ka)^{-1}), \quad (177)$$

$$\mathcal{J}'_l(ka) = \cos(ka - l\pi/2) + O((ka)^{-1}). \quad (178)$$

So the Riccati-Bessel function and its derivative have zeros, and as a consequence, there are an infinite number of frequencies such that one of the $\sigma_{l,X}$ or $\sigma_{l,U}$ is 0. Therefore the condition number of the EFIE is unbounded around these resonant frequencies, which makes the EFIE unsolvable in practice without further treatment. The MFIE suffers from the same problem.

The resonant frequencies can be read on Figure 16 with the first one appearing at $ka \approx 2.74$. Also, it is clear from Figure 16 that the higher $ka$ is, the more common the resonant frequencies are. The result of the resonances on the spectrum is clear on Figure 17, where there is a finite number
<figure><img src="image_19.png" /><figcaption>FIGURE 17. Singular values of the operator of the EFIE on a sphere for three different frequencies. The maximum singular value grows as h⁻¹ and the minimum singular value decreases as h, resulting in a condition number growth proportional to h⁻².</figcaption></figure>

of singular values that are unbounded from below at a fixed frequency.

The common way to solve this problem is to discretize the CFIE, which combines the EFIE (6) scaled by α and the MFIE (11) on which is applied (1 − α) ˆn ×, where 0 < α < 1 is a dimensionless constant. The CFIE is

$$
\begin{align}
\alpha \hat{\mathbf{n}} \times \mathbf{e}^i + (1 - \alpha) \hat{\mathbf{n}} \times \hat{\mathbf{n}} \times \mathbf{h}^i & \nonumber \\
& = (-\alpha \mathcal{T}_k + (1 - \alpha) \hat{\mathbf{n}} \times (\mathcal{I}/2 + \mathcal{K}_k)) \mathbf{j}_\Gamma. \tag{179}
\end{align}
$$

The singular values of this CFIE operator associated with
the right singular vector $X_{lm}$ and $U_{lm}$ are, respectively,

$$
|\alpha \mathcal{J}_l(ka) + i(1-\alpha) \mathcal{J}'_l(ka)| |\mathcal{H}_l(ka)|, \quad (180)
$$

$$
|\alpha \mathcal{J}'_l(ka) - i(1-\alpha) \mathcal{J}_l(ka)| |\mathcal{H}'_l(ka)|, \quad (181)
$$

which are never 0 since the zeros of J_l and J'_l never coin-
cide (due to the fact that the zeros of the Bessel functions
are simple with an exception at ka = 0, which we do
not consider). In other words, the CFIE can always be
inverted.

Note that the CFIE only fixes the problem of the res-
onances that make the MFIE operator and EFIE operator
non-invertible for some frequencies. Around these resonant
frequencies, the condition number is unbounded. However,
the CFIE still suffer from the dense-mesh breakdown.
The Calderón-Yukawa CFIE is dense-mesh stable, low-
frequency stable and free from resonances. But it still
suffers from the high-frequency breakdown, which is an
asymptotic growing of the condition number as k → ∞
while keeping kh constant. It is discussed in the next
section.

<figure><img src="image_20.png" /><figcaption>FIGURE 18. Singular values of the operator of the Calderón-Yukawa CFIE (182) on a sphere for three different frequencies. The maximum singular value grows as (ka)¹/³ and the minimum singular value decreases as (ka)⁻¹/³, resulting in a condition number growth proportional to (ka)²/³.</figcaption></figure>

D. HIGH-FREQUENCY BREAKDOWN

In this section we study the asymptotic behavior of the
singular values of the Calderón-Yukawa CFIE operator [96]

$$
- \mathcal{T}_{ik} \mathcal{T}_k + (\mathcal{I}/2 - \mathcal{K}_{ik})(\mathcal{I}/2 + \mathcal{K}_k). \quad (182)
$$

Recall that the CFIE is used to avoid the interior reso-
nances, while Calderón preconditioning avoids the dense-
discretization breakdown. The Yukawa part denotes the use
of pure imaginary wavenumbers in the preconditioning oper-
ators so that the preconditioners do not reintroduce interior
resonances.

We show in this section that the maximum singular value,
which is associated with the solenoidal functions, scales at
least proportionally to k¹/³, while the minimum singular
value, which is associated with the non-solenoidal functions,
scales at most as O(k⁻¹/³). The result is a growth of the
condition number faster than a constant times k²/³, where
numerical evidences on the spherical harmonics show that
this limit is in practice attained. The singular values of the
Calderón-Yukawa CFIE have been represented on Figure 18
as well as the asymptotic bounds that are derived in the
following sections.

1) *Minimum singular value:* The singular value associated with the nonsolenoidal VSH $U_{lm}$ are

$$
\sigma_{l,U} = |\mathcal{H}_l(i ka) \mathcal{H}'_l(ka) (\mathcal{J}'_l(ika) \mathcal{J}_l(ka) + \mathcal{J}_l(ika) \mathcal{J}'_l(ka))|. \quad (183)
$$

It is useful to treat the order *l* as a real variable instead
of an integer. We obtain an upper bound for the mini-
mum singular values by choosing as *l* the value for which
$J_l(ka) = \sqrt{\pi ka/2} J_{l+1/2}(ka) = 0$, that is, for which $ka$ is a zero of $J_{l+1/2}$ (J is the Bessel function and $J$ is the Ricatti-Bessel function). For the first zero of $J_{l+1/2}$, the asymptotic expansion is given by [193, eq. (10.21.40)]

$$ka \sim l + \alpha l^{1/3} + 1/2, \quad (184)$$

or, equivalently, when

$$l \sim ka - \alpha(ka)^{1/3} - 1/2 \quad (185)$$

with $\alpha = -2^{-1/3}a_1 \approx 1.856$, where $a_1 \approx -2.34$ is the first zero of Ai, which denotes the Airy function of the first kind. Using the uniform asymptotic expansion of the derivatives of the Bessel and Hankel functions for large orders [193, eq. (10.20.7-9)], we find that the upper bound for the minimum singular value is asymptotically

$$\frac{\pi \operatorname{Ai}'(a_1) |\operatorname{Ai}'(a_1) - i \operatorname{Bi}'(a_1)|}{2^{7/6}} (ka)^{-1/3} \approx 0.69(ka)^{-1/3}, \quad (186)$$

where Bi denotes the Airy function of the second kind.

2) **Maximum singular value:** The singular values associated with the solenoidal VSH $X_{lm}$ are

$$\sigma_{l,X} = |\mathbb{H}'_l(ika) \mathbb{H}_l(ka)(\mathbb{J}'_l(ika) \mathbb{J}_l(ka) + \mathbb{J}_l(ika) \mathbb{J}'_l(ka))|.$$

Similarly to the previous case, there is another constant $\tilde{\alpha}$ such that the index of a lower bound for the maximum singular value is asymptotically

$$l \sim ka - \tilde{\alpha}(ka)^{1/3} - 1/2. \quad (187)$$

Using the asymptotic expansion of the Bessel functions in the transition region [193, eq. (10.19.8-9)] we find that the dominant term in $\sigma_{l,X}$ is proportional to $Ai(-2^{1/3}\tilde{\alpha})|Ai(e^{-i\pi/3}2^{1/3}\tilde{\alpha})|$. So $\tilde{\alpha}$ has to maximize this quantity, and a derivative of it shows that $\tilde{\alpha}$ is solution of

$$\operatorname{Ai}'(z) (2 \operatorname{Ai}(z)^2 + \operatorname{Bi}(z)^2) + \operatorname{Ai}(z) \operatorname{Bi}(z) \operatorname{Bi}'(z) = 0, \quad (188)$$

where $z = -2^{1/3}\tilde{\alpha}$. Numerically, we obtain $\tilde{\alpha} \approx 0.623$, and the lower bound for the maximum singular value is asymptotically

$$\frac{\pi \operatorname{Ai}(z) |\operatorname{Ai}(z) - i \operatorname{Bi}(z)|}{2^{5/6}} (ka)^{1/3} \approx 0.524(ka)^{1/3}. \quad (189)$$

3) **High-frequency stabilization:** Notice that the Calderón-EFIE (and the MFIE, respectively) does, in some sense, not suffer from the high-frequency breakdown since the singular values of its operator are bounded from above by a constant independent from $k$, but it is unstable due to the interior resonances. In contrast, the Calderón-Yukawa CFIE is free from interior resonances, but it suffers from the high-frequency breakdown as the maximum singular values of its operator is unbounded in $k$. To solve this issue one can construct a Calderón-CFIE, but instead of using the wavenumber $k$ in the preconditioner, one uses a modified wavenumber $k_m$ that

<figure><img src="image_21.png" /><figcaption>FIGURE 19. Singular values of the Calderón-CFIE with modified wavenumber $k_m$ in the preconditioning operators (190) on a sphere for three different frequencies. The singular values are bounded from above and from below by a constant resulting in a stable condition number when the frequency increases.</figcaption></figure>

is shifted in the complex plane to avoid reintroducing the resonances resulting in [166]

$$-\mathcal{T}_{k_m} \mathcal{T}_k + (\mathcal{I}/2 - \mathcal{K}_{k_m}) (\mathcal{I}/2 + \mathcal{K}_k). \quad (190)$$

This modified wavenumber for the sphere has the form [123], [162]

$$k_m = k + ick^{1/3}a^{-2/3}, \quad (191)$$

where $c \in \mathbb{R}$ is a constant. A common choice that gives a good condition number is $c \approx 0.4$. The singular values of this operator are similar to those of the Calderón-Yukawa CFIE with the pure imaginary wavenumber $ik$ replaced by $k_m$, namely

$$\sigma_{l,U} = |\mathbb{H}_l(k_m a) \mathbb{H}'_l(ka) (\mathbb{J}'_l(k_m a) \mathbb{J}_l(ka) + \mathbb{J}_l(k_m a) \mathbb{J}'_l(ka))| \quad (192)$$

$$\sigma_{l,X} = |\mathbb{H}'_l(k_m a) \mathbb{H}_l(ka) (\mathbb{J}'_l(k_m a) \mathbb{J}_l(ka) + \mathbb{J}_l(k_m a) \mathbb{J}'_l(ka))|. \quad (193)$$

These singular values have been plotted in Figure 19 for increasing frequencies. It is evident that the condition number is bounded.

## APPENDIX B
### INVERSE GRAM MATRIX PROPERTIES

**Lemma 1:** For the inverse mixed Gram matrix $G_{\hat{n}\times\tilde{f},f}$, we have

$$P_{\Sigma H} G_{\hat{n}\times\tilde{f},f}^{-1} P_{\Lambda} = \mathbf{0}, \quad (194)$$

$$P_{\Sigma} G_{\hat{n}\times\tilde{f},f}^{-1} P_{\Lambda H} = \mathbf{0}. \quad (195)$$

*Proof:* First, we define the transformation matrix $Q = [\Lambda H \Sigma]$, where the global loops are constructed such
$$ (Q^T G Q)^{-1} = \begin{bmatrix} (\Lambda^T G \Lambda)^{-1} - (\Lambda^T G \Lambda)^{-1} \Lambda^T G H (H^T G H)^{-1} & (\Lambda^T G \Lambda)^{-1} \Lambda^T G H (H^T G H)^{-1} H^T G \Sigma (\Sigma^T G \Sigma)^{-1} - (\Lambda^T G \Lambda)^{-1} \Lambda^T G \Sigma (\Sigma^T G \Sigma)^{-1} \\ \mathbf{0} & (H^T G H)^{-1} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} - (H^T G H)^{-1} H^T G \Sigma (\Sigma^T G \Sigma)^{-1} (\Sigma^T G \Sigma)^{-1} \quad (200) $$

$$ Q(Q^T G Q)^{-1} Q^T = \Lambda(\Lambda^T G \Lambda)^{-1} \Lambda^T - \Lambda(\Lambda^T G \Lambda)^{-1} \Lambda^T G H (H^T G H)^{-1} H^T + \Lambda(\Lambda^T G \Lambda)^{-1} \Lambda^T G H (H^T G H)^{-1} H^T G \Sigma (\Sigma^T G \Sigma)^{-1} \Sigma^T \\ - \Lambda(\Lambda^T G \Lambda)^{-1} \Lambda^T G \Sigma (\Sigma^T G \Sigma)^{-1} \Sigma^T + H(H^T G H)^{-1} H^T - H(H^T G H)^{-1} H^T G \Sigma (\Sigma^T G \Sigma)^{-1} \Sigma^T + \Sigma (\Sigma^T G \Sigma)^{-1} \Sigma^T \quad (201) $$

that $\Lambda^T H = \mathbf{0}$. Furthermore, without loss of generality, we assume that loop and star functions have been eliminated such that they are linearly independent. In this case, we can denote the projectors as

$$ P_{\Lambda H} = \Lambda(\Lambda^T \Lambda)^{-1} \Lambda^T + H(H^T H)^{-1} H^T, \quad (196) $$

$$ P_{\Sigma H} = \Sigma (\Sigma^T \Sigma)^{-1} \Sigma^T + H(H^T H)^{-1} H^T. \quad (197) $$

To reveal how $G_{\hat{n} \times \tilde{f} f}^{-1}$ acts on the different quasi-Helmholtz subspaces, we derive an explicit expression for $(Q^{-T}Q^T G Q Q^{-1})^{-1} = Q(Q^T G Q)^{-1}Q^T$, where in the following we abbreviate $G_{\hat{n} \times \tilde{f} f} = G$. First, we recall the block structure of

$$ Q^T G Q = \begin{bmatrix} \Lambda^T G \Lambda & \Lambda^T G H & \Lambda^T G \Sigma \\ H^T G \Lambda & H^T G H & H^T G \Sigma \\ \Sigma^T G \Lambda & \Sigma^T G H & \Sigma^T G \Sigma \end{bmatrix} = \begin{bmatrix} \Lambda^T G \Lambda & \Lambda^T G H & \Lambda^T G \Sigma \\ \mathbf{0} & H^T G H & H^T G \Sigma \\ \mathbf{0} & \mathbf{0} & \Sigma^T G \Sigma \end{bmatrix}, \quad (198) $$

which follows from the orthogonality of solenoidal (i.e., $\Lambda$) and irrotational functions (i.e., $\hat{n} \times \tilde{\Sigma}$) / harmonic functions with irrotational perturbation (i.e., $\hat{n} \times \tilde{H}$).

Using the Schur complement, we have

$$ \begin{bmatrix} A & B \\ \mathbf{0} & D \end{bmatrix}^{-1} = \begin{bmatrix} A^{-1} & -A^{-1}BD^{-1} \\ \mathbf{0} & D^{-1} \end{bmatrix}, \quad (199) $$

Recursively application of (199) to (198) yields (200) shown at the top of the page; explicitly computing $Q(Q^T G Q)^{-1}Q^T$ results in (201), shown at the top of the page. Evidently, we have

$$ H^T G^{-1} \Lambda = \mathbf{0}, \quad (202) $$

$$ \Sigma^T G^{-1} \Lambda = \mathbf{0}, \quad (203) $$

$$ \Sigma^T G^{-1} H = \mathbf{0}, \quad (204) $$

and considering the definition of projectors, the lemma follows. ■

## REFERENCES

[1] A. F. Peterson, S. L. Ray, and R. Mittra, *Computational Methods for Electromagnetics*. New York, NY, USA: IEEE Press/Oxford Univ. Press. 1998.

[2] J.-M. Jin, *Theory and Computation of Electromagnetic Fields*, 2nd ed. New York, NY, USA: Wiley, 2015, p. 715.

[3] W. C. Gibson, *The Method of Moments in Electromagnetics*, 2nd ed. Boca Raton, FL, USA: CRC Press/Taylor Francis, 2014, p. 425.

[4] D. Colton and R. Kress, *Integral Equation Methods in Scattering Theory*. Philadelphia, PA, USA: SIAM, 2013.

[5] P. Monk, *Finite Element Methods for Maxwell's Equations*. New York, NY, USA: Oxford Univ. Press, 2003.

[6] W. C. Chew, E. Michielssen, J. Song, and J.-M. Jin, *Fast and Efficient Algorithms in Computational Electromagnetics*. Boston, MA, USA: Artech House, 2001.

[7] E. Michielssen and A. Boag, “A multilevel matrix decomposition algorithm for analyzing scattering from large structures,” *IEEE Trans. Antennas Propag.*, vol. 44, no. 8, pp. 1086–1093, Aug. 1996, doi: 10.1109/8.511816.

[8] E. Candès, L. Demanet, and L. Ying, “A fast butterfly algorithm for the computation of Fourier integral operators,” *Multiscale Model. Simulat.*, vol. 7, no. 4, pp. 1727–1750, Jan. 2009, doi: 10.1137/080734339.

[9] M. Bebendorf, C. Kuske, and R. Venn, “Wideband nested cross approximation for Helmholtz problems,” *Numerische Mathematik*, vol. 130, no. 1, pp. 1–34, May 2015, doi: 10.1007/s00211-014-0656-7.

[10] S. Börm, “Directional $H^2$-matrix compression for high-frequency problems,” *Numer. Linear Algebra Appl.*, vol. 24, no. 6, 2017, Art. no. e2112, doi: 10.1002/nla.2112.

[11] M. Bebendorf and R. Venn, “Constructing nested bases approximations from the entries of non-local operators,” *Numerische Mathematik*, vol. 121, no. 4, pp. 609–635, 2012.

[12] M. Li, M. A. Francavilla, F. Vipiana, G. Vecchi, and R. Chen, “Nested equivalent source approximation for the modeling of multiscale structures,” *IEEE Trans. Antennas Propag.*, vol. 62, no. 7, pp. 3664–3678, Jul. 2014.

[13] W. Hackbusch and Z. P. Nowak, “On the fast matrix multiplication in the boundary element method by panel clustering,” *Numerische Mathematik*, vol. 54, no. 4, pp. 463–491, 1989, doi: 10.1007/BF01396324.

[14] M. Bebendorf and S. Rjasanow, “Adaptive low-rank approximation of collocation matrices,” *Computing*, vol. 70, no. 1, pp. 1–24, 2003, doi: 10.1007/s00607-002-1469-6.

[15] 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, doi: 10.1109/TEMC.2005.857898.

[16] W. Chai and D. Jiao, “An $H^2$-matrix-based integral-equation solver of reduced complexity and controlled accuracy for solving electro-dynamic problems,” *IEEE Trans. Antennas Propag.*, vol. 57, no. 10, pp. 3147–3159, Oct. 2009, doi: 10.1109/TAP.2009.2028665.

[17] I. Fierro and C. Jerez-Hanckes, “Fast Calderón preconditioning for helmholtz boundary integral equations,” *J. Comput. Phys.*, vol. 409, no. 5, 2020, Art. no. 109355.

[18] F. Kpadonou, S. Chaillat, and P. Ciarlet, “On the efficiency of nested GMRES preconditioners for 3D acoustic and elastodynamic H-matrix accelerated boundary element methods,” *Comput. Math. Appl.*, vol. 80, no. 3, pp. 471–489, 2020.

[19] G. H. Golub and C. F. Van Loan, *Matrix Computations*. Baltimore, MD, USA: Johns Hopkins Univ. Press, 1996.
<figure><img src="image_22.png" /></figure>

[20] J. Mautz and R. Harrington, “An E-field solution for a conducting surface small or comparable to the wavelength,” *IEEE Trans. Antennas Propag.*, vol. 32, no. 4, pp. 330–339, Apr. 1984, doi: 10.1109/TAP.1984.1143316.

[21] J. R. Mautz and R. F. Harrington, “H-field, E-field, and combined-field solutions for conducting bodies of revolution,” *Archiv Elektronik Übertragungstechnik*, vol. 32, pp. 157–163, Apr. 1978.

[22] X. Antoine and M. Darbas, “An introduction to operator preconditioning for the fast iterative integral equation solution of time-harmonic scattering problems.” *Multiscale Sci. Eng.*, vol. 3, pp. 1–35, Feb. 2021, doi: 10.1007/s42493-021-00057-6.

[23] P. Yla-Oijala, S. P. Kiminki, J. Markkanen, and S. Jarvenpaa, “Error-controllable and well-conditioned mom solutions in computational electromagnetics: Ultimate surface integral-equation formulation [open problems in CEM],” *IEEE Antennas Propag. Mag.*, vol. 55, no. 6, pp. 310–331, Dec. 2013, doi: 10.1109/MAP.2013.6781749.

[24] B. Carpentieri, “Preconditioning for large-scale boundary integral equations in electromagnetics,” *IEEE Antennas Propag. Mag.*, vol. 56, no. 6, pp. 338–345, Dec. 2014, doi: 10.1109/MAP.2014.7011082.

[25] B. Carpentieri, “New trends in algebraic preconditioning,” in *New Trends in Computational Electromagnetics*, O. Ergul, Ed. London, U.K.: Inst. Eng. Technol., Jan. 2019, pp. 535–566, doi: 10.1049/SBEW533E ch13.

[26] K. Niino and N. Nishimura, “Calderón preconditioning approaches for PMCHWT formulations for Maxwell’s equations,” *Int. J. Numer. Model. Electron. Netw. Devices Fields*, vol. 25, nos. 5–6, pp. 558–572, 2012.

[27] D. Dobbelaere, D. De Zutter, J. Van Hese, J. Sercu, T. Boonen, and H. Rogier, “A Calderón multiplicative preconditioner for the electromagnetic Poincaré–Steklov operator of a heterogeneous domain with scattering applications,” *J. Comput. Phys.*, vol. 303, pp. 355–371, Dec. 2015.

[28] M. Gossye, M. Huynen, D. V. Ginste, D. De Zutter, and H. Rogier, “A Calderón preconditioner for high dielectric contrast media,” *IEEE Trans. Antennas Propag.*, vol. 66, no. 2, pp. 808–818, Feb. 2018.

[29] T. L. Chhim, A. Merlini, L. Rahmouni, J. E. O. Guzman, and F. P. Andriulli, “Eddy current modeling in multiply connected regions via a full-wave solver based on the quasi-Helmholtz projectors,” *IEEE Open J. Antennas Propag.*, vol. 1, pp. 534–548, 2020.

[30] Z. Peng, X.-c. Wang, and J.-F. Lee, “Integral equation based domain decomposition method for solving electromagnetic wave scattering from non-penetrable objects,” *IEEE Trans. Antennas Propag.*, vol. 59, no. 9, pp. 3328–3338, Sep. 2011.

[31] Z. Peng, R. Hiptmair, Y. Shao, and B. MacKie-Mason, “Domain decomposition preconditioning for surface integral equations in solving challenging electromagnetic scattering problems,” *IEEE Trans. Antennas Propag.*, vol. 64, no. 1, pp. 210–223, Jan. 2016.

[32] M. A. E. Bautista, F. Vipiana, M. A. Francavilla, J. A. T. Vasquez, and G. Vecchi, “A nonconformal domain decomposition scheme for the analysis of multiscale structures,” *IEEE Trans. Antennas Propag.*, vol. 63, no. 8, pp. 3548–3560, Aug. 2015.

[33] Z. Peng and J.-F. Lee, “Non-conformal domain decomposition method with second-order transmission conditions for time-harmonic electromagnetics,” *J. Comput. Phys.*, vol. 229, no. 16, pp. 5615–5629, 2010.

[34] Z. Peng, V. Rawat, and J.-F. Lee, “One way domain decomposition method with second order transmission conditions for solving electromagnetic wave problems,” *J. Comput. Phys.*, vol. 229, no. 4, pp. 1181–1197, 2010.

[35] Y.-J. Li and J.-M. Jin, “A new dual-primal domain decomposition approach for finite element simulation of 3-D large-scale electromagnetic problems,” *IEEE Trans. Antennas Propag.*, vol. 55, no. 10, pp. 2803–2810, Oct. 2007.

[36] W. K. Anderson, L. Wang, S. Kapadia, C. Tanis, and B. Hilbert, “Petrov–Galerkin and discontinuous-Galerkin methods for time-domain and frequency-domain electromagnetic simulations,” *J. Comput. Phys.*, vol. 230, no. 23, pp. 8360–8385, 2011.

[37] Z. Peng, K.-H. Lim, and J.-F. Lee, “A discontinuous Galerkin surface integral equation method for electromagnetic wave scattering from nonpenetrable targets,” *IEEE Trans. Antennas Propag.*, vol. 61, no. 7, pp. 3617–3628, Jul. 2013, doi: 10.1109/TAP.2013.2258394.

[38] S. Silver, *Microwave Antenna Theory and Design* (IEE Electromagnetic Waves Series), vol. 19. London, U.K.: P. Peregrinus Behalf Inst. Electr. Eng., 1984, p. 623.

[39] C. Müller, “Die Grundzüge einer mathematischen Theorie elektromagnetischer Schwingungen,” *Archiv Mathematik*, vol. 1, no. 4, pp. 296–302, Jul. 1948, doi: 10.1007/BF02038758.

[40] J. van Bladel, *Electromagnetic Fields*, 2nd ed. Hoboken, NJ, USA: Wiley, 2007.

[41] J.-C. Nedelec, *Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems* (Applied Mathematical Sciences), vol. 144. New York, NY, USA: Springer, 2001, p. 316.

[42] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” *IEEE Trans. Antennas Propag.*, vol. 30, no. 3, pp. 409–418, May 1982, doi: 10.1109/TAP.1982.1142818.

[43] K. Cools, F. P. Andriulli, D. De Zutter, and E. Michielssen, “Accurate and conforming mixed discretization of the MFIE,” *IEEE Antennas Wireless Propag. Lett.*, vol. 10, pp. 528–531, 2011, doi: 10.1109/LAWP.2011.2155022.

[44] Q. Chen and D. Wilton, “Electromagnetic scattering by three-dimensional arbitrary complex material/conducting bodies,” in *Proc. IEEE Antennas Propag. Soc. Int. Symp.*, 1990, pp. 590–593, doi: 10.1109/APS.1990.115179.

[45] A. Buffa and S. H. Christiansen, “A dual finite element complex on the barycentric refinement,” *Math. Comput.*, vol. 76, no. 260, pp. 1743–1770, Oct. 2007, doi: 10.1090/S0025-5718-07-01965-5.

[46] F. P. Andriulli et al., “A multiplicative Calderon preconditioner for the electric field integral equation,” *IEEE Trans. Antennas Propag.*, vol. 56, no. 8, pp. 2398–2412, Aug. 2008, doi: 10.1109/TAP.2008.926788.

[47] A. Gillman, P. M. Young, and P.-G. Martinsson, “A direct solver with O(N) complexity for integral equations on one-dimensional domains,” *Front. Math. China*, vol. 7, no. 2, pp. 217–247, Apr. 2012, doi: 10.1007/s11464-012-0188-3.

[48] H. Guo, Y. Liu, J. Hu, and E. Michielssen, “A butterfly-based direct integral equation solver using hierarchical LU factorization for analyzing scattering from electrically large conducting objects,” *IEEE Trans. Antennas and Propag.*, vol. 65, no. 9, pp. 4742–4750, 2017, doi: 10.1109/TAP.2017.2727511.

[49] Y. Brick, V. Lomakin, and A. Boag, “Fast direct solver for essentially convex scatterers using multilevel non-uniform grids,” *IEEE Trans. Antennas Propag.*, vol. 62, no. 8, pp. 4314–4324, Aug. 2014, doi: 10.1109/TAP.2014.2327651.

[50] P.-G. Martinsson, *Fast Direct Solvers for Elliptic PDEs* (CBMSNSF Regional Conference Series in Applied Mathematics), vol. 96. Philadelphia, PA, USA: Soc. Ind. Appl. Math., 2020, p. 315.

[51] J. R. Shewchuk, *An Introduction to the Conjugate Gradient Method Without the Agonizing Pain*, Carnegie Mellon Univ., Pittsburgh, PA, USA, 1994.

[52] M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving linear systems,” *J. Res. Nat. Bureau Stand.*, vol. 49, no. 1, pp. 409–436, 1952.

[53] O. Axelsson, *Iterative Solution Methods*. Cambridge, U.K.: Cambridge Univ. Press, 1996.

[54] D. R. Wilton and A. W. Glisson, “On improving the stability of the electric field integral equation at low frequencies,” in *Proc. USNC/URSI Spring Meeting Dig.*, 1981, p. 24.

[55] G. Vecchi, “Loop-star decomposition of basis functions in the discretization of the EFIE,” *IEEE Trans. Antennas Propag.*, vol. 47, no. 2, pp. 339–346, Feb. 1999, doi: 10.1109/8.761074.

[56] F. P. Andriulli and E. Michielssen, “Recovering the global loops by randomized projections,” in *Proc. IEEE Int. Symp. Antennas Propag.*, Jul. 2011, pp. 3197–3199, doi: 10.1109/APS.2011.5997213.

[57] D. Wilton, “Topological consideration in surface patch and volume cell modeling of electromagnetic scatterers,” in *Proc. URSI Int. Symp. Electromagn. Theory*, 1983, pp. 65–68.

[58] Z.-G. Qian and W. C. Chew, “Enhanced A-EFIE with perturbation method,” *IEEE Trans. Antennas Propag.*, vol. 58, no. 10, pp. 3256–3264, Oct. 2010, doi: 10.1109/TAP.2010.2055795.

[59] H. V. Henderson and S. R. Searle, “On deriving the inverse of a sum of matrices,” *SIAM Rev.*, vol. 23, no. 1, pp. 53–60, Jan. 1981. [Online]. Available: http://epubs.siam.org/doi/abs/10.1137/1023004
[60] F. P. Andriulli, I. Bogaert, K. Cools, and E. Michielssen, "A well-conditioned combined field integral equation based on quasi-Helmholtz projectors," in *Proc. IEEE ICEAA*, Sep. 2013, pp. 644–647, doi: 10.1109/ICEAA.2013.6632323.

[61] J.-S. Zhao and W. C. Chew, "Integral equation solution of Maxwell's equations from zero frequency to microwave frequencies," *IEEE Trans. Antennas Propag.*, vol. 48, no. 10, pp. 1635–1645, Oct. 2000, doi: 10.1109/8.899680.

[62] F. P. Andriulli, "Loop-star and loop-tree decompositions: Analysis and efficient algorithms," *IEEE Trans. Antennas Propag.*, vol. 60, no. 5, pp. 2347–2356, May 2012, doi: 10.1109/TAP.2012.2189723.

[63] T. F. Eibert, "Iterative-solver convergence for loop-star and loop-tree decompositions in method-of-moments solutions of the electric-field integral equation," *IEEE Antennas Propag. Mag.*, vol. 46, no. 3, pp. 80–85, Jun. 2004, doi: 10.1109/MAP.2004.1374101.

[64] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, "On a well-conditioned electric field integral operator for multiply connected geometries," *IEEE Trans. Antennas Propag.*, vol. 61, no. 4, pp. 2077–2087, Apr. 2013, doi: 10.1109/TAP.2012.2234072.

[65] L. N. Olson, and J. B. Schroder, "PyAMG: Algebraic multigrid solvers in python v4.0." 2018. [Online]. Available: https://github.com/pyamg/pyamg

[66] O. E. Livne and A. Brandt, "Lean algebraic multigrid (LAMG): Fast graph Laplacian linear solver," *SIAM J. Sci. Comput.*, vol. 34, no. 4, pp. 499–522, Jan. 2012, doi: 10.1137/110843563.

[67] W.-L. Wu, A. W. Glisson, and D. Kajfez, "A study of two numerical solution procedures for the electric field integral equation," *Appl. Comput. Electromagn. Soc. J.*, vol. 10, no. 3, pp. 69–80, 1995.

[68] M. Burton and S. Kashyap, "A study of a recent, moment-method algorithm that is accurate to very low frequencies," *Appl. Comput. Electromagn. Soc. J.*, vol. 10, no. 3, pp. 58–68, 1995.

[69] A. Bendali, "Numerical analysis of the exterior boundary value problem for the time-harmonic Maxwell equations by a boundary finite element method. I. The continuous problem," *Math. Comput.*, vol. 43, no. 167, pp. 29–46, Sep. 1984, doi: 10.1090/S0025-5718-1984-0744923-1.

[70] A. Bendali, "Numerical analysis of the exterior boundary value problem for the time-harmonic Maxwell equations by a boundary finite element method. II. The discrete problem," *Math. Comput.*, vol. 43, no. 167, pp. 47–68, Sep. 1984, doi: 10.1090/S0025-5718-1984-0744924-3.

[71] M. Taskinen and P. Ylä-Oijala, "Current and charge integral equation formulation," *IEEE Trans. Antennas Propag.*, vol. 54, no. 1, pp. 58–67, Jan. 2006, doi: 10.1109/TAP.2005.861580.

[72] Z. G. Qian and W. C. Chew, "An augmented electric field integral equation for high-speed interconnect analysis," *Microw. Opt. Technol. Lett.*, vol. 50, no. 10, pp. 2658–2662, Oct. 2008, doi: 10.1002/mop.23736.

[73] J. Cheng, R. J. Adams, J. C. Young, and M. A. Khayat, "Augmented EFIE with normally constrained magnetic field and static charge extraction," *IEEE Trans. Antennas Propag.*, vol. 63, no. 11, pp. 4952–4963, Nov. 2015, doi: 10.1109/TAP.2015.2478936.

[74] Y. Liu, X.-Y. Guo, and H.-J. Zhou, "A memory saving augmented EFIE with modified basis functions for low-frequency problems," *IEEE Trans. Antennas Propag.*, vol. 66, no. 3, pp. 1359–1365, Mar. 2018, doi: 10.1109/TAP.2018.2789983.

[75] J. Zhu, S. Omar, and D. Jiao, "Solution of the electric field integral equation when it breaks down," *IEEE Trans. Antennas Propag.*, vol. 62, no. 8, pp. 4122–4134, Aug. 2014, doi: 10.1109/TAP.2014.2322899.

[76] F. Vico, M. Ferrando, L. Greengard, and Z. Gimbutas, "The decoupled potential integral equation for time-harmonic electromagnetic scattering," *Commun. Pure Appl. Math.*, vol. 69, no. 4, pp. 771–812, 2016, doi: 10.1002/cpa.21585/abstract.

[77] Q. S. Liu, S. Sun, and W. C. Chew, "A potential-based integral equation method for low-frequency electromagnetic problems," *IEEE Trans. Antennas Propag.*, vol. 66, no. 3, pp. 1413–1426, Mar. 2018, doi: 10.1109/TAP.2018.2794388.

[78] W. C. Bandlow, G. Schneider, and A. F. Jacob, "Vector-valued wavelets with triangular support for method of moments applications," *IEEE Trans. Antennas Propag.*, vol. 53, no. 10, pp. 3340–3346, Oct. 2005, doi: 10.1109/TAP.2005.856309.

[79] F. P. Andriulli, A. Tobacco, and G. Vecchi, "A multiresolution approach to the electric field integral equation in antenna problems," *SIAM J. Sci. Comput.*, vol. 29, no. 1, pp. 1–21, Jan. 2007, doi: 10.1137/050634943.

[80] F. Vipiana, G. Vecchi, and P. Pirinoli, "A multiresolution system of Rao-Wilton-Glisson functions," *IEEE Trans. Antennas Propag.*, vol. 55, no. 3, pp. 924–930, Mar. 2007, doi: 10.1109/TAP.2007.891864.

[81] F. P. Andriulli, F. Vipiana, and G. Vecchi, "Hierarchical bases for nonhierarchical 3-D triangular meshes," *IEEE Trans. Antennas Propag.*, vol. 56, no. 8, pp. 2288–2297, Aug. 2008, doi: 10.1109/TAP.2008.926756.

[82] R.-S. Chen, J. Ding, D. Z.Ding, Z. H. Fan, and D. Wang, "A multiresolution curvilinear Rao-Wilton-Glisson basis function for fast analysis of electromagnetic scattering," *IEEE Trans. Antennas Propag.*, vol. 57, no. 10, pp. 3179–3188, Oct. 2009, doi: 10.1109/TAP.2009.2028599.

[83] F. P. Andriulli, A. Tobacco, and G. Vecchi, "Solving the EFIE at low frequencies with a conditioning that grows only logarithmically with the number of unknowns," *IEEE Trans. Antennas Propag.*, vol. 58, no. 5, pp. 1614–1624, May 2010, doi: 10.1109/TAP.2010.2044325.

[84] S. B. Adrian, F. P. Andriulli, and T. F. Eibert, "A hierarchical preconditioner for the electric field integral equation on unstructured meshes based on primal and dual Haar bases," *J. Comput. Phys.*, vol. 330, pp. 365–379, Feb. 2017, doi: 10.1016/j.jcp.2016.11.013.

[85] S. B. Adrian, F. P. Andriulli, and T. F. Eibert, "Hierarchical bases preconditioners for the electric field integral equation on multiply connected geometries," *IEEE Trans. Antennas Propag.*, vol. 62, no. 11, pp. 5856–5861, Nov. 2014, doi: 10.1109/TAP.2014.2347392.

[86] M. B. Stephanson and J.-F. Lee, "Preconditioned electric field integral equation using Calderón identities and dual loop/star basis functions," *IEEE Trans. Antennas Propag.*, vol. 57, no. 4, pp. 1274–1279, Apr. 2009, doi: 10.1109/TAP.2009.2016173.

[87] S. Yan, J.-M. Jin, and Z. Nie, "EFIE analysis of low-frequency problems with loop-star decomposition and Calderón multiplicative preconditioner," *IEEE Trans. Antennas Propag.*, vol. 58, no. 3, pp. 857–867, Mar. 2010, doi: 10.1109/TAP.2009.2039336.

[88] J.-S. Zhao, W. C. Chew, T. Cui, and Y. Zhang, "Cancellations of surface loop basis functions," in *Proc. Int. Symp. Antennas Propag. Soc.*, vol. 1, 2002, pp. 58–61.

[89] P. Werner, "On an integral equation in electromagnetic diffraction theory," *J. Math. Anal. Appl.*, vol. 14, no. 3, pp. 445–462, 1966.

[90] K. Cools, F. Andriulli, F. Olyslager, and E. Michielssen, "Nullspaces of MFIE and Calderón preconditioned EFIE operators applied to toroidal surfaces," *IEEE Trans. Antennas Propag.*, vol. 57, no. 10, pp. 3205–3215, Oct. 2009, doi: 10.1109/TAP.2009.2028669.

[91] I. Bogaert, K. Cools, F. P. Andriulli, and D. De Zutter, "Low frequency scaling of the mixed MFIE for scatterers with a non-simply connected surface," in *Proc. Int. Conf. Electromagn. Adv. Appl.*, Sep. 2011, pp. 951–954, doi: 10.1109/ICEAA.2011.6046467.

[92] A. Bossavit, "Magnetostatic problems in multiply connected regions: Some properties of the curl operator," *IEEE Proc. A Phys. Sci. Meas. Instrum. Manage. Educ. Rev.*, vol. 135, no. 3, pp. 179–187, Mar. 1988, doi: 10.1049/ip-a.1.19880028.

[93] M. Cessenat, *Mathematical Methods in Electromagnetism: Linear Theory and Applications* (Series on Advances in Mathematics for Applied Sciences), vol. 41. Singapore: World Sci., 1996, p. 376.

[94] Y. Zhang, T. J. Cui, W. C. Chew, and J.-S. Zhao, "Magnetic field integral equation at very low frequencies," *IEEE Trans. Antennas Propag.*, vol. 51, no. 8, pp. 1864–1871, Aug. 2003, doi: 10.1109/TAP.2003.814753.

[95] I. Bogaert, K. Cools, F. Andriulli, J. Peeters, and D. De Zutter, "Low frequency stability of the mixed discretization of the MFIE," in *Proc. 5th Eur. Conf. Antennas Propag. (EUCAP)*, Rome, Italy, Apr. 2011, pp. 2463–2465.

[96] A. Merlini, Y. Beghein, K. Cools, E. Michielssen, and F. P. Andriulli, "Magnetic and combined field integral equations based on the quasi-Helmholtz projectors," *IEEE Trans. Antennas Propag.*, vol. 68, no. 5, pp. 3834–3846, May 2020, doi: 10.1109/TAP.2020.2964941.

[97] S. B. Adrian, F. P. Andriulli, and T. F. Eibert, "On the hierarchical preconditioning of the combined field integral equation," *IEEE Antennas Wireless Propag. Lett.*, vol. 15, pp. 1897–1900, 2016, doi: 10.1109/LAWP.2016.2542828.
<figure><img src="image_23.png" /></figure>

[98] C. L. Epstein, Z. Gimbutas, L. Greengard, A. Klockner, and M. O'Neil, "A consistency condition for the vector potential in multiply-connected domains," *IEEE Trans. Magn.*, vol. 49, no. 3, pp. 1072-1076, Mar. 2013, doi: 10.1109/TMAG.2012.2223480.

[99] F. Vico, Z. Gimbutas, L. Greengard, and M. Ferrando-Bataller, "Overcoming low-frequency breakdown of the magnetic field integral equation," *IEEE Trans. Antennas Propag.*, vol. 61, no. 3, pp. 1285-1290, Mar. 2013, doi: 10.1109/TAP.2012.2230232.

[100] D. Wen, Z. Nie, L. Liu, and J. Tan, "Augmented combined field integral equation for low frequency scattering analysis," *IET Microw. Antennas Propag.*, vol. 11, no. 11, pp. 1510-1515, Sep. 2017, doi: 10.1049/iet-map.2016.0894.

[101] M. A. E. Bautista, M. A. Francavilla, F. Vipiana, and G. Vecchi, "A hierarchical fast solver for EFIE-MoM analysis of multiscale structures at very low frequencies," *IEEE Trans. Antennas Propag.*, vol. 62, no. 3, pp. 1523-1528, Mar. 2014, doi: 10.1109/TAP.2013.2296327.

[102] F. P. Andriulli and G. Vecchi, "A Helmholtz-stable fast solution of the electric field integral equation," *IEEE Trans. Antennas Propag.*, vol. 60, no. 5, pp. 2357-2366, May 2012, doi: 10.1109/TAP.2012.2189693.

[103] S. H. Christiansen and J.-C. Nédélec, "A preconditioner for the electric field integral equation based on Calderon formulas," *SIAM J. Numer. Anal.*, vol. 40, no. 3, pp. 1100-1135, Jan. 2002, doi: 10.1137/S0036142901388731.

[104] R. C. Kirby, "From functional analysis to iterative methods," *SIAM Rev.*, vol. 52, no. 2, pp. 269-293, Jan. 2010, doi: 10.1137/070706914.

[105] W. McLean and W. C. H. McLean, *Strongly Elliptic Systems and Boundary Integral Equations*. Cambridge, U.K.: Cambridge Univ. Press, 2000.

[106] K. E. Atkinson, *The Numerical Solution of Integral Equations of the Second Kind*. Cambridge, U.K.: Cambridge Univ. Press, 1999.

[107] R. J. Adams and N. J. Champagne, "A numerical implementation of a modified form of the electric field integral equation," *IEEE Trans. Antennas Propag.*, vol. 52, no. 9, pp. 2262-2266, Sep. 2004, doi: 10.1109/TAP.2004.834112.

[108] H. Contopanagos et al., "Well-conditioned boundary integral equations for three-dimensional electromagnetic scattering," *IEEE Trans. Antennas Propag.*, vol. 50, no. 12, pp. 1824-1830, Dec. 2002, doi: 10.1109/TAP.2002.803956.

[109] S. Borel, D. Levadoux, and F. Alouges, "A new well-conditioned integral formulation for Maxwell equations in three dimensions," *IEEE Trans. Antennas Propag.*, vol. 53, no. 9, pp. 2995-3004, Sep. 2005, doi: 10.1109/TAP.2005.854561.

[110] O. Steinbach and W. L. Wendland, "The construction of some efficient preconditioners in the boundary element method," *Adv. Comput. Math.*, vol. 9, nos. 1-2, pp. 191-216, 1998, doi: 10.1023/A:1018937506719.

[111] J. Peeters, K. Cools, I. Bogaert, F. Olyslager, and D. De Zutter, "Embedding Calderón multiplicative preconditioners in multilevel fast multipole algorithms," *IEEE Trans. Antennas Propag.*, vol. 58, no. 4, pp. 1236-1250, Apr. 2010.

[112] S. Sun, Y. G. Liu, W. C. Chew, and Z. Ma, "Calderón multiplicative preconditioned EFIE with perturbation method," *IEEE Trans. Antennas Propag.*, vol. 61, no. 1, pp. 247-255, Jan. 2013.

[113] Y. P. Chen, S. Sun, L. Jiang, and W. C. Chew, "A Calderón preconditioner for the electric field integral equation with layered medium green's function," *IEEE Trans. Antennas Propag.*, vol. 62, no. 4, pp. 2022-2030, Apr. 2014.

[114] A. C. Woo, H. T. G. Wang, M. J. Schuh, and M. L. Sanders, "EM programmer's notebook-benchmark radar targets for the validation of computational electromagnetics programs," *IEEE Antennas Propag. Mag.*, vol. 35, no. 1, pp. 84-89, Feb. 1993, doi: 10.1109/74.210840.

[115] P. Escapil-Inchauspé and C. Jerez-Hanckes, "Fast Calderón preconditioning for the electric field integral equation," *IEEE Trans. Antennas Propag.*, vol. 67, no. 4, pp. 2555–2564, Apr. 2019, doi: 10.1109/TAP.2019.2891608.

[116] P. Escapil-Inchauspé and C. Jerez-Hanckes, "Bi-parametric operator preconditioning," *Comput. Math. Applicat.*, vol. 102, pp. 220-232, 2021.

[117] S. B. Adrian, F. P. Andriulli, and T. F. Eibert, "On a refinement-free Calderón multiplicative preconditioner for the electric field integral equation," *J. Comput. Phys.*, vol. 376, pp. 1232-1252, Jan. 2019, doi: 10.1016/j.jcp.2018.10.009.

[118] R. Adams and G. Brown, "Stabilisation procedure for electric field integral equation," *Electron. Lett.*, vol. 35, no. 23, p. 2015, 1999, doi: 10.1049/el:19991339.

[119] S. H. Christiansen and J.-C. Nédélec, "Des préconditionneurs pour la résolution numérique des équations intégrales de frontière de l'électromagnétisme," *Comptes Rendus l'Académie des Sci. I Math.*, vol. 331, no. 9, pp. 733-738, Nov. 2000, doi: 10.1016/S0764-4442(00)01717-1.

[120] D. P. Levadoux, "Some preconditioners for the CFIE equation of electromagnetism," *Math. Methods Appl. Sci.*, vol. 31, no. 17, pp. 2015-2028, 2008, doi: 10.1002/mma.1004.

[121] O. Bruno, T. Elling, R. Paffenroth, and C. Turc, "Electromagnetic integral equations requiring small numbers of Krylov-subspace iterations," *J. Comput. Phys.*, vol. 228, no. 17, pp. 6169-6183, Sep. 2009, doi: 10.1016/j.jcp.2009.05.020.

[122] R. J. Adams, "Combined field integral equation formulations for electromagnetic scattering from convex geometries," *IEEE Trans. Antennas Propag.*, vol. 52, no. 5, pp. 1294-1303, May 2004, doi: 10.1109/TAP.2004.827246.

[123] M. Darbas, "Generalized combined field integral equations for the iterative solution of the three-dimensional Maxwell equations," *Appl. Math. Lett.*, vol. 19, no. 8, pp. 834-839, Aug. 2006, doi: 10.1016/j.aml.2005.11.005.

[124] F. Alouges, S. Borel, and D. P. Levadoux, "A stable well-conditioned integral equation for electromagnetism scattering," *J. Comput. Appl. Math.*, vol. 204, no. 2, pp. 440-451, Jul. 2007, doi: 10.1016/j.cam.2006.02.049.

[125] M. El Bouajaji, X. Antoine, and C. Geuzaine, "Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations," *J. Comput. Phys.*, vol. 279, pp. 241-260, Dec. 2014, doi: 10.1016/j.jcp.2014.09.011.

[126] S. Sun, Y. G. Liu, W. C. Chew, and Z. Ma, "Calderón multiplicative preconditioned EFIE with perturbation method," *IEEE Trans. Antennas Propag.*, vol. 61, no. 1, pp. 247-255, Jan. 2013, doi: 10.1109/TAP.2012.2220099.

[127] K. Niino, S. Akagi, and N. Nishimura, "A discretization method with the $H_{\text{div}}$ inner product for electric field integral equations," *IEEE Trans. Antennas Propag.*, vol. 65, no. 6, pp. 3102-3113, Jun. 2017, doi: 10.1109/TAP.2017.2696422.

[128] R. Hiptmair and C. Urzúa-Torres, "Compact equivalent inverse of the electric field integral operator on screens," *Integr. Equations. Oper. Theory*, vol. 92, no. 1, p. 9, Feb. 2020, doi: 10.1007/s00020-020-2564-9.

[129] W. Dahmen, S. Prößdorf, and R. Schneider, "Wavelet approximation methods for pseudodifferential equations I: Stability and convergence," to be published, doi: 10.20347/WIAS.PREPRINT.7.

[130] W. Dahmen, S. Prößdorf, and R. Schneider, "Wavelet approximation methods for pseudodifferential equations II: Matrix compression and fast solution," *Adv. Comput. Math.*, vol. 1, no. 3, pp. 259-335, Oct. 1993, doi: 10.1007/BF02072014.

[131] B. Alpert, G. Beylkin, R. Coifman, and V. Rokhlin, "Wavelet-like bases for the fast solution of second-kind integral equations," *SIAM J. Sci. Comput.*, vol. 14, no. 1, pp. 159-184, 1993, doi: 10.1137/0914010.

[132] W. Dahmen and R. Schneider, "Wavelets on manifolds I: Construction and domain decomposition," *SIAM J. Math. Anal.*, vol. 31, no. 1, pp. 184-230, Jan. 1999, doi: 10.1137/S0036141098333451.

[133] W. Dahmen and R. Stevenson, "Element-by-element construction of wavelets satisfying stability and moment conditions," *SIAM J. Numer. Anal.*, vol. 37, no. 1, pp. 319-352, Jan. 1999, doi: 10.1137/S0036142997330949.

[134] R. Stevenson, "Stable three-point wavelet bases on general meshes," *Numerische Mathematik*, vol. 80, no. 1, pp. 131-158, Jul. 1998, doi: 10.1007/s002110050363.

[135] J. Tausch and J. White, "Multiscale bases for the sparse representation of boundary integral operators on complex geometry," *SIAM J. Sci. Comput.*, vol. 24, no. 5, pp. 1610-1629, 2003, doi: 10.1137/S1064827500369451.

[136] H. Harbrecht and R. Schneider, "Biorthogonal wavelet bases for the boundary element method," *Mathematische Nachrichten*, vols. 269-270, no. 1, pp. 167-188, Jun. 2004, doi: 10.1002/mana.200310171.
[137] J. Tausch, "A variable order wavelet method for the sparse representation of layer potentials in the non-standard form," *J. Numer. Math.*, vol. 12, no. 3, pp. 233-254, 2004, doi: 10.1515/1569395041931473.

[138] H. Harbrecht, U. Kähler, and R. Schneider, "Wavelet Galerkin BEM on unstructured meshes," *Comput. Visualization Sci.*, vol. 8, nos. 3-4, pp. 189-199, 2005.

[139] W. Dahmen, H. Harbrecht, and R. Schneider, "Compression techniques for boundary integral equations—Asymptotically optimal complexity estimates," *SIAM J. Numer. Anal.*, vol. 43, no. 6, pp. 2251-2271, Jan. 2006, doi: 10.1137/S0036142903428852.

[140] H. Harbrecht and R. Schneider, "Wavelet Galerkin schemes for boundary integral equations—Implementation and quadrature," *SIAM J. Sci. Comput.*, vol. 27, no. 4, pp. 1347-1370, Jan. 2006, doi: 10.1137/S1064827503429387.

[141] F. Andriulli, G. Vecchi, F. Vipiana, P. Pirinoli, and A. Tabacco, "Optimal a priori clipping estimation for wavelet-based method of moments matrices," *IEEE Trans. Antennas Propag.*, vol. 53, no. 11, pp. 3726-3734, Nov. 2005, doi: 10.1109/TAP.2005.858855.

[142] F. Vipiana, P. Pirinoli, and G. Vecchi, "A multiresolution method of moments for triangular meshes," *IEEE Trans. Antennas Propag.*, vol. 53, no. 7, pp. 2247-2258, Jul. 2005, doi: 10.1109/TAP.2005.850710.

[143] F. Vipiana and G. Vecchi, "A novel, symmetrical solenoidal basis for the MoM analysis of closed surfaces," *IEEE Trans. Antennas Propag.*, vol. 57, no. 4, pp. 1294-1299, Apr. 2009, doi: 10.1109/TAP.2009.2015857.

[144] S. Sun, L. Jiang, and W. C. Chew, "Enhanced A-EFIE with Calderón multiplicative preconditioner," in *Proc. IEEE Antennas Propag. Soc. Int. Symp. (APSURSI)*, Jul. 2013, pp. 77-78, doi: 10.1109/APS.2013.6710699.

[145] G. Alléon, M. Benzi, and L. Giraud, "Sparse approximate inverse preconditioning for dense linear systems arising in computational electromagnetics," *Numer. Algorithms*, vol. 16, no. 1, pp. 1-15, Feb. 1997, doi: 10.1023/A:1019170609950.

[146] B. Carpentieri, I. S. Duff, and L. Giraud, "Sparse pattern selection strategies for robust Frobenius-norm minimization preconditioners in electromagnetism," *Numer. Linear Algebra Appl.*, vol. 7, nos. 7-8, pp. 667-685, 2000, doi: 10.1002/1099-1506(200010/12)7:7<8:667>:AID-NLA218>3.0.CO;2-X.

[147] B. Carpentieri, I. S. Duff, L. Giraud, and M. M. M. Made, "Sparse symmetric preconditioners for dense linear systems in electromagnetism," *Numer. Linear Algebra Appl.*, vol. 11, nos. 8-9, pp. 753-771, 2004, doi: 10.1002/nla.345.

[148] J. Lee, J. Zhang, and C.-C. Lu, "Sparse inverse preconditioning of multilevel fast multipole algorithm for hybrid integral equations in electromagnetics," *IEEE Trans. Antennas Propag.*, vol. 52, no. 9, pp. 2277-2287, Sep. 2004, doi: 10.1109/TAP.2004.834084.

[149] B. Carpentieri, I. S. Duff, L. Giraud, and G. Sylvand, "Combining fast multipole techniques and an approximate inverse preconditioner for large electromagnetism calculations," *SIAM J. Sci. Comput.*, vol. 27, no. 3, pp. 774-792, Jan. 2005, doi: 10.1137/040603917.

[150] T. Malas and L. Gürel, "Accelerating the multilevel fast multipole algorithm with the sparse-approximate-inverse (SAI) preconditioning," *SIAM J. Sci. Comput.*, vol. 31, no. 3, pp. 1968-1984, Jan. 2009, doi: 10.1137/070711098.

[151] C. Delgado and M. F. Cátedra, "Sparse approximate inverse preconditioner with parametric sparsity pattern applied to the macrobasis function methods," *IEEE Antennas Wireless Propag. Lett.*, vol. 17, no. 5, pp. 849-852, May 2018, doi: 10.1109/LAWP.2018.2819085.

[152] X.-M. Pan and X.-Q. Sheng, "Sparse approximate inverse preconditioner for multiscale dynamic electromagnetic problems," *Radio Sci.*, vol. 49, no. 11, pp. 1041-1051, Nov. 2014, doi: 10.1002/2014RS005387.

[153] C. Delgado, E. García, I. Somolinos, and M. F. Cátedra, "Hybrid parallelisation scheme for the application of distributed near-field sparse approximate inverse preconditioners on high-performance computing clusters," *IET Microw. Antennas Propag.*, vol. 14, no. 4, pp. 320-328, Jan. 2020, doi: 10.1049/iet-map.2019.0789.

[154] B. Carpentieri, "A matrix-free two-grid preconditioner for solving boundary integral equations in electromagnetism," *Computing*, vol. 77, no. 3, pp. 275-296, May 2006, doi: 10.1007/s00607-006-0161-7.

[155] K. Sertel and J. L. Volakis, "Incomplete LU preconditioner for FMM implementation," *Microw. Opt. Technol. Lett.*, vol. 26, no. 4, pp. 265-267, Aug. 2000, doi: 10.1002/1098-2760(20000820)26:4<265::AID-MOP18>3.0.CO;2-O.

[156] J. Lee, J. Zhang, and C.-C. Lu, "Incomplete LU preconditioning for large scale dense complex linear systems from electromagnetic wave scattering problems," *J. Comput. Phys.*, vol. 185, no. 1, pp. 158-175, Feb. 2003, doi: 10.1016/S0021-9991(02)00052-9.

[157] T. Malas and L. Gürel, "Incomplete LU preconditioning with the multilevel fast multipole algorithm for electromagnetic scattering," *SIAM J. Sci. Comput.*, vol. 29, no. 4, pp. 1476-1494, Jan. 2007, doi: 10.1137/060659107.

[158] B. Carpentieri and M. Bollhöfer, "Symmetric inverse-based multilevel ILU preconditioning for solving dense complex non-hermitian systems in electromagnetics," *Progr. Electromagn. Res.*, vol. 128, pp. 55-74, May 2012, doi: 10.2528/PIER12041006.

[159] T. Eibert, "Iterative near-zone preconditioning of iterative method of moments electric field integral equation solutions," *IEEE Antennas Wireless Propag. Lett.*, vol. 2, pp. 101-102, 2003, doi: 10.1109/LAWP.2003.814777.

[160] J. Mautz and R. Harrington, "A combined-source solution for radiation and scattering from a perfectly conducting body," *IEEE Trans. Antennas Propag.*, vol. 27, no. 4, pp. 445-454, Jul. 1979, doi: 10.1109/TAP.1979.1142115.

[161] H. Bagci, F. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, "A Calderón multiplicative preconditioner for the combined field integral equation," *IEEE Trans. Antennas Propag.*, vol. 57, no. 10, pp. 3387-3392, Oct. 2009.

[162] Y. Boubendir and C. Turc, "Well-conditioned boundary integral equation formulations for the solution of high-frequency electromagnetic scattering problems," *Comput. Math. Appl.*, vol. 67, no. 10, pp. 1772-1805, 2014, doi: 10.1016/j.camwa.2014.04.003.

[163] Y. Boubendir and C. Turc, "Wave-number estimates for regularized combined field boundary integral operators in acoustic scattering problems with Neumann boundary conditions," *IMA J. Numer. Anal.*, vol. 33, no. 4, pp. 1176-1225, Oct. 2013, doi: 10.1093/imanum/drs038.

[164] S. N. Chandler-Wilde, I. G. Graham, S. Langdon, and M. Lindner, "Condition number estimates for combined potential boundary integral operators in acoustic scattering," *J. Integr. Equ. Appl.*, vol. 21, no. 2, pp. 229-279, 2009.

[165] M. Darbas, "Préconditionneurs analytiques de type Calderon pour les formulations intégrales des problèmes de diffraction d'ondes," Ph.D. dissertation, Département de Génie Mathématique, Nat. Inst. Appl. Sci. Toulouse, Toulouse, France, 2004.

[166] F. P. Andriulli, I. Bogaert, and K. Cools, "On the high frequency behavior and stabilization of a preconditioned and resonance-free formulation," in *Proc. Int. Conf. Electromagn. Adv. Appl. (ICEAA)*, Sep. 2015, pp. 1321-1324, doi: 10.1109/ICEAA.2015.7297331.

[167] A. Dely, A. Merlini, S. B. Adrian, and F. P. Andriulli, "On preconditioning electromagnetic integral equations in the high frequency regime via Helmholtz operators and quasi-Helmholtz projectors," in *Proc. Int. Conf. Electromagn. Adv. Appl. (ICEAA)*, Granada, Spain, Sep. 2019, pp. 1338-1341, doi: 10.1109/ICEAA.2019.8879366.

[168] J. R. Shewchuk, *What Is a Good Linear Finite Element? Interpolation, Conditioning, Anisotropy, and Quality Measures*, Dept. Electr. Eng. Comput. Sci., Univ. California, Berkeley, CA, USA, 2002.

[169] M. B. Stephanson and J.-F. Lee, "Automatic basis function regularization for integral equations in the presence of ill-shaped mesh elements," *IEEE Trans. Antennas Propag.*, vol. 61, no. 8, pp. 4139-4147, Aug. 2013, doi: 10.1109/TAP.2013.2260118.

[170] J. Li, D. Dault, B. Liu, Y. Tong, and B. Shanker, "Subdivision based isogeometric analysis technique for electric field integral equations for simply connected structures," *J. Comput. Phys.*, vol. 319, pp. 145-162, Aug. 2016, doi: 10.1016/j.jcp.2016.04.008.

[171] J. Dölz, S. Kurz, S. Schöps, and F. Wolf, "Isogeometric boundary elements in electromagnetism: Rigorous analysis, fast methods, and examples," *SIAM J. Sci. Comput.*, vol. 41, no. 5, pp. B983-B1010, Jan. 2019, doi: 10.1137/18M1227251.

[172] J. Dölz, H. Harbrecht, S. Kurz, M. Multerer, S. Schöps, and F. Wolf, "Bembel: The fast isogeometric boundary element C++ library for Laplace, Helmholtz, and electric wave equation," *SoftwareX*, vol. 11, Jan.-Jun. 2020, Art. no. 100476, doi: 10.1016/j.softx.2020.100476.
<figure><img src="image_24.png" /></figure>

[173] J. Dolz, S. Kurz, S. Schops, and F. Wolf, "A numerical comparison of an isogeometric and a parametric higher order Raviart-Thomas approach to the electric field integral equation," *IEEE Trans. Antennas Propag.*, vol. 68, no. 1, pp. 593-597, Jan. 2020, doi: 10.1109/TAP.2019.2935778.

[174] J. C. Nédélec, "Mixed finite elements in $R^3$," *Numerische Mathematik*, vol. 35, no. 3, pp. 315-341, Sep. 1980, doi: 10.1007/BF01396415.

[175] J.-C. Nédélec, "A new family of mixed finite elements in $R^3$," *Numerische Mathematik*, vol. 50, no. 1, pp. 57-81, 1986.

[176] R. Graglia, D. Wilton, and A. Peterson, "Higher order interpolatory vector bases for computational electromagnetics," *IEEE Trans. Antennas Propag.*, vol. 45, no. 3, pp. 329-342, Mar. 1997, doi: 10.1109/8.558649.

[177] R. D. Graglia, A. F. Peterson, and F. P. Andriulli, "Curl-conforming hierarchical vector bases for triangles and tetrahedra," *IEEE Trans. Antennas Propag.*, vol. 59, no. 3, pp. 950-959, Mar. 2011, doi: 10.1109/TAP.2010.2103012.

[178] R. D. Graglia and A. F. Peterson, "Hierarchical curl-conforming Nédélec elements for quadrilateral and brick cells," *IEEE Trans. Antennas Propag.*, vol. 59, no. 8, pp. 2766-2773, Aug. 2011, doi: 10.1109/TAP.2011.2158789.

[179] B. M. Notaros, "Higher order frequency-domain computational electromagnetics," *IEEE Trans. Antennas Propag.*, vol. 56, no. 8, pp. 2251-2276, Aug. 2008, doi: 10.1109/TAP.2008.926784.

[180] R. D. Graglia and A. F. Peterson, *Higher-Order Techniques in Computational Electromagnetics* (Mario Boella Series on Electromagnetism in Information & Communication). Edison, NJ, USA: SciTech Publ., 2015, p. 392.

[181] J. Wang and J. P. Webb, "Hierarchical vector boundary elements and p-adaption for 3-D electromagnetic scattering," *IEEE Trans. Antennas Propag.*, vol. 45, no. 12, pp. 1869-1879, Dec. 1997, doi: 10.1109/8.650207.

[182] J. P. Webb, "Hierarchical vector basis functions of arbitrary order for triangular and tetrahedral finite elements," *IEEE Trans. Antennas Propag.*, vol. 47, no. 8, pp. 1244-1253, Aug. 1999, doi: 10.1109/8.791939.

[183] D.-K. Sun, J.-F. Lee, and Z. Cendes, "Construction of nearly orthogonal Nedelec bases for rapid convergence with multilevel preconditioned solvers," *SIAM J. Sci. Comput.*, vol. 23, no. 4, pp. 1053-1076, Jan. 2001, doi: 10.1137/S1064827500367531.

[184] P. Ingelstrom, "A new set of H(curl)-conforming hierarchical basis functions for tetrahedral meshes," *IEEE Trans. Microw. Theory Techn.*, vol. 54, no. 1, pp. 106-114, Jan. 2006, doi: 10.1109/TMTT.2005.860295.

[185] E. Jorgensen, J. Volakis, P. Meincke, and O. Breinbjerg, "Higher order hierarchical Legendre basis functions for electromagnetic modeling," *IEEE Trans. Antennas Propag.*, vol. 52, no. 11, pp. 2985-2995, Nov. 2004, doi: 10.1109/TAP.2004.835279.

[186] M. M. Kostić, and B. M. Kolundžija, "Maximally orthogonalized higher order bases over generalized wires, quadrilaterals, and hexahedra," *IEEE Trans. Antennas Propag.*, vol. 61, no. 6, pp. 3135-3148, Jun. 2013, doi: 10.1109/TAP.2013.2249036.

[187] I. Bogaert and F. Andriulli, "Maximally orthogonal high-order basis functions have a well-conditioned Gram matrix," *IEEE Trans. Antennas Propag.*, vol. 62, no. 8, pp. 4096-4104, Aug. 2014, doi: 10.1109/TAP.2014.2323081.

[188] S. V. Savić, M. M. Iilć, and B. M. Kolundžija, "Maximally orthogonalized higher order basis functions in large-domain finite element modeling in electromagnetics," *IEEE Trans. Antennas Propag.*, vol. 68, no. 8, pp. 6455-6460, Aug. 2020, doi: 10.1109/TAP.2020.2970038.

[189] F. Valdés, F. P. Andriulli, K. Cools, and E. Michielssen, "High-order div- and quasi curl-conforming basis functions for Calderón multiplicative preconditioning of the EFIE," *IEEE Trans. Antennas Propag.*, vol. 59, no. 4, pp. 1321-1337, Apr. 2011, doi: 10.1109/TAP.2011.2109692.

[190] R. Hiptmair, "Operator preconditioning," *Comput. Math. Appl.*, vol. 52, no. 5, pp. 699-706, Sep. 2006, doi: 10.1016/j.camwa.2006.10.008.

[191] G. C. Hsiao and R. E. Kleinman, "Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics," *IEEE Trans. Antennas Propag.*, vol. 45, no. 3, pp. 316-328, Mar. 1997, doi: 10.1109/8.558648.

[192] F. Vico, L. Greengard, and Z. Gimbutas, "Boundary integral equation analysis on the sphere," *Numerische Mathematik*, vol. 128, no. 3, pp. 463-487, Nov. 2014, doi: 10.1007/s00211-014-0619-z.

[193] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, *NIST Handbook of Mathematical Functions*. Cambridge, U.K.: Cambridge Univ. Press, 2010, p. 951.

[194] J. D. Jackson, *Classical Electrodynamics*, 3rd ed. New York, NY, USA: Wiley, 1999, p. 808.

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

**SIMON B. ADRIAN** (Member, IEEE) received the Bachelor of Science degree in electrical engineering and information technology from the Technical University of Munich (TUM), Munich, Germany, in 2009, the Master of Science degree in electrical and computer engineering from the Georgia Institute of Technology, Atlanta, GA, USA, in 2010, the Diplom-Ingenieur degree in electrical engineering and information technology from the TUM in 2012, and the Doktor-Ingenieur degree from the TUM and the École Nationale supérieure

Mines-Télécom Atlantique Bretagne Pays de la Loire (IMT Atlantique), Brest, France, in 2018.

From 2012 to 2019, he was Research Assistant with TUM. In 2018, he was a Visiting Professor with the Politecnico di Torino, Turin, Italy. From 2019 to 2020, he was a Senior Engineer with Infineon AG, Neubiberg, Germany. Since 2020, he is an Assistant Professor with the University of Rostock. Areas of application include antenna modeling and bio-electromagnetic problems. His research interest is in computational electromagnetics with a focus on integral equation solvers. He is interested in preconditioning techniques, low-frequency stable formulations, and fast solvers.

Dr. Adrian is laureate of an IEEE Antennas and Propagation Society Pre-Doctoral Research Award and an IEEE Antennas and Propagation Society Doctoral Research Award. He was awarded with the Kurt-Fischer-Preis 2012 for his diploma thesis. His contribution to the student paper competition at IEEE APS/URSI Symposium in 2013 was selected as an honorable mention paper. He received the Second Prize in the Third International Union Radio-Scientifique Internationale (URSI) Student Prize Paper Competition of the URSI General Assembly and Scientific Symposium 2014, the Best Young-Scientist-Paper Award of the Kleinheubacher Tagung 2014, the Third Prize Student Paper Award of the IEEE APS/URSI Symposium 2015, the First Prize EMTS 2016 Young Scientist Best Paper Award of the URSI International Symposium on Electromagnetic Theory 2016, the Third Prize Student Paper Award of the IEEE APS/URSI Symposium 2018, and the Dr.-Georg-Spinner-Hochfrequenz-Preis for the doctoral thesis in 2018. He has been an Associate Editor of the IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION since 2020 and the IEEE Antennas and Propagation Magazine since 2021. He has been a member of the IEEE Antennas and Propagation Society Education Committee since 2016, where he has been chair of the Student Activities Subcommittee since 2021. He has been member of URSI since 2018.

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

**ALEXANDRE DÉLY** received the M.Sc.Eng. degree from the École Nationale Supérieure des Télécommunications de Bretagne, France, in 2015, and the Ph.D. degree from the École Nationale Supérieure Mines-Télécom Atlantique, France, and from the University of Nottingham, U.K., in 2019.

He is currently working with the Politecnico di Torino, Turin, Italy. His research focuses on preconditioned and fast solution of boundary element methods, frequency domain, and time domain integral equations.
<figure><img src="image_27.png" /></figure>

**DAVIDE CONSOLI** (Graduate Student Member, IEEE) received the M.Sc. (Laurea) degree from the Politecnico di Torino, Italy, in 2018, where he is currently pursuing the Ph.D. degree. His research interests are in computational electromagnetics with a focus on fast and preconditioned integral equation formulations.

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

**FRANCESCO P. ANDRIULLI** (Senior Member, IEEE) received the Laurea degree in electrical engineering from the Politecnico di Torino, Turin, Italy, in 2004, the M.Sc. degree in electrical engineering and computer science from the University of Illinois at Chicago in 2004, and the Ph.D. degree in electrical engineering from the University of Michigan at Ann Arbor in 2008.

From 2008 to 2010, he was a Research Associate with the Politecnico di Torino. From 2010 to 2017, he was an Associate Professor (2010–2014) and then a Full Professor with the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique, formerly ENST Bretagne), Brest, France. Since 2017, he has been a Full Professor with the Politecnico di Torino. His research interests are in computational electromagnetics with focus on frequency- and time-domain integral equation solvers, well-conditioned formulations, fast solvers, low-frequency electromagnetic analyses, and modeling techniques for antennas, wireless components, microwave circuits, and biomedical applications with a special focus on Brain Imaging.

Prof. Andriulli was the recipient of the Best Student Paper Award at the 2007 URSI North American Radio Science Meeting. He received the First Place Prize of the student paper context of the 2008 IEEE Antennas and Propagation Society International Symposium. He was the recipient of the 2009 RMTG Award for junior researchers and was awarded two International Union of Radio Science (URSI) Young Scientist Awards at the International Symposium on Electromagnetic Theory in 2010 and 2013, where he was also awarded the Second Prize in the Best Paper Contest. He also received the 2015 ICEAA IEEE-APWC Best Paper Award. In addition, he coauthored with his students and collaborators other three first prize conference papers (EMTS 2016, URSI-DE Meeting 2014, ICEAA 2009), a second prize conference paper (URSI GASS 2014), two third prize conference papers (IEEE-APS 2015, IEEE-APS 2018), three honorable mention conference papers (ICEAA 2011, URSI/IEEE-APS 2013, URSI/IEEE-APS 20) and other three finalist conference papers (URSI/IEEE-APS 2012, URSI/IEEE-APS 2007, URSI/IEEE-APS 2006). He received the 2014 IEEE AP-S Donald G. Dudley Jr. Undergraduate Teaching Award, the triennium 2014–2016 URSI Issac Koga Gold Medal, and the 2015 L. B. Felsen Award for Excellence in Electrodynamics. He is the Editor in Chief of the *IEEE Antennas and Propagation Magazine*. He serves as a Track Editor for the IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, and as an Associate Editor for IEEE ACCESS, *URSI Radio Science Letters*, and *IET Microwaves, Antennas & Propagation*. He is the PI of the ERC Consolidator Grant “321.” He is a member of Eta Kappa Nu, Tau Beta Pi, Phi Kappa Phi, and URSI.

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

**ADRIEN MERLINI** (Member, IEEE) received the M.Sc.Eng. degree from the École Nationale Supérieure des Télécommunications de Bretagne, France, in 2015, and the Ph.D. degree from the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique), France, in 2019.

From 2018 to 2019, he was a visiting Ph.D. student with the Politecnico di Torino, Italy, which he then joined as a Research Associate. Since 2019, he has been an Associate Professor with the Microwave Department, IMT Atlantique. His

research interests include preconditioning and acceleration of integral equation solvers for electromagnetic simulations and their application in brain imaging.

Dr. Merlini received the Young Scientist Award at the URSI GASS 2020 meeting. His coauthored conference paper recipient of an honorable mention at the URSI/IEEE-APS 2020. He is a member of IEEE-HKN, the IEEE Antennas and Propagation Society, URSI France, and of the Lab-STICC Laboratory.