A Low-Frequency Stable, Excitation Agnostic Discretization of the Right-Hand Side for the Electric Field Integral Equation on Multiply-Connected Geometries

In order to accurately compute scattered and radiated fields in the presence of arbitrary excitations, a low-frequency stable discretization of the right-hand side (RHS) of a quasi-Helmholtz preconditioned electric field integral equation (EFIE) on multiply-connected geometries is introduced, which avoids an ad hoc extraction of the static contribution of the RHS when tested with solenoidal functions. To obtain an excitation agnostic approach, our approach generalizes a technique to multiply-connected geometries where the testing of the RHS with loop functions is replaced by a testing of the normal component of the magnetic field with a scalar function. To this end, we leverage orientable global loop functions that are formed by a chain of Rao–Wilton–Glisson (RWG) functions around the holes and handles of the geometry, for which we introduce cap surfaces that allow to uniquely define a suitable scalar function. We show that this approach works with open and closed, orientable, and nonorientable geometries. The numerical results demonstrate the effectiveness of this approach.


I. Introduction
T HE electric field integral equation (EFIE) is widely used to model electromagnetic scattering and radiation problems.It is well known, however, that the EFIE becomes unstable in the low-frequency regime: when the frequency decreases, the condition number of the system matrix resulting from its discretization increases [1]- [4].Historically, this lowfrequency breakdown has been cured with explicit quasi-Helmholtz decompositions of the surface current density such as the loop-star or the loop-tree basis [1], [3]- [7].More recently, an implicit decomposition based on quasi-Helmholtz projectors has been presented [8]: unlike loop-star or looptree functions, quasi-Helmholtz projectors maintain the 2stability of Rao-Wilton-Glisson (RWG) [9] functions; thus, no additional ill-conditioning with respect to the average edge length ℎ of the triangulation is introduced [7].
However, for the accurate computation of surface currents and fields, a well-conditioned system matrix alone is not sufficient.One issue that must be resolved is that testing the incident field directly with solenoidal and, on multiply-connected geometries, with solenoidal functions associated with the handles and holes of the structure can lead to catastrophic roundoff errors.In consequence, the physically correct scaling in frequency of the right-hand side (RHS) cannot be maintained resulting in incorrect surface current densities [10].
One strategy to obtain an accurate discretization of the RHS leverages a Taylor series expansion to set the static contribution to zero when tested with solenoidal functions [8], [11], [12].Notably, this approach is independent of the topology of the underlying geometry.Since it is not agnostic of the RHS, the extraction of the static part must be derived and implemented for each excitation.For an impinging wave that is modeled based on measurement data, however, this might be impractical.
In the case that information on the magnetic field is available, an alternative solution is to replace the testing of the electric field with a solenoidal function by testing the normal component of the magnetic field with a corresponding scalar function.One of the earliest accounts of this method is by Mautz and Harrington [2], [13], where they used a decomposition similar to the loop-star decomposition, but adapted to a body of revolution formulation.Moreover, the method was used in the context of magnetostatics [14], [15], and later applied to the RWG discretized, loop-star decomposed EFIE [3].In these works, the treatment of multiply-connected geometries was only partial.For example, [3] notes that global loops must be incorporated into the quasi-Helmholtz basis, but the discretization of the RHS is only studied for local loops, for which a corresponding scalar function is straightforward to define.In [2], open surfaces with holes are treated and it is suggested to introduce a cap surface in order to close the hole.However, while the analysis was general, the basis functions provided were limited to bodies of revolution.Closed bodies, such as a torus, or non-orientable surfaces were not treated.Hence it is unclear, how such an approach can be extended to global loops on surfaces of arbitrary topology.
In this work, we generalize the approach of Mautz and IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION Harrington to the case of an RWG-based discretization, where the underlying geometry may be open or closed, simply-or multiply-connected, and orientable or non-orientable, thus, enabling the treatment of arbitrary excitations in an RHS agnostic approach.Specifically, the contribution is two-fold: i) We show how a scalar function can be efficiently derived for global RWG loops which are constructed as a chain of RWG functions forming an orientable strip.This allows to introduce a cap surface such that a scalar function can be defined.The obtained scalar function is then used to test the normal component of the magnetic field resulting in a stable testing procedure for arbitrary RHS excitations.ii) We demonstrate how this generalized excitation agnostic RHS scheme can be combined with the quasi-Helmholtz projectors of [8], thus, facilitating their use with arbitrary excitations.Numerical results corroborate the presented theory and demonstrate the effectiveness of our approach, including cases of non-orientable and multiplyconnected geometries such as the Möbius strip.Note, that some preliminary results were presented in [16]- [18].
To this end, this article is organized as follows: Section II introduces the background material including the corresponding quasi-Helmholtz decomposition and fixing the notation.Section III elaborates on the numerically accurate testing scheme for arbitrary incident waves, including an analysis for the root causes of the occurring round-off errors.The necessary adaption of the projector based approach is deduced in Section IV, and numerical studies on the influence of the RHS evaluation on the solution accuracy are presented in Section VI.

II. Quasi-Helmholtz Decomposed Electric Field
Integral Equations Consider a perfectly electrically conducting (PEC) object described by a bounded, two-dimensional Lipschitz manifold ⊂ R 3 which can be open or closed, simply or multiply connected, where we assume that is orientable in the case that is closed.1The object shall be embedded in a homogeneous background medium with permittivity and permeability .A time-harmonic excitation field ( ex , ex ) induces a surface current density on from which the scattered or radiated fields can be computed.The current density satisfies the EFIE [21] ( where with wavenumber = √ , angular frequency , imaginary unit j 2 = −1, and an assumed but suppressed time dependency of e j ; moreover, an implicit normalization of the current with respect to the wave impedance is assumed.In order 1A Möbius strip is an example of an open, multiply-connected, and nonorientable .As noted in [19], [20], the EFIE can be solved on such . to determine , the surface is approximated by a triangular mesh and is expanded with RWG functions as ≈ =1 [j ] ( ), where j ∈ C contains the expansion coefficients.Employing a Galerkin scheme by using as testing functions, the linear system of equations is obtained with the matrix T ∈ C × exhibiting the entries for more details on the discretization see, e.g., [9]) and the RHS e ex ∈ C exhibiting the entries The first approach, which we consider here to address the low-frequency breakdown of the EFIE, is the loop-star decomposition.This decomposition will be used to analyze and to highlight different low-frequency issues related to the excitations.In a second step, we will then show how to adapt the solution strategies such that they become applicable to the quasi-Helmholtz projectors-the actual decomposition we are interested in.
As is well known [3]- [5], the basic idea of the loopstar decomposition is to express the current density as a superposition of local loops and global loops (associated with the handles and holes of ) representing the solenoidal components of , as well as stars representing the non-solenoidal component of resulting in the decomposition j = Λj Λ + Hj H + Σj Σ (5) with the loop to RWG expansion coefficient mapping matrix Λ ∈ R × , the star to RWG expansion coefficient mapping matrix Σ ∈ R × , both as defined in [7], and the global loop to RWG expansion coefficient mapping matrix H ∈ R × .As further detailed in Section III-A, we consider in the following matrices H that are constructed such that the corresponding global loops are formed by a chain of RWG functions with orientable support.Applying the transformation matrix where j = Qj ′′ , allows to remove the ill-conditioning of (3) with respect to by introducing suitable diagonal block normalization matrices D 1 and D 2 yielding the stabilized system with j = QD 2 j ′ .In literature, many choices for D 1 and D 2 have been presented [2], [8], [22]- [26]; we employ the normalization coefficients proposed in [25], [26].

III. Stable Evaluation of the Right-Hand Side On Multiply-Connected Geometries
In the following, we are interested in stabilizing the evaluation of the RHS of (6), which contains the products Λ T e ex , H T e ex , and Σ T e ex .We note that the -th entry of Λ T e ex corresponds to testing the incident field ex with the -th loop function , that is, The same is true for the -th entry of H T e ex , where the corresponding loop is a global one.As indicated in the introduction, evaluating (7) by numerical integration, either by adding up the contributions from each test function or by directly testing with , can lead to catastrophic round-off errors, the precise behavior depending on the RHS as will be discussed in the following subsection.Remedies are known for that avoid ad-hoc extractions of the static part of the RHS.In the following, we generalize these to , after an analysis for which RHSs we expect round-off errors.

As local loops
and global loops will be treated in a similar way, we call their union ℓ .That is, ℓ can represent or .Since all ℓ are solenoidal, i.e., ∇ • ℓ = 0 with ∇ • denoting the surface divergence, they can be expressed by the surface gradient of a scalar function2 as [27], [28, (A1.65) While for general this can lead to multi-valued [29], we consider here only which are constructed as a chain of RWG functions forming an orientable surface with surface normal ˆ (note that a surface normal for the whole surface is neither needed for the EFIE nor introduced, as for nonorientable surfaces there is no canonical normal vector-field on [19], [20]).An example of such a chain is shown in Fig. 1.The construction of global loops is a common problem in many branches of finite and boundary element method problems, see, for example [30].Specifically, algorithms [31]- [36] are suitable for finding global loops on a surface.3 Naturally, this poses the question if for any for which the EFIE is admissible we can find such an orientable surface.
3The reader should be reminded that the found loops are not uniquely determined.Even when the sequence is optimized such that the shortest possible loops are determined, the found global loops are not unique.E.g., for a torus several poloidal loops of equal length can exist.As noted in Section II, we limit ourselves to manifolds (e.g., this excludes geometries with T-junctions) in order to leverage results from topology, though we conjecture that the method would work even for more general geometries.More specifically, invoking the classification theorem for surfaces [37, p. 181], the topological constraints on denoted in Section II imply that is either homeomorphic to a sphere, to the connected sum of tori, or to a non-orientable Möbius strip, where each of these may have a finite number of holes cut out.It turns out that on such , we can always find with the desired properties: for homeomorphic to a sphere or to the connected sum of tori, this follows from the fact that they are orientable surfaces; for homeomorphic to a Möbius strip, we recall that when the strip is cut along its center line, the resulting strip is no longer a Möbius strip but a longer orientable strip with two half twists [37, pp. 162ff].On this orientable strip, which is homeomorphic to a sphere with two holes, we can construct a global loop with the desired properties.For example, the green triangles shown in Fig. 2 form an orientable surface that could serve as support for a global loop (evidently, the Möbius strip must be meshed finely enough so that there is a center line).Hence, no additional restrictions on are imposed.
Constructing the in such a manner (as chain of RWG functions) ensures that the boundary S of the support S of any consists of two disjoint curves C ,0 and C ,1 (see Fig. 1).As will be shown in Section III-B, it allows to choose the corresponding to be zero on one of the boundary curves, and it ensures to be single-valued by implicitly considering one of the curves as a cut, where otherwise would be discontinuous.
Inserting ( 8) into (7) and using the identity of [13, (C-1)], the representation is obtained.Here the line integral is along the boundary S of the support S of ℓ , and ˆ is the unit vector tangent to S forming a right-handed system with ˆ .For , the contribution of the line integral is zero as = 0 on ∂S .For , we have = 0 only on one of the two disjoint curves forming ∂S , for example, C ,0 , whereas we have = 1 on the other curve C ,1 .Thus, the line integral does not vanish on C ,1 .To remove this line integral, we introduce an arbitrary cap surface S with ∂ S = C ,1 (see Fig. 1), apply Stokes' IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION theorem with respect to S , and leverage that is constant on This representation shows that testing the incident field ex with a loop function (i.e., when evaluating ( 7)) corresponds to computing the curl of ex perpendicular to ˆ weighted with a scalar function (see also the remark for local loops in [3]).Such a computation is stable whenever |∇ × ex | and | ex | have a similar order of magnitude.Otherwise, the limitation of the numerical precision can lead to a loss of significant digits resulting in an erroneous RHS.In general, |∇ × ex | and | ex | will asymptotically scale differently in for → 0, and thus result in unstable discretizations of the RHS.However, it is interesting to note that if is satisfied (i.e., |∇ × ex | and | ex | have the same asymptotic scaling in ) and if a sufficient machine precision is used, then the testing with loop functions is stable independent from (and thus, also the evaluation in ( 7)).To obtain a numerical test that detects whether or not this condition is satisfied, we recast (11) by leveraging Faraday's law ∇ × ex = j ex resulting in the equivalent condition Again, it should be stressed that, in general, an excitation will not satisfy this condition and thus its discretization will become unstable for → 0. This will be confirmed by the numerical results in Section V.

B. Handling Arbitrary Excitations-Construction of Scalar
Functions for Global RWG Loops Equation ( 10) can be leveraged to obtain a stable RHS evaluation scheme for : We forgo the implicit computation of ∇ × ex by using its analytical expression ∇ × ex = j ex (as was originally suggested in [2] in the context of simplyconnected body of revolution problems) resulting in To derive an expression for corresponding to , we first consider the well-known local loop case.Here, the scalar function can be expressed in closed form on the -th triangle as [22], [38]  which corresponds to a pyramid shape around the center node as shown in Fig. 3.The vectors 1 and 2 are defined per triangle with area .Inserting ( 14) into (13), the products Λ T e ex can be evaluated as with = [Λ T ] 1 the number of RWGs spanning , since S vanishes.
To derive for a global RWG loop , we exemplify the discussion by considering the with support S around the hole shown in Fig. 4 (a).First, we introduce a vertex at an arbitrary point (e.g., in the center of the hole) such that a triangulated cap surface is introduced (see Fig. 4 (b)) corresponding to in Section III-A. 4 The global loop can then be expressed by introducing local loops of equal strength around each vertex inside the global loop.Figure 4 (b) shows that the overall current density of the global loop is left unaltered, since all other contributions cancel each other.In consequence, is the superposition of the scalar functions of the local loops such that is constant on the parts of lying inside C ,1 (the inner part of S ) and on the cap surface, which is illustrated in Fig. 4 (c).Hence, the support of partially lies outside of .This corresponds well with the expression derived in (10), where by applying Stokes' theorem a cap surface was introduced.
Expressing through (14), using = 1 on the cap surface elements S , the stable evaluation scheme for global loops is obtained with As a consequence, the numerically stable testing requires the evaluation of the incident fields across the holes and handles of a multiply-connected surface .
It should be stressed that the support of does not necessarily have to touch , that is, it is not required to find the shortest global loop.Also, the cap surface does not necessarily have to be introduced over the inner boundary (in the example C ,1 ), but also the outer boundary (C ,0 ) can be chosen.As noted before, this approach can also be applied to multiply-connected surfaces which are closed, such as a torus, where both toroidal and poloidal loops can be treated in the same manner.
It is also interesting to note that the introduction of cap surfaces strikes some resemblance to the introduction of socalled cutting surfaces in the context of eddy current problems using a magneto-quasistatic approximation (see, e.g., [39], [40]).In the latter case, the cutting surfaces are needed to render the involved potentials single-valued as a fundamental requirement to obtain correct solutions [39]- [44].In contrast to the cap surfaces employed in this work, however, the introduced surfaces can have a finite thickness, they cut the domain instead of closing handles and holes, the potential itself is part of the unknowns, and the setting is commonly based on a tetrahedral discretization of the volume [45], [46].

IV. Stable Right-Hand Side Evaluation for
Quasi-Helmholtz Projectors Instead of directly using decomposition (5), it is advantageous to utilize the quasi-Helmholtz projectors introduced in [8], in particular, when combined with the (refinement-free) Calderón multiplicative preconditioner [47], even though the 4Alternatively, several vertices can be introduced or one of the vertices on one of the boundary curves can be selected and connected to all other vertices of the same boundary curve.
required search for global loops in our scheme loses one of the advantages of quasi-Helmholtz projectors, which do not require this search.One could think of avoiding the construction of global loops; however, introducing a cap surface that closes the handles and holes requires a cut in the surface.Finding such a cut corresponds to finding a global loop.
The fundamental idea of the quasi-Helmholtz projectors is to decompose the surface current into solenoidal and nonsolenoidal components by employing projectors, that is, j = P ΛH j + P Σ j = j sol,hsol + j nsol (17) with the orthogonal projectors P ΛH ∈ R × and P Σ ∈ R × defined by P Σ = Σ (Σ T Σ) + Σ T and where (•) + denotes the Moore-Penrose pseudoinverse.Defining the decomposition operator with and as defined in [25], [26, ( 26)-( 28)] and applying it to the discretized EFIE T j = e ex in (3) as yields a stable linear system of equations (LSE) (which is better conditioned than ( 6)) with j = P 2 j ′ .Yet, analogous issues for the stable evaluation of the RHS arise and we need to establish how ( 15) and ( 16) can be used in this context.

A. Stable Evaluation of the Right-Hand Side
A direct evaluation of the RHS P 1 e ex in (20) as results in similar round-off errors as for the loop-star decomposition.To see why, we construct P ΛH explicitly.This can be done by applying [Λ H] T to (5) yielding where the orthogonality relations Λ T Σ = 0 and H T Σ = 0 were used.After solving for j Λ and j H and mapping back to Λj Λ and Hj H , the projector is obtained.In consequence, evaluating P ΛH e ex involves the products Λ T e ex and H T e ex exhibiting the round-off errors discussed in the previous section.At the same time, using (23) instead of (18) for the RHS product allows to employ the stable evaluation schemes given in ( 15) and ( 16) for the quasi-Helmholtz projectors rendering the product P ΛH e ex numerically stable.However, evaluating the pseudoinverse in ( 23) involves matrices Λ T Λ which exhibit a condition number which grows with the average edge length ℎ of the triangulation as 1/ℎ 2 [7].Hence, an efficient inversion strategy is required.

B. Efficient Evaluation of the Projected RHS
As Λ T Λ is a graph Laplacian, it can be inverted with a constant number of iterations independent of ℎ by employing an algebraic multigrid (AMG) preconditioner [7], [48].In order to maintain this property for (23), we propose to employ the Schur complement (see, e.g., [49, p. 650]).More precisely, is evaluated, where y Λ and y H are the solutions to These solutions are obtained by first solving the LSE containing the Schur complement with P Λ = Λ(Λ T Λ) + Λ T and subsequently solving This allows to use an AMG preconditioner each time an inverse of Λ T Λ is formed, and thus, fully removes the illconditioning with respect to ℎ in (28).We do not expect an ℎ-ill-conditioning of S in (26) due to the global nature of the .As an example consider the rectangular torus shown in Fig. 5 for which 1/ℎ is varied between 2 and 100.The corresponding number of iterations for conjugate gradient (CG) to compute the pseudoinverse when evaluating P ΛH e ex is depicted at the top of Fig. 6.Clearly, if no preconditioner is employed, the iteration count increases with 1/ℎ.This is improved by applying an AMG preconditioner directly to the whole matrix in (25).However, a detrimental dependency on ℎ remains.Only when computing the pseudoinverse via the Schur complement, the number of iterations (summed up for both LSEs to be solved) is constant and overall the smallest for all considered cases.
As the pseudoinverse depends also on H, it is of interest to investigate the influence of the number of global loops .To this end, we vary the number of handles of the plate with finite thickness depicted in Fig. 7, where there are two global loops per handle.The corresponding number of iterations is shown at the bottom of Fig. 6.Again, the overall fewest iterations are required when the proposed Schur complement scheme is used.The remaining dependency on for the LSE in ( 26) is acceptable in most realistic cases, since the involved matrix is only of size × and the LSE has to be solved only once in the overall solution process.
In addition, it should be noted that the proposed approach is fully compatible with a Calderón preconditioner.The latter is commonly combined with the stabilized system (20) to prevent a dense-discretization breakdown, that is, an ill-conditioning of the system matrix with respect to ℎ [8], [50].As this does neither affect the problems with the RHS (the same kind of RHS stabilization is needed) nor the proposed solution strategy all benefits of the Calderón preconditioning are maintained.

V. Numerical Results
In order to show the impact of the stabilized RHS evaluation on the accuracy of the computed scattered and radiated fields, several scenarios are investigated.For all of them, a GMRES solver is employed without restarts and a relative residual of res = 1 × 10 −6 as stopping criterion.Furthermore, the quasi-Helmholtz projectors are employed in all scenarios.The implementation of (3) is based on [51].

A. Scattering from a Sphere
We begin with a study of the scattering from a sphere since a semi-analytic reference is available.The sphere has radius s = 1 m and is excited by a plane wave, a Hertzian dipole HD [52, pp. 69 f.], a magnetic dipole FD [52, pp. 76 ff.], an electric ring current rc [52, pp. 362 ff.], and a magnetic ring current rc (employing duality) as depicted in Fig. 8.For all of these excitations, the scattered fields can be determined semianalytically in form of a series expansion (see [52, pp. 368 ff.] and [53] for implementation details), which serves as reference solution.The error with respect to this solution is determined for the far field (FF), the electric, and the magnetic near field (NF) at a distance of = 5 m.To this end, the fields are computed on a spherical 5 • grid and the relative logarithmic worst case error with ∈ , ℎ, FF is determined.
For the plane-wave excitation, the results over frequency are depicted at the top of Fig. 9 (b) for the cases where P Λ e ex is evaluated in a naïve manner and where P Λ e ex is evaluated via the scalar functions .Clearly, for the naïve RHS evaluation, all fields are erroneous when going below 1 × 10 −4 Hz.Using the RHS evaluation via solves this for all fields.The second plot of Fig. 9 (b) shows the norms of the solenoidal components of the RHS e sol,hsol = P ΛH e ex and the non-solenoidal components of the RHS e nsol = P Σ e ex .From these it can be seen that the physically correct scalings (see Table I which generalizes the findings of [10]) are only obtained for the evaluation via .Otherwise a deviation is observed below 1 × 10 −8 Hz, showing that the scalings reflect the breakdown in accuracy.The same is true for the norms of the solenoidal and non-solenoidal components of the solution current j as shown in the third plot of Fig. 9 (b).
For the Hertzian dipole, it can be seen from Fig. 9 (c) that the breakdown in accuracy occurs already in the kHz region.This can be explained by the fact that the solenoidal and the non-solenoidal component of the RHS differ in their scaling by 2 (see Table I), which is in contrast to the plane wave, where the difference is only .Consequently, the naïve RHS evaluation breaks down earlier.Evaluating the RHS via leads to correct FFs and NFs over the whole frequency range.Notably, other than the FF and the magnetic NF, the electric NF is determined correctly also with the naïve RHS evaluation.The reason is that the non-solenoidal current component, which represents the charge and thus determines the electric NF, dominates asymptotically and is thus recovered sufficiently accurately (see third plot of Fig. 9 (c)).At the same time, for this specific excitation the electric NF is dominated by the non-solenoidal current in the low-frequency regime as can be shown by an asymptotic analysis [26,Table II].One should not jump, however, to the conclusion that the electric NF is always recovered accurately.In fact, for the plane-wave excitation both current components are recovered erroneous below 1×10 −8 Hz leading to all fields being erroneous as well.
A similar behavior is observed for the magnetic ring current as depicted in Fig. 9 (d).As it behaves like a Hertzian dipole in the low-frequency regime, it shows the same breakdown frequencies.Only the scalings of the RHS components are different, however, the ratio between solenoidal and nonsolenoidal components stays the same.
The results for the electric ring current shown in Fig. 9 (e) are particularly insightful.In Section III-A, we predicted that if O (| |) = O (| ˆ • | ), then the -based evaluation of the RHS is unnecessary.We observe that, indeed, there is no further improvement in the error and the physically correct scalings for the RHS are recovered in both approaches.As an interesting observation, we note that this is directly reflected by the scalings of the RHS components: Considering (10) together with Faraday's law and using e sol,hsol = P ΛH e ex , we find O (| ˆ • ex | ) = O (∥e sol,hsol ∥).Analogously, with e nsol = P Σ e ex , we get O (| ex |) = O (∥e nsol ∥).Hence, condition (12) can be expressed equivalently as showing that no numerical round-off errors occur if the nonsolenoidal component scales as the solenoidal components of the RHS.As the asymptotic behavior of the magnetic dipole differs from the electric ring current only by a constant factor, we omit it here.
As another class of excitations, we consider spherical vector waves travelling towards the origin.Again the scattered fields can be determined analytically by enforcing the boundary conditions providing a reference solution.Following the definitions in [54, pp. 325ff], the incident fields are the transverse electric (TE) modes (3)   1 and the transverse magnetic (TM) modes (4)     1) ) as long as the scatterer is not a sphere placed in the origin.In the latter case, we have e TM,sol,hsol = 0 as the surface curl of the TM modes vanishes in Hz (g) spherical mode TM 11 Fig. 9. Scattering from a sphere: worst case errors of the FF, the electric, and the magnetic NF, as well as the scaling of the RHS components for different excitations with and without stabilized RHS.

Table I
Scaling of the RHS components for → 0. plane wave Hertzian dipole el.ring current mag.dipole mag.ring current TE TM on the surface of a sphere [55].The results for a TM 11 excitation in Fig. 9 (g) show that without the stabilization scheme the error level increases already below a frequency of 10 MHz.Notably, despite e TM,sol,hsol not fully vanishing with the stabilized RHS, it is kept small enough such that the fields are computed correctly.This property is maintained if (1)   2 waves (with spherical Bessel instead of Hankel function in the radial dependency) are used as TM excitation modes.The increase of the error at = 1 × 10 8 Hz, on the other hand, stems from the (geometrical) approximation error.Increasing the number of triangles from 624 to only 1526 restores an error level of about −80 dB.

B. Scattering from Multiply-Connected Geometries
As a first example of a multiply-connected geometry, we consider the double torus depicted in Fig. 10.The results for a plane-wave excitation are shown in Fig. 11 (a).The reference solution is determined by stabilizing the RHS with a Taylor series expansion following [12].The comparison with the stabilization via scalar functions shows good agreement over the whole frequency range, clearly demonstrating the correctness of the proposed scalar function formulation for multiply-connected geometries.The same applies for a Hertzian dipole excitation as depicted in Fig. 11 (b).For the latter, the error is determined in form of a manufactured solution [56], [57]: The dipole is placed inside the scatterer, and we leverage that the total field sc + ex = 0 outside the scatterer.The latter condition is then checked on the spherical 5 • grid.Note, that this is not the same approach as in [58], where the current density itself is compared to a manufactured solution.The obtained error ℵ shows that the stabilization of the RHS via scalar functions yields accurate solutions.The strong relevance of the stabilized RHS evaluation in this case can be explained by the close proximity of the excitation to the scatterer.
As a more challenging scatterer, the geometric model of a (PEC) Fokker Dr.I shown in Fig. 12 with = 390 global loops excited by a Hertzian dipole is considered.In order to accelerate the computation, an adaptive cross approximation (ACA) algorithm is used [59], [60] with compression rate 1 × 10 −4 .The error of the solution is determined in the same manner as for the double torus.The results in Fig. 13 prove again the necessity of the proposed scheme.To demonstrate that also topologically challenging geometries can be successfully handled by the proposed approach, the trefoil knot depicted in Fig. 14 is excited by a Hertzian dipole.The hole of the toroidal loop is closed by introducing a point in the center of the knot.Consequently, the resulting surface is intersecting itself as shown in Fig. 15.While in principle this can be avoided by introducing a more complex (Seifert) surface [61], the results in Fig. 16 confirm that the correct fields are obtained.Figure 17 confirms that also the non-orientable Möbius strip depicted in Fig. 2 can be handled.Again, we use as reference solution a Taylor series expansion for the plane-wave excitation.Similar to the trefoil knot, the cap surface closing the hole is self-intersecting as shown in Fig. 18.

VI. Conclusion
We showed how to evaluate the RHS of the quasi-Helmholtz decomposed EFIE in an excitation agnostic manner such that no catastrophic round-off errors occur in the presence of a multiply-connected geometry due to the testing of the incident field with global loops.The numerical results show that, depending on the scattering object and the specific excitation, without the adaptions inaccurate fields are obtained already for frequencies below the MHz region, but with the stabilized evaluation the scattered and radiated FFs can be determined correctly down to the static limit.

S C , 1 C , 0 SFig. 1 .
Fig. 1.Example of a global loop as a chain of RWG functions with support S and boundary ∂S = C ,0 ∪ C ,1 .

Fig. 2 .
Fig.2.Mesh of a non-orientable Möbius strip discretized with 852 triangles and 1180 RWGs.The green triangles form an orientable surface that could serve as support for a global loop.

2 1 Fig. 3 .Fig. 4 .
Fig. 3. Representation of a loop formed by RWG functions around the center node by the scalar function .

Fig. 6 .Fig. 7 .
Fig. 6.Number of iterations for the iterative computation of the pseudoinverse in the projector P ΛH via CG.

Fig. 8 .
Fig. 8. Scattering from a sphere of radius s = 1 m (a) excited by dipoles and ring currents at position = 2 m and radius rc = 0.5 m; (b) discretized with 624 triangles and 936 RWGs.

Fig. 18 .
Fig. 18.Mesh to close the global loop of the Moebius strip.