Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning

Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning / Adrian, Simon B.; Alexandre, Dély; Consoli, Davide; Merlini, Adrien; Andriulli, Francesco P.. In: IEEE OPEN JOURNAL OF ANTENNAS AND PROPAGATION. ISSN 2637-6431. ELETTRONICO. (In corso di stampa). [10.1109/OJAP.2021.3121097] Original Electromagnetic Integral Equations: Insights in Conditioning and Preconditioning

products in O ( log ) or even ( ) complexity, where 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 ( iter log ) (or ( iter ) in the low-frequency regime), where iter is the number of iterations.
The number of iterations 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 and, when the BEM formulation is set in the frequency domain, of the wavenumber . This can potentially result in a solution complexity greater, and sometimes much greater, than ( 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) highfrequency issues (including internal resonances and the highfrequency 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 highfrequency issues. A few years ago, Ylä-Oĳala 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 not that a substantial amount of literature and quite effective preconditioned methods have been presented for modelling 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 literature [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 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 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 low of mesh and basis functions quality on the overall conditioning and Section VIII presents the conclusions and final considerations.

II. N B
We are interested in solving the electromagnetic scattering problem where a time-harmonic, electromagnetic wave ( i , i ) in a space with permittivity and permeability impinges on a connected domain − ⊂ R 3 with PEC boundary ≔ − resulting in the scattered wave ( s , s ). The total electric ≔ i + s and magnetic ≔ i + s fields satisfy Maxwell's equations where + ≔ − c , ≔ √ is the wave number, the angular frequency, and , must satisfy the boundary conditions for PEC boundarieŝ where is the induced electric surface current density. In addition, s and s must satisfy the Silver-Müller radiation condition [38], [39] lim →∞ s × − s = 0 .
We assumed (and suppressed) a time dependency of e −i and normalized with the wave impedance ≔ / . To find ( s , s ), we can solve the EFIE for , whereˆ is the surface normal vector directed into + and is the EFIE operator composed of the vector potential operator and the scalar potential operator where is the free-space Green's function. A definition of the surface differential operators grad and div can be found in [40,Appendix 3] or [41,Chapter 2]. Once is obtained, s , s can be computed using the free-space radiation operators.
Alternatively, one can solve the MFIE for the exterior scattering problem where I is identity operator, M + is the MFIE operator for the exterior scattering problem, and The MFIE operator for the interior scattering problem is M − ≔ −I/2+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] − T + (1 − )ˆ × M + = ˆ × i + (1 − )ˆ ×ˆ × i (13) which is uniquely solvable for all frequencies.
For the discretization of the EFIE, we employ Rao-Wilton-Glisson (RWG) basis functions ∈ which are here-in contrast to their original definition in [42]-not normalized with the edge length, that is, using the convention depicted in Figure 1. Following a Petrov-Galerkin approach, we obtain the system of equations that can be solved to obtain an approximation of the solution in the form and where Even though we are testing withˆ × , the resulting system matrix T 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ˆ × 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 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 ∈ , 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 where and where the Gram matrix for any two function spaces and is defined as with ∈ and ∈ . For the discretization of the CFIE, we have with the combination parameter 0 < < 1.
III. C N , I S , C C 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 an iterative solver should stop when where > 0 is the solver tolerance and x ( ) the approximate solution after the 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 for obtaining an approximation of x within the tolerance , a relation between iter and is needed. One way to obtain such a relationship is via the condition number of the matrix, which is defined as where 2 is the spectral norm, and 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 ( ) ≔ x ( ) − x given by [51] where A is the energy norm defined by x A ≔ x † Ax 1/2 and x † denotes the conjugate transpose of x. If the objective is to reduce the relative error e ( ) / e (0) below and by considering limits for cond A ≫ 1, one notes [51] that iterations are at most needed (assuming an exact arithmetic). If the condition number grows linearly in , as observed for the EFIE when the mesh is uniformly refined, this implies that the complexity is at most O ( 1.5 log ).
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 The price for this, however, is that the condition number of the resulting system matrix is (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 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 and or imaginary axis and are either strictly positive or negative (if they are on the imaginary axis, then positive or negative with respect to Im ( )). 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 , which is independent from the average edge length ℎ of the mesh. For dynamic problems, however, we also need to study the condition number as a function of the frequency ≔ /(2π), and one must specify if a formulation is well-conditioned with respect to ℎ, to , to both, or only in a particular regime, for example, for frequencies where the corresponding wavelength is larger then the diameter of .
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 where x = P R y and the matrices should be chosen such that, if possible, where is a constant both independent of ℎ and (in which case the preconditioner is optimal). Normally, the matrixmatrix 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 CFIE operator, discuss the causes of their ill-conditioning as well as potential remedies.
IV. L -F S The low-frequency breakdown of the EFIE, that is, the growth of the condition number when the frequency 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 can be decomposed into a solenoidal, irrotational, and a harmonic vector field, which in the case of a tangential surface vector field such as leads to [41, p. (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. a quasi-Helmholtz decomposition is possible, where is decomposed into a solenoidal, a non-solenoidal, and a quasiharmonic current density. It is not possible to obtain irrotational or harmonic current densities, since the curl of divconforming (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 , a quasi-Helmholtz decomposition decomposes the discrete solution j as  [54] and where j , j , and j are the vectors containing the associated expansion coefficients; moreover, V is the number of vertices and 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 , , and can be represented in terms of RWG functions [54]. Thus the expansion coefficients are linked by linear transformation matrices Λ, H, and Σ. For the loop transformation matrix, we have where is the th vertex of the mesh (inner vertex if is open), and for the star transformation matrix where is the 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 The linear combinations of RWGs implied by the coefficient vectors j sol , j nsol , and j 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 sol , j nsol , and j 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 Λ and Σ a meaning that goes beyond merely interpreting them as basis transformation matrices. The matrices Λ and Σ are edge-node and edge-cell incidence matrices of the graph defined by the mesh and they are orthogonal, that is, Σ T Λ = 0. It follows that j sol and j nsol are 2 -orthogonal, that is, j T nsol j sol = 0. We find this noteworthy for two reasons: i) the loop and star functions are, in general, not 2 -orthogonal (after all, is not irrotational); ii) the 2 -orthogonality is not true for other quasi-Helmholtz decompositions such as the loop-tree basis. In light of this consideration, the matrices Λ and Σ could be interpreted as the graph curl (Λ) and graph gradient (Σ) of the standard mesh, an interpretation that further increases the correspondence with the continuous decomposition (32).
For global loops , 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 Σ T H = 0; however, Λ T H = 0 is, in general, not true. This property can be enforced by constructing H as the right nullspace of Λ Σ T . Such a construction is possible, for example, via a full singular value decomposition (SVD), or, via more computationally efficient randomized projections [56]. However, the computational cost is higher, in general, compared with using a global loopfinding algorithm, in particular, since H will be a dense matrix.
A similar decomposition can be obtained for dual functions Note that the same matrices Σ and Λ are present both in the decomposition of RWG functions and in the one of dual functions. However, while for RWGs the transformation matrix Λ describes solenoidal functions and the transformation matrix Σ describes non-solenoidal functions, the opposite is true for the dual functions: it is Σ that describes solenoidal functions, 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ˆ × are a rotation by 90°, which is not div-confirming); given that curl ≔ ∇ ×ˆ , 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 H = H. This equality is only true if H is the nullspace of Λ Σ T , a condition, which evidently leads to the aforementioned unique definition of H. Even though the construction of H as the nullspace of Λ Σ 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 analytic 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 nonsolenoidal 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 EFIEglobal 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 T LS ≔ Λ Σ T T Λ Σ can be represented in block matrix form as and the corresponding matrix equation now reads T LS LS = Λ Σ T i , where = Λ Σ LS . In these definitions, the Λ and Σ matrices refer to the full-rank transformation matrices in which linearly dependent columns have been removed: for each connected component of one star basis function (column of Σ) must always be removed and one loop basis function must be removed (column of Λ) 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 Λ T T Φ, = 0 and T Φ, Λ = 0, which follow directly from the divergence-free nature of solenoidal functions, 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 T A, and T Φ, for → 0. For instance, We can deduce that, in general, This process can be repeated for both T A, and T Φ, when both expansion and testing functions are non-solenoidal and when at least one of the two is solenoidal. In summary, By symmetry, both ˆ × , T and ˆ × , T have the same low-frequency behavior as ˆ × , T . The scaling of the behavior of the block matrix is now straightforward to obtain and the dominant behavior of T 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 = diag −1/2 1/2 in which the block dimensions are consistent with that of the loop star decomposition matrix. Clearly, is a well-conditioned matrix, in the sense that lim →0 cond D T LS D ≕ is finite. It then follows that and thus lim →0 cond T LS = O ( −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 T LS , which proves that Considering these results and that the loop-star transformation matrix is invertible and frequency independent, we conclude that cond (T ) ∼ −2 when → 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 where i PW is the electric field of the incident plane-wave. It is crucial to remember that when the standard EFIEwith 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 ( 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 planewave 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 which, in combination with the right-hand side results yields the behavior of the solution coefficients 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 nonsolenoidal 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 where j = Λ Σ D 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 = Λ 0 D j DLS and j sol = 0 Σ D 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 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 in 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 Λ Σ is replaced by an RWG-to-loop-tree mapping Λ Θ and T LS becomes is the general term of the RWG-to-tree transformation matrix. The resulting preconditioned equation is where j = Λ Θ D j DLT . At first glance, the computational overhead of the two methods seems low, since Λ, Σ, Θ, and D are sparse matrices. However, while both methods adequately address the lowfrequency 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 faster than the loop-star preconditioned [63], which can be explained by the fact that Λ and Σ 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 conditioning of the of the loop-star and looptree 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 loopstar/tree decompositions are helpful in analyzing the reasons behind of 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) the quasi-Helmholtz projector for the non-solenoidal part is the operator that maps j into Σj . Since the looked for projector is where + denotes the Moore-Penrose pseudoinverse. The projector for the solenoidal plus harmonic components can be obtained out of complementarity as The same reasoning for dual functions leads to the dual definitions of the projector which is the non-solenoidal projector for dual functions. The solenoidal plus harmonic projector for dual functions is, again, obtained by complementarity as It is important to note that, even though the projectors presented so far include a pseudo-inverse in their definition, they can be applied to arbitrary vectors in quasi-linear complexity by leveraging algebraic multigrid preconditioning [62], [65], [66] and, as such, are fully compatible with standard fast solvers.
Quasi-Helmholtz projectors can be used to cure the different 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 projectors have been used to cure the low-frequency breakdowns of several formulations, however for the sake of readability 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 where, following a frequency analysis similar to the one used for loop-star/tree decompositions, an optimal coefficient choice can be found to be ∝ − 1 2 and ∝ 1 2 , that is,  resulting in a new system of equations where P y = j . The constant 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 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 where we used that P Σ T Φ, P Σ = T Φ, and P Λ T Φ, = T Φ, P Λ = 0 and thus lim where is a frequency independent constant. This approach can be proved to simultaneously solve the problem of catastrophic round-off errors in both the current and the righthand 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.

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-infrequency, 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 nonsolenoidal basis, for example, the conditioning of a loopstar 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 scalarvalued functions. While they yield an overall improved condition number with respect to classical loop-star and looptree 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].

B. 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 extracted , where the static contribution is extracted. We obtain this by replacing e i ˆ · with e i ˆ · − 1, whereˆ denotes the direction of propagation. Then 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 ˆ · − 1 should be replaced by a Taylor series, where the static part is omitted. Similarly, the far-field cannot be computed by simply evaluating ∫

=1
[j ] where j = P y from (72),ˆ = /| |. On the one hand, by computing the unknown vector of the unpreconditioned formulation j = P y , the different asymptotic behavior in of the quasi-Helmholtz components of j as denoted in (57)-(60) would lead to a loss of the solenoidal/quasiharmonic 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 sol-qhar = / P ΛH y and j nsol = i √ P Σ 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 nsol to the far-field by evaluating and the contribution of j sol-qhar by Also for the near-field computation, the separation in j sol-qhar and in j 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.

C. Magnetic Field Integral Equation
The MFIE has, other than in the Green's function kernel, no explicit dependency on and should thus be expected to remain well-conditioned in frequency for → 0. Indeed, for simply-connected geometries , we have cond(M ) = O (1) when → 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 illconditioned system matrix [91]; in the following, we are going to show that cond (M ) ≥ / 2 for some constant ∈ 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 − has genus , then the space ( − ) defined by the harmonic functions in − and the space ( + ) defined by the harmonic functions in + have both dimension . The space defined by ˆ P ( ) ≔ˆ × ( − ) | are the poloidal loops and the space defined by ˆ T ( ) ≔ˆ × ( + ) | are the toroidal loops [92], [93]. ˆ P ( ) has been show [89] to be the nullspace of M − 0 and ˆ T ( ) the nullspace of M + 0 and thatˆ ×ˆ T ∈ ˆ P ( ) andˆ ×ˆ P ∈ ˆ T ( ). 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 ˆ P ( ) or in ˆ T ( ), a consequence of the fact that these functions are not curl-conforming, as mentioned in Section IV. On the other hand, quasi-harmonic functions (and dual quasi-harmonic functions , respectively) are associated with the holes and handles of the geometry. They are not, unlike the locally defined loop functions , derived from a continuous scalar potential on . Together with the fact that and are solenoidal but not irrotational (since the RWG/BC functions are not curl-conforming), (32) implies that quasiharmonic loops are linear combinations of solenoidal and harmonic functions (i.e., quasi-harmonic functions are harmonic functions with solenoidal perturbation). Clearly, any quasiharmonic basis { } 2 =1 can be rearranged into two bases, where one basis { T } =1 is orthogonal to poloidal loops and the other basis { P } =1 is orthogonal to toroidal loops. In the following,ˆ T ∈ ˆ T ( ) denote the harmonic toroidal and P ∈ ˆ P ( ) denote the harmonic poloidal basis functions on , while T , P ∈ and T , P ∈ are their quasi-harmonic counterparts.
In [91], scalings were reported of the blocks of the system matrix of M + in terms of a quasi-Helmholtz decomposition An analogous result can be obtained for M − 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 → 0 of the block matrices, we start by considering the Taylor series of the Green's function kernel of the MFIE where = | − ′ | and = − ′ . An O ( 2 )-scaling is observed for a block matrix if the contribution due to the static term in the Taylor series vanishes. For as expansion function orˆ × 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 .
We now consider the static MFIE, which models the magnetostatic problem, where is either solenoidal or harmonic. In this case, we have This can be seen by considering that a solenoidal function has a corresponding scalar potential Furthermore, we have the equality Inserting (87) for the testing function in the right-hand side of (88) and using the fact that curl is the adjoint operator of curl , we obtain where ∓ T ≔ˆ ×(M ± 0 ) is the (rotated) tangential component of the magnetic field that is, − T is the (rotated) tangential component of when is approached from within − , and + T when is approach from within + . We recall that (90) is obtained by finding a vector potential such that = curl and noting that curl curl = under the assumption that div = 0 [40, see Chapter 6.1]. From and from [40, (6.17)] where curly braces {} mean that this part is evaluated only in ± and δ is the surface Dirac delta function, together with (90), we have Next we will establish that ˆ × P , M ± 0 = 0 and ˆ × T , M ± 0 = 0. First, note that the exterior MFIE operator has the mapping properties that [89] M + 0ˆ T = 0 and while for the interior MFIE operator, we have the mappings and Furthermore, for any two surface functions , , we have [89, Section 5] Then we find where we used (94)-(97) and the orthogonality of harmonic and irrotational functionsˆ × . Now consider that P = P + P and T =ˆ T + T , where P/T is the respective perturbation. Thus by taking into account (84), we have From now on, we will only consider the harmonic functionŝ T andˆ P instead of their quasi-harmonic counterpart since, as we have seen, the solenoidal pertubation will always vanish. Then for M + 0 , we have asˆ ×ˆ P ∈ ˆ T ( ) and this space is orthogonal to ˆ P ( ). Likewise, we can conclude that there is at least one such that Analogously, we obtain This has established the scalings of the different block matrices composing M ± .
To prove the ill-conditioning in , we consider the static limit → 0, where j = 0 and following the argumentation in [95]. Forming the matrix-vector product, we find Thus for → 0, we have the reduced system where is the dimensionality of the respective finite element space of the basis functions . Since = , this implies that M + 0 has an P = -dimensional nullspace. To obtain a lower bound on the condition number dependency of M ± on , we consider that for all

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 → 0, for example, we should observe P Σ M ± 0 P Λ = 0; this is not observed in practice due to errors, for example, in numerical integration. To obtain accurate solutions even for → 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 Σ M ± 0 P Λ to zero. On multiplyconnected 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, where is the counterpart of P 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 → 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 and its discretization as where we note that for small | − ′ | a Taylor-series should be used to remove the static contribution. Then the overall system matrix is which evidently does not contain the static kernel contribution

D. Handling of the right-hand side
Similarly to the case of EFIE, the right-hand side requires special treatment. First, we note where h i extracted is computed analogously to e i extracted . 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 K − will dominate. The second term, however, needs to be further stabilized. We note where we used P ΣH G −1 × , P Λ = 0 from Lemma 1. Using (81), we have Summarizing, we obtain for the right-hand side which can be rearranged into

1) Other Low-Frequency Regularizations of the MFIE
When it comes to the low-frequency regularization of the MFIE, the body of literature falls shorter compared with 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 ℎ [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 quasistatic 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 lowfrequency behavior, a word should be said regarding fast methods. Typically, a fast method is employed to accelerate the matrix-vector products. On 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 either 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. D -D S 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 ℎ → 0.
The condition number of the EFIE is known to grow as 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 space C and the domain and range of the respective operator is considered [104]. Based on inverse inequalities, the norm of a Sobolev space is bounded by the norm of , where < . In the case of the EFIE operator, this is approach needs a non-trivial extension as the −1/2 div -space necessitates a Helmholtz decomposition.
To avoid the introduction of either Sobolev space or pseudodifferential 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 T A has singular values that accumulate at zero with a rate of ℎ, while T Φ has singular values that grow to infinity with a rate of 1/ℎ resulting in an overall condition number scaling proportional to 1/ℎ 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 −1 , T . Figure 9 shows the spectrum of the EFIE operator discretized on a sphere, confirming that indeed the condition number scales as 1/ℎ 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 wellconditioned system matrices if the expansion and testing functions are 2 -stable (i.e., that the Gram matrix is wellconditioned) [106].
For the CFIE, the conditioning improves to 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 T A . Hence, only T Φ contributes to the growing condition number.

A. Calderón Preconditioning
The Calderón identity suggests that the EFIE operator can be turned into a secondkind integral operator, which, as other second-kind integral equation like the MFIE, is well-conditioned in ℎ. 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 where The preconditioning strategies has been used extensively and successfully in several application scenarios with effective improvements [111]- [113]. However, when the frequency is extremely low (when 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 −I/4 + K 2 = (I/2 + K ) (−I/2 + 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 The projectors ensure that the overall system matrix is wellconditioned 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 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] where the outer matrix P o is with where G , is the Gram matrix of piecewise linear dual basis functions as defined in [45] and patch basis functions with  patch height 1/ , where is the area of cell of the mesh; furthermore, the middle matrix is P m ≔ P mΛ + P mΣ (132) with the definitions where G and G 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 where the norms can be estimated with the power iteration method for computing the largest singular value. We note that for → 0 we have which are similar scalings as encounter 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 guaranties convergence, and has a lower computational cost with respect to CGS or GMRES [53]. Figure 13 shows the number of iterations for the refinementfree preconditioned compared with (128), a standard EFIE and a loop-tree preconditioned EFIE. It displays that the number of iterations become bounded independently from ℎ, this 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 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 Nyström 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 highfrequency 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 divinner 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 ℎ. Recently, a formulation has been proposed that yields an ℎ-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 Φ, -part of T [81]. Recently, a scheme based on generalized primal and dual Haar prewavelets has been presented, where the condition grows only logarithmically in , and whichin 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 lowfrequency 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 ℎ. 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. 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 nearfield 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 · 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 Operator preconditioning, in this logic, strives to find B as an approximation of A −1 .
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 a local in nature, instead, to consider global interactions, [152] proposes to use interpolative 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 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 nearzone preconditioning has been proposed [159].

VI. H -F I
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 in T needs to be made a complex number. One possibility for a complex is to use the Yukawa potential resulting in 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óntype 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 with a suitably chosen 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/ℎ. 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 ℎ, the condition number can be asymptotically bounded for ℎ → 0), the condition number cannot be bounded independently from . 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 increases and the mesh parameter ℎ 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 ℎ = const while → ∞.
It can be shown that in general both the standard CFIE and the Calderón-Yukawa CFIE have singular values that are unbounded in . As, however, the problem is that the singular values of the analytic operator are unbounded in , the condition number of the Calderón-Yukawa CFIE grows when → ∞. On spherical geometries, it has been proved that the condition number grows as O ( 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 ( 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 ( 1/3 ) and a minimum singular value in the non-solenoidal part of the spectrum that decreases as O ( −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 ( 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 m . This modified wavenumber m is the original wavenumber plus a shift in the complex plane. This modified wavenumber has the form [165] where is a dimensionless constant and is a length that depends on the geometry. For example, on spherical geometries, is the radius of the sphere and is often chosen to be ≈ 0.4, which is actually based on optimizations for the conditioning of the acoustic operators [166], 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 = 1/ .
Among the solutions found in the literature, [162] proposes a Calderón-like preconditioning with a modified wavenumber m in the preconditioner of the EFIE operator, then an unchanged MFIE is added in a CFIE fashion. The work [167] also uses a Calderón-like preconditioning of the MFIE with the same 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 [168], 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 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.
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 [169]). 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 discretized I, that is, G , : 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 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 [170]. 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 > REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER < 20 Fig. 15. Example of ill-shaped triangles.
isogeometric analysis. For the electric field integral equation such approaches have recently been presented [171]- [174]. Improvements of the Calderón preconditioning in (128) can be expected by placing an additional Gram matrix resulting in . It ensures that the condition number can be bounded independently from ℎ 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-homogeneous which occurs, for instance, in the case of meshes resulting from adaptive refinement. The cost for inverting G −T × , are 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 T are often more expensive than inverting G −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 [175]- [186].
Hierarchical higher-order basis functions should not be confused with hierarchical bases used for curing the densediscretization 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 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 has been reported, in contrast to a polynomial growth for hierarchical bases [181, p. 187].
Strategies to overcome the -breakdown rely on constructing orthogonal higher-order basis functions. In fact, hierarchical higher-order basis functions enforce orthogonality between functions on the same level [178], [179]. 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 [190]. As a consequence from operator preconditioning theory [110], [191], it is possible to obtain a higher-order discretization of the EFIE, where the condition number is bounded independently from ℎ and . 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 .
VIII. C 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 different strategies.
is a sphere of radius centered at the origin. The scalar spherical harmonics (SH) are noted with ≥ 0 and | | < . The vector spherical harmonics (VSH) are defined as They form an orthonormal basis for vector fields on with respect to the weighted inner product Note that and are tangential, while are radial vector fields. Therefore, only and are used in the analysis of the surface integral operators. By construction, this basis is Helmholtz decomposed: are solenoidal and are irrotational.
In the following, J is the Riccati-Bessel function and H is the Riccati-Hankel function of first kind. They are related to the Bessel function and the Hankel function of first kind by J ( ) = π /2 J +1/2 ( ) and H ( ) = π /2 +1/2 ( ). Also, ′ denotes their derivative with respect to their argument. We thus obtain the analytic expressions of the EFIE operator [108], [192], [193] the exterior MFIE operator and the interior MFIE operator The tangential vector fields and are eigenvectors of the MFIE operator, and they are singular vectors of the EFIE operators. Assume we discretize the operator T with the VSH, then it is clear that the singular values associated to the right singular vectors are , = | J ( ) H ( )|, and the singular values associated to the right singular vectors are , = | J ′ ( ) H ′ ( )|. 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 ∇ · = 0, the study of the operators applied to gives the behavior of the operators applied to solenoidal functions. Whereas the study of the operators applied to gives the behavior of the operators applied to nonsolenoidal functions. The following asymptotic expressions can be obtained easily from the asymptotic behavior of the special functions in different regimes [194].

A. Low-Frequency Breakdown and other Issues 1) Low-Frequency Breakdown
We study the condition number of the discretized EFIE for → 0, where we keep the number of basis functions constant. In this case, it corresponds to considering only the orders less than some constant max . Thus in the asymptotic expansion is fixed and independent of . In the following, we will use the asymptotic behavior of the Riccati-Bessel, Riccati-Hankel, and their derivatives for small arguments given by where !! is the double factorial (or semifactorial) that is defined by the recurrence relation !! = ( − 2)!! and 1!! = 0!! = 1. Therefore, as → 0, the singular values scale as As a consequence, the condition number grows as O ( 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 −1 and by scaling the non-solenoidal part by . 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 (T ) 2 , whose singular values are as → 0. As a result, the condition number of the Calderón EFIE remains bounded as → 0.

2) Low-Frequency Round-Off Errors in the Excitation
We consider the plane wave pw ( ) = 0 e i ˆ · ˆ . Its tangential component expands as [195] ,±1 , pw = 2πi 0 2 + 1 4π ,±1 , pw = ±2πi 0 2 + 1 4π For → 0, the dominant terms of the solenoidal and nonsolenoidal parts are 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 T = −ˆ × i when the excitation is the plane wave i ( ) = 0 e i ˆ · ˆ . From (149) and (162), the solution has the following expansion For → 0, the dominant terms of the solenoidal and nonsolenoidal parts are respectively 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 ℎ decreases and the frequency remains constant. A decrease in ℎ 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 . The number of VSH whose order is less or equal to grows quadratically with , due to their multiplicity on the order (| | ≤ ). When we use local basis functions such as RWG to discretize the operators, the number of unknowns also grows quadratically with /ℎ. Therefore, to study the dense-discretization breakdown, we can use the correspondence = O ( /ℎ). In the following, we will use the asymptotic behavior of the Riccati-Bessel, Riccati-Hankel, and their derivatives for large order given by Therefore, as → +∞, the singular values scale as As a consequence, the condition number grows as O (ℎ −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 → +∞ (or equivalently as ℎ → 0), the singular values of (T ) 2 are As a result the condition number of the Calderón EFIE remains bounded as ℎ → 0.

C. High-Frequency Resonances
We study the condition number for ≫ 1. Here the point to notice is that is kept constant as → +∞, and 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 , or , 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 ≈ 2.74. Also, it is clear from Figure 16 that the higher 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 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 − )ˆ ×, where 0 < < 1 is a dimensionless constant. The CFIE is The singular values of this CFIE operator associated to the right singular vector and are, respectively, which are never 0 since the zeros of J and J ′ never coincide (due to the fact that the zeros of the Bessel functions are simple with an exception at = 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 resonances 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 → ∞ while keeping ℎ constant. It is discussed in the next section.

D. High-Frequency Breakdown
In this section we study the asymptotic behavior of the singular values of the Calderón-Yukawa CFIE operator [96] −T i T + (I/2 − K i ) (I/2 + K ) .
Recall that the CFIE is used to avoid the interior resonances, while Calderón preconditioning avoids the dense-discretization breakdown. The Yukawa part denotes the use of pure imaginary wavenumbers in the preconditioning operators 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 1/3 , while the minimum singular value, which is associated with the non-solenoidal functions, scales at most as O ( −1/3 ). The result is a growth of the condition number faster than a constant times 2/3 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.
a) Minimum singular value: The singular value associated with the nonsolenoidal VSH are It is useful to treat the order as a real variable instead of an integer. We obtain an upper bound for the minimum singular values by, choosing as the value for which J ( ) = π /2 J +1/2 ( ) = 0, that is, for which is a zero of J +1/2 (J is the Bessel function and J is the Ricatti-Bessel function). For the first zero of J +1/2 , the asymptotic expansion is given by [194, eq. 10.21.40] ∼ + 1/3 + 1/2 , or, equivalently, when with = −2 −1/3 1 ≈ 1.856, where 1 ≈ −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 [194, eq. 10.20.7-9], we find that the upper bound for the minimum singular value is asymptotically π Ai ′ ( 1 )| Ai ′ ( 1 ) − i Bi ′ ( 1 )| 2 7/6 ( ) −1/3 ≈ 0.69( ) −1/3 , Similarly to the previous case, there is another constant such that the index of a lower bound for the maximum singular value is asymptotically Using the asymptotic expansion of the Bessel functions in the transition region [194, eq. 10.19.8-9] we find that the dominant term in , is proportional to Ai(−2 1/3 )|Ai(e −iπ/3 2 1/3 )|. So has to maximize this quantity, and a derivative of it shows that is solution of c) 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 , 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 . To solve this issue one can to construct a Calderón-CFIE, but instead of using the wavenumber in the preconditioner, one uses a modified wavenumber m that is shifted in the complex plane to avoid reintroducing the resonances resulting in [167] −T m T + (I/2 − K ) (I/2 + K ) .
This modified wavenumber for the sphere has the form [123], [162] m = + i 1/3 −2/3 (191) where ∈ R is a constant. A common choice that gives a good condition number is These singular values have been plot in Figure 19 for increasing frequencies. It is evident that the condition number is bounded.
A B I G Lemma 1. For the inverse mixed Gram matrix Gˆ × , , we have Proof. First, we define the transformation matrix Q = Λ H Σ , where the global loops are constructed such that Λ T H = 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 To reveal how G −1 × , acts on the different quasi-Helmholtz subspaces, we derive an explicit expression for Q T , where in the following we abbreviate Gˆ × , = G. First, we recall the block structure of which follows from the orthogonality of solenoidal (i.e., ) and irrotational functions (i.e.,ˆ × ) / harmonic functions with irrotational perturbation (i.e.,ˆ × ). Using the Schur complement, we have Recursively application of (199) to (198) and considering the definition of projectors, the lemma follows.