An Excitation-Aware and Self-Adaptive Frequency Normalization for Low-Frequency Stabilized Electric Field Integral Equation Formulations

The accurate solution of quasi-Helmholtz decomposed electric field integral equations (EFIEs) in the presence of arbitrary excitations is addressed: Depending on the specific excitation, the quasi-Helmholtz components of the induced current density do not have the same asymptotic scaling in frequency, and thus, the current components are solved for with, in general, different relative accuracies. In order to ensure the same asymptotic scaling, we propose a frequency normalization scheme of quasi-Helmholtz decomposed EFIEs, which adapts itself to the excitation and which is valid irrespective of the specific excitation and irrespective of the underlying topology of the structure. Specifically, neither an ad hoc adaption nor a priori information about the excitation is needed as the scaling factors are derived based on the norms of the right-hand side (RHS) components and the frequency. Numerical results corroborate the presented theory and show the effectiveness of our approach.


I. INTRODUCTION
T HE electric field integral equation (EFIE) constitutes a flexible and accurate formulation for electromagnetic radiation and scattering problems. Preconditioning the discretized equation based on quasi-Helmholtz decompositions is among the most mature approaches to cure an otherwise occurring low-frequency breakdown [1], [2], [3], [4] (i.e., the condition number increases with decreasing frequency) and enables accurate solutions also for low frequencies down to the static limit. For example, this has been demonstrated for the approach based on quasi-Helmholtz projectors proposed in [5], an enhancement to the explicit decompositions based on a loop-star or loop-tree basis [1], [2], [3], [4], [6], [7], [8]. Manuscript  Yet, ensuring that the system matrix is well-conditioned is insufficient if one desires to accurately compute the fields. Two more aspects must be carefully addressed. First, testing the right-hand side (RHS) with solenoidal functions, which is needed for the considered decompositions, can lead to catastrophic round-off errors. Different means to overcome this issue have been proposed [9], [10], [11], [12], and in the following, it is assumed that a corresponding stabilization is applied since all further derivations rely on an accurately discretized RHS reflecting the physically correct asymptotic scalings in frequency of the quasi-Helmholtz components.
The second problem is addressed in this contribution. It has to be ensured that all quasi-Helmholtz components of the decomposed surface current density, that is, the solenoidal and the nonsolenoidal component, are recovered accurately enough to obtain accurate scattered and radiated fields [5], [13]. More precisely, the unknown vector of the linear system of equations (LSE) to be solved has two contributions (i.e., the solenoidal and the nonsolenoidal quasi-Helmholtz component), which can differ largely in their order of magnitude. This is due to the, in general, different asymptotic scaling of the quasi-Helmholtz components in the wavenumber k as k → 0, where the specific behavior is dependent on the excitation source. Still, both contributions need to be determined accurately. 1 As the employed preconditioners entail a rescaling of the quasi-Helmholtz decomposed expansion and basis functions, they enforce a different scaling of the quasi-Helmholtz components of the preconditioned current (i.e., the unknown vector) compared to the physical (nonrescaled) one. The same is true for the RHS components. Different choices for the corresponding scaling coefficients have been employed in the past. Typical choices are, for example, given in [2], [5], [19], [20], and [21]. While they are suitable for plane-wave excitations, for arbitrary excitations, it is in general not ensured that all components can be recovered as the choices in these existing schemes do not reflect the scaling of the quasi-Helmholtz components of the excitation [22]. Yet, excitations, such as voltage gaps, are commonly used, and thus, the EFIE should be solved accurately for these as well. The same is true for spherical waves, which are able to represent the field radiated by antennas, or line currents representing radiation from wire structures.
A remedy is using ad hoc techniques as done in [5]. However, it constitutes an approximation, which limits the overall accuracy, and the behavior of the quasi-Helmholtz components of the excitation needs to be known a priori. Moreover, in [5], only three types of excitations were studied (plane wave, capacitive gap excitation, and inductive gap excitation). For other excitations such as spherical TM or TE waves, we find, however, that the approach is not applicable.
In this work, we propose a frequency normalization scheme, which adapts itself to the excitation requiring neither a priori information nor ad hoc adaptions. To this end, the scaling coefficients incorporate the norms of the quasi-Helmholtz components of the discretized RHS in a black-box-like manner. Hence, the approach is agnostic to the specific excitation (and the topology of the structure). More precisely, the scaling coefficients are determined such that all quasi-Helmholtz components of the preconditioned current and the RHS have the same asymptotic scaling. Thereby, all components are recovered with a similar relative accuracy, ensuring that all radiated and scattered fields can be computed accurately. Numerical results corroborate the effectiveness of our approach both for canonical and complex geometries. Note that some preliminary results were presented in [22] and [23].
This article is organized as follows. Section II introduces background material and fixes the notation. In Section III, the implications of the frequency normalization in quasi-Helmholtz decompositions of the EFIE are analyzed; from the gained insights, the adaptive normalization scheme is derived for the loop-star decomposition as well as for the quasi-Helmholtz projectors. In Section IV, we show how to achieve positive eigenvalues of the preconditioned system matrix, again for both decompositions, and implementationspecific aspects are addressed. Numerical studies are presented in Section V.

II. QUASI-HELMHOLTZ DECOMPOSED EFIES
Let a perfectly electrically conducting (PEC) object with surface Γ be embedded in a homogeneous background medium with permittivity ε and permeability µ. The structure is excited by a time-harmonic field (e ex , h ex ) resulting in an induced surface current density j on Γ from which the scattered or radiated fields can be computed. The current density satisfies the EFIE [24] (T j ) tan = e ex tan (1) where 2 T j = jkT A j + jk −1 T Φ j contains the vector potential operator 2 Note that the choice of defining without a surface normal vectorn is independent of the considered low-frequency stabilization and the proposed scheme. However, it has the advantage that nonorientable surfaces (which do not possess a canonically definable normal vector field) can be considered as scattering objects in a natural manner [25], [26]. and the scalar potential operator with the free-space Green's function the wavenumber k = ω √ µε, the angular frequency ω, and the imaginary unit j 2 = −1. Moreover, an implicit normalization of the current with respect to the wave impedance is assumed as well as a suppressed time dependency of e jωt . To solve for j , the surface Γ is triangulated and j is expanded with Rao-Wilton-Glisson (RWG) functions f n as j ≈ N n=1 [j ] n f n (r), where j ∈ C N contains the unknown expansion coefficients. Employing f n as testing functions in a Galerkin scheme results in the LSE where the matrix n dS(r) (to be evaluated in a weak sense), and the RHS vector e ex ∈ C N exhibits the entries [e ex ] m = Γ f m · e ex dS(r).
To address the low-frequency breakdown of the EFIE, we consider as the first approach the loop-star decomposition. It allows to carry out the analysis in a demonstrative way before adapting it to the second and actual decomposition of interest: quasi-Helmholtz projectors.

A. Loop-Star Decomposition
The well-known [3], [4], [6] basic idea of the loop-star decomposition is to express the current density j as a superposition of N Λ local loops Λ m and N H global loops H m (associated with the handles and holes of Γ ) representing the solenoidal component of j , as well as N Σ stars Σ m representing the nonsolenoidal component of j . To do so, the loop to RWG expansion coefficient mapping matrix Λ ∈ R N ×N Λ , the global loop to RWG expansion coefficient mapping matrix H ∈ R N ×N H , and the star to RWG expansion coefficient mapping matrix Σ ∈ R N ×N Σ all as defined in [8] are introduced. Applying the transformation matrix Q = [Λ H Σ] to (5) as T ΛHΣ j ′′ = Q T T Qj ′′ = Q T e ex , where j = Qj ′′ , allows to remove the ill-conditioning of (5) with respect to k by introducing suitable diagonal block normalization matrices D 1 and D 2 yielding the stabilized system with j = QD 2 j ′ . Common choices in the literature for these matrices are [5], or D 1 = diag(1/k, 1/k, 1) and D 2 = diag(1, 1, k) [19]. Here, we define it in a general form as and systematically derive the coefficients α i , γ i , and β i to suit arbitrary excitations. Consequently, the preconditioned which is related to the actual current vector as

B. Quasi-Helmholtz Projectors
As an improvement over the loop-star decomposition, the quasi-Helmholtz projectors introduced in [5] can be employed. The orthogonal projectors P Σ ∈ R N ×N defined by P Σ = Σ(Σ T Σ) + Σ T , where (·) + denotes the Moore-Penrose pseudoinverse, P Λ ∈ R N ×N defined by P Λ = Λ(Λ T Λ) + Λ T , and P H ∈ R N ×N defined by P H = I − P Λ − P Σ are introduced (note that the pseudoinverses can be computed efficiently via algebraic multigrid (AMG) preconditioning as detailed in [5, Section V]) as well as the decomposition operator Applying the latter to (5) as yields a well-conditioned LSE, which is better conditioned than (6), with if α i , β i , and γ i , i = 1, 2 are suitably chosen; specifically, in [5], a combined projector, P ΛH = P Λ + P H , is used with In contrast, we chose with our definition of P i to treat P Λ and P H separately as we need the general form for our theoretical apparatus. Our final formulation uses again the combined projector P ΛH .
For both decompositions, we define the quasi-Helmholtz components j ′ sol , j ′ hsol , and j ′ nsol , that is, the solenoidal component, the solenoidal component associated with the handles and holes of Γ , and the nonsolenoidal component, respectively, as shown in Table I. This leads to the preconditioner-agnostic definitions of j sol , j hsol , and j nsol . Moreover, we introduce j sol,hsol = j sol + j hsol . The RHS Helmholtz components e sol , e hsol , and e nsol , as well as their primed variants are defined analogously.
Moreover, note that in order to also cure the dense-discretization breakdown of the EFIE (i.e., the ill-conditioning with respect to the average edge length of the triangulation of Γ ), formulation (10) could be combined with a Calderón preconditioner similar to what has been done using the standard choice of scaling coefficients in [5], [27]. While such a combination should be compatible with the here proposed frequency normalization scheme, further investigations would be necessary, which are beyond the scope of this work.

III. ANALYSIS OF THE FREQUENCY NORMALIZATION
AND SELF-ADAPTIVE SCHEME Before deriving our new formulation, we start by investigating the influence of the different quasi-Helmholtz components of j on the scattered or radiated near field (NF) and far field (FF).

A. Relevant Current Components
To assess the influence of a current component on the FF, we utilize the stabilized evaluation [28] where T(·) denotes a Taylor series expansion around k · r = 0, and we divided by a normalization factor −jke −jkr /(4πr ).
Since T(e −jk·r − 1) = O(k), the combined current component j sol,hsol is scaled with an additional factor k. Adopting the Analogously, we can determine whether a current component is relevant for the computation of the electric or the magnetic NF. In the case of the electric NF, we use its stabilized evaluation which introduces a scaling of j sol,hsol with k and of j nsol with 1/k. In the case of the magnetic NF evaluation, both current components are weighed identically.
As an example, consider a plane-wave excitation, where we have ∥j sol,hsol ∥ = O(1) and ∥j nsol ∥ = O(k): here, only j nsol is relevant for the electric NF and only j sol,hsol is relevant for the magnetic NF, but both current components are relevant for the FF. Analogously, the relevant current components that are required to recover the NFs and the FF can be derived for the excitations shown in Table II, where we used the physical current scalings given in [29] and [5]. In addition, we have included the cases of the spherical TE mn and TM mn modes (corresponding to the definitions in [30, p. 325ff] with n ∈ {1, 2, . . . } and |m| ≤ n) based on the physically correct scalings of the current components ∥j TE,sol,hsol ∥ = O(k −(n+1) ) and ∥j TE,nsol ∥ = O(k −(n−1) ) for TE modes as well as ∥j TM,sol,hsol ∥ = O(k −n ) and ∥j TM,nsol ∥ = O(k −n ) for TM modes. The asymptotic scalings can be derived from the RHS component scalings given in [10] and a Schur complement analysis of j ′′ = T −1 ΛHΣ Q T e ex (see also the following). Clearly, for all of the excitations in Table II, both j nsol and j sol,hsol must be recovered accurately if one desires to compute both FF and NF accurately. Even more, both current components have to be recovered with a similar relative accuracy as one component can solely determine the accuracy of a field even though its magnitude might be small compared to the other current component. To ensure accurate recovery, we study, in a first step, the frequency scalings of the loop-star decomposed system (6).

B. Self-Adaptive Normalization for Loop-Star Basis
As will be shown in the following, while for a planewave excitation, the normalization matrices in (6) can be chosen as , or D 1 = diag(1/k, 1/k, 1) and D 2 = diag(1, 1, k) [19], for arbitrary excitations, this yields, in general, wrong solutions; instead, a more flexible approach is required. For the analysis, we use D i from (7), which allows for an asymmetric normalization. Performing a Schur-complement analysis for the inverse of D 1 T ΛHΣ D 2 for k → 0, we obtain since the multiplication with the mapping matrices Λ, H, and Σ does not change the scaling in k. The blocks of the matrix (D 1 T ΛHΣ D 2 ) −1 exhibit the scaling which is a direct generalization of the findings in [29]. In order to obtain a well-conditioned matrix for k → 0, we enforce that the blocks on the main diagonal scale as O(1); thus, the coefficients must obey and Inserting (16) and (17) into (15) and forming the matrix-vector product in (14) result in which shows that, dependent on the scalings of the incident field, the Helmholtz components of j ′ will differ in their asymptotic behavior in k.
While the standard choices (16) and (17), the resulting quasi-Helmholtz components of j ′ will, in general, not have the same asymptotic scaling. Consequently, a large relative error in one of the Helmholtz components remains undetected if, for example, the relative residual error of D 1 Q T e ex −D 1 T ΛHΣ D 2 j ′ is studied (as done by typical iterative solvers for the LSE). As noted, however, in Section III-A, all current components must be recovered with a similar relative accuracy when solving (6) in order to accurately compute all fields. These considerations lead us to impose a second requirement on α i , γ i , and β i : all current components should scale identically, where we, due to the constraints of finite-precision arithmetic, opt for an O(1)-type of scaling.
With the aim of satisfying this additional requirement, we first note that the scalings derived in (18) show that the solenoidal components associated with the local and the ones associated with the global loops can be treated identically. Consequently, in all of the following, we set: and define Moreover, we introduce the ratio of the norms of the solenoidal components over the nonsolenoidal component Excluding the trivial case ∥e sol,hsol ∥ = ∥e nsol ∥ = 0, the scalings in (18) then dictate Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. and Note that the case w > O(k 2 ) includes ∥e sol,hsol ∥ = 0 when interpreting it as the result of w → 0 and the case w < O(1) includes the case ∥e nsol ∥ = 0 when interpreting it as the result of w → ∞. Moreover, definition (20) covers also the cases where O(∥e sol ∥) ̸ = O(∥e hsol ∥), including the cases where not both but either e sol or e hsol vanishes (e.g., as encountered for gap excitations). Together with (16) and (17) for the remaining coefficients, we find α 2 = 1/(kα 1 ) and β 2 = k/β 1 (24) arriving at a scheme that adapts itself to the incident field based on the norms of the RHS components. Hence, there is no need to have a priori information about the nature of the incident field.

C. Quasi-Helmholtz Projector Normalization Factors
When employing the quasi-Helmholtz projector stabilized system (10) instead of the loop-star decomposition, the very same normalization factors (22)-(24) can be employed 3 with the corresponding definitions for e sol,hsol and e nsol . This follows from the intrinsic relation to the loop-star decomposition, that is, . However, it should be verified that α i and β i in (9) properly take care of another aspect only relevant for the projectors: since solenoidal and nonsolenoidal components are stored in the same vector, that is, j ′ = j ′ sol,hsol +j ′ nsol and P 1 e ex = e ′ sol,hsol +e ′ nsol , significant digits can be lost if the asymptotic scaling in k of j ′ sol,hsol and j ′ nsol (or e ′ sol,hsol and e ′ nsol ) is different. This is precisely the reason why an imaginary constant was included in the definition of the projectors in [5], where P 1 = P 2 = P ΛH / √ k +j √ kP Σ was employed. Due to the imaginary constant, for example, for the inductive gap excitation, the dominant contribution of j ′ sol,hsol was stored in the imaginary part and the dominant contribution of j ′ nsol was stored in the real part of j ′ , where we consider a component as dominant if its asymptotic behavior O(k a ) < O(k b ) with k b as the scaling for the other component. 4 This approach relies on the circumstance that real and imaginary parts of each quasi-Helmholtz component often exhibit a different asymptotic scaling so that, in fact, only one of the parts is relevant for an accurate field computation in lowfrequency scenarios. However, using the factors α i and β i as defined in (22)- (24) for P i , the inclusion of an imaginary constant is in the first place no longer necessary.
In order to illustrate this problem of the significant digits in more detail, we consider the physical frequency scalings of the RHS for seven different types of excitations. Those can be derived as shown in Appendix A, where the findings of [29] are generalized, resulting in the first three rows of Table III. Note that in contrast to the analysis so far, we study now real and imaginary parts separately to show that indeed no dominant contribution is lost. The scalings show that directly storing the components in one vector would lead to a loss of significant digits for a Hertzian dipole or a magnetic ring current excitation since the solenoidal components scale differently than the nonsolenoidal components. The next three rows show the scalings using the standard normalization from [5]. Here, the real and the imaginary part of the nonsolenoidal components are interchanged, which allows to store the components in one vector for all the considered excitations. However, for the excitation by spherical TE and TM modes, this approach fails, as shown in Table IV. Both real and imaginary parts have the same asymptotic scaling such that none can be neglected. At the same time, the here proposed normalization shown in the last rows (Tables III and IV) 1), and thus, they can be added without the danger of losing significant digits.
To complete the picture, we study also the scalings of the solution current components in Table V: the first three rows show the physical scalings as obtained by separately handling real and imaginary parts in (14). Evidently, storing the physical current components in one vector fails for the electric ring current, the magnetic (Fitzgerald) dipole, and the inductive gap excitation. From j = j sol,hsol + j nsol = P 2 j ′ = α 2 j ′ sol,hsol + β 2 j ′ nsol , it can be seen that the solution current of the preconditioned LSE scales as j ′ sol,hsol = j sol,hsol /α 2 and j ′ nsol = j nsol /β 2 .
Using the normalization from [5], rows 2-6 in Table V show that, at least for the considered excitations, all current components can be stored in one vector. However, again, for the TE mn and TM mn modes, this fails as, for example, ∥j TE,sol,hsol ∥ ⋆ = O(k −(n+1) ) + jO(k −(n+1) ) and ∥j TE,nsol ∥ ⋆ = O(k −(n−1) ) + jO(k −(n−1) ). Even more, for all the excitations (except for the plane wave), the dominant components of the current still scale differently, hinting that they cannot be recovered with a similar relative accuracy in a straightforward manner, but an ad hoc technique is required. This can be avoided when the proposed normalization scheme is employed as shown in the last three rows.
As an example of which current components can be recovered, we consider the double torus in Fig. 1 excited by a magnetic diople. The real and imaginary parts of the determined current components are shown in Fig. 2(a) for the standard normalization and in Fig. 2(b) for the proposed normalization. In order to ensure a high accuracy, a direct solver is employed. While the standard normalization, at first glance, seems to approximate the theoretical scalings (indicated by dashed lines) better, only the adaptive normalization maintains the correct scalings of the dominant components for all frequencies. This is, however, the decisive property. The worse agreement of the nondominant components is not relevant for the result-   ing fields as the nondominant components are (in this case 13) orders of magnitude smaller than the dominant components. A similar behavior can be observed for an excitation by an electric ring current as shown in Fig. 2(c) and (d), which confirms, again, the preceding analysis.

IV. EIGENVALUE DISTRIBUTION AND IMPLEMENTATION DETAILS
The definition of the normalization constants α i and β i in (22)-(24) allows for two more improvements leading typically to faster convergence of iterative solvers. By reintroducing an imaginary constant j and a constant C (to be defined for the loop-star basis and the quasi-Helmholtz projectors in the following), we ultimately propose to choose the normalization factors: where the simplified conditions for the case distinctions (not involving w) are sufficient for practical purposes as will be detailed in Section IV-A. Together with (16) and (17), we then find for the remaining coefficients Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  (19) is employed for γ i in the context of the loopstar decomposition. The motivation for these modifications will again be discussed first for the loop-star decomposition and then for the quasi-Helmholtz projectors.

A. On the Simplification of the Conditions to Be Checked
The suggested simplified conditions for the case distinctions in (26) and (27) are based on the observation that except for the capacitive gap excitation for all standard excitations listed in Table III, the ratio in (21) satisfies as can be seen from the first three rows. The capacitive gap excitation, on the other hand, is accounted for by the case ∥e sol,hsol ∥ = 0. Even when combining arbitrary TE mn and TM mn modes with arbitrary orders n, the relation (29) is satisfied with the sole exception when Γ is a sphere placed in the origin, all shown in Appendix B. In this particular case, ∥e TM,sol ∥ = 0. Consequently, if no other excitations than the ones listed are present, the case distinctions only require to check whether one of the norms is zero. For cases not covered by (29), the more general conditions on w given in (22) and (23) can be checked by evaluating the norms at least at two frequencies.

B. Loop-Star Decomposition
The idea of introducing the constant C in (26) and (27) is to improve the overall condition number of the system matrix.
Choosing it similar to the one suggested in [28] as enforces an equal contribution of the scalar and the vector potential. The matrix norms in this expression can be determined efficiently as ∥Υ ∥ = λ max (Υ H Υ ) for a square matrix Υ with λ max (Υ H Υ ) denoting the largest eigenvalue of Υ H Υ . This eigenvalue can be estimated, for example, via a few Arnoldi iterations or the power iteration method (which are both compatible with matrix-free methods) [31], [32]. The factor j, on the other hand, causes D 1 T ΛHΣ D 2 to have only positive eigenvalues in the static limit. For In order for C not to vanish in the static limit, ν ≤ 0 has to hold. For B not to vanish, on the other hand, ν ≥ 2 has to hold. Since both conditions cannot be satisfied at the same time, either C or B or both vanish for k → 0 independent of the RHS. For the case w > O(k 2 ), we find that B = O(1) and C = O(k 2 ), and for the case w < O(1), we find that B = O(k 2 ) and C = O(1). Thus, for all possible excitation scalings, D 1 T ΛHΣ D 2 is at least block triangular (if not block diagonal). In consequence, for the determinant, det(D 1 T ΛHΣ D 2 ) = det(A) det(D) holds, which shows that the eigenvalues are solely determined by A and D. As A and D are known to have positive eigenvalues [33], so does D 1 T ΛHΣ D 2 (for k → 0). This property is, in general, beneficial as it leads to faster convergence for typical Krylov subspace methods [34, p. 205ff].

C. Quasi-Helmholtz Projectors
Similarly, for the quasi-Helmholtz projectors, we introduce [28] to improve the overall condition number. Since the imaginary constants were demonstrated to be irrelevant for the significant digits in Section III-C (due to the new normalization scheme), imaginary constants can again be placed to improve the eigenvalue distribution. This leads to a positive eigenvalue in the static limit also for the quasi-Helmholtz projected EFIE, which can be shown by using the block triangular nature of D 1 T ΛHΣ D 2 for k → 0 as starting point. Defining the matrix ] and the matrix G i = GD i , the matrix G T 1 T G 2 is also block triangular with positive eigenvalues. Furthermore, we have P 1 T P 2 = GG T 1 T G 2 G T , which allows to deduce properties about the eigenvalues λ i of P 1 T P 2 . Due to the properties of similarity from which the positive eigenvalues of P 1 T P 2 follow for k → 0. In order to illustrate the impact, the eigenvalues λ i of the projected matrix P 1 T P 2 for the example of the double torus in Fig. 1 are shown in Fig. 3. Only if the imaginary constant is included, we have Re{λ i } > 0 and Im{λ i }, which are close to being zero. If the frequency is decreased further down to the kHz region as shown in Fig. 4, we see that the clustering on the real axis increases further. The influence on the number of iterations for the generalized minimum residual (GMRES) [35], the induced dimension reduction (IDR) [36], and the stabilized biconjugate gradient (BiCGstab) [37] algorithm are summarized in Table VI for the case that the double torus is excited by a Hertzian dipole and a relative residual of ϵ res = 1 × 10 −6 as a stopping criterion. Clearly, the iteration count (reflected by the number of matrix-vector products) is reduced if either the imaginary constant or the constant C is included. The best result is obtained if both constants are used simultaneously as is done in the proposed formulation.

V. NUMERICAL RESULTS
To demonstrate the impact of frequency normalization on the accuracy of the computed scattered and radiated fields, several scenarios are investigated. To this end, we employ the quasi-Helmholtz projectors and compare the standard choice  (22)  of scaling coefficients, α 1 = α 2 = √ C/k and β 1 = β 2 = j √ k/C, where we have included the constant C for a fairer comparison with the adaptive ones in (26) and (27). Again, if not stated otherwise, a GMRES solver is used without restarts and a relative residual of ϵ res = 1 × 10 −6 as a stopping criterion. For the RHS stabilization, the approach of [10] is employed, and in order to accelerate the computation, an adaptive cross approximation (ACA) algorithm is used [38], [39] with compression rate 1 × 10 −4 . The relative worst case errors with a ∈ {e, h, e FF } are computed based on a spherical 5 • grid and NF distances of r = 8 m.

A. Scattering From a Sphere
We start with scattering from the sphere shown in Fig. 5, where series expansions are used as a reference (see [40, pp. 368 ff.] and [41] for implementation details). The results for different excitations are shown in Fig. 6. For a plane-wave excitation, the standard normalization is sufficient to determine all fields for the whole frequency range correctly. Fig. 6(b) shows that the adaptive normalization maintains this property. This is as expected: for the plane wave, both schemes recover all current components with a similar accuracy.
Considering next the excitation by a Hertzian dipole, Fig. 6(c) evidences that the adaptive frequency normalization removes an otherwise occurring breakdown in accuracy of the magnetic NF below 1 MHz. Specifically, the magnetic NF is affected, which corresponds well with the preceding analysis (see especially Section III-A): as an incorrect solenoidal current component is recovered, only the magnetic NF is affected in accordance with Table II. Also depicted is the accuracy of an ad hoc adaption technique using the projectors P = P 1 = P 2 from [5] with α 1 = α 2 = √ C/k and β 1 = β 2 = j √ k/C. The ad hoc technique requires that two LSEs are solved. More precisely, as ∥Re{P T P }∥ ≪ ∥Im{P T P }∥ for low frequencies, the two real-valued systems −Im{P T P }Im{j ′ } = Re{P e ex } and Im{P T P }Re{j ′ } = Im{P e ex } are solved. Depending on the specific dominant components of the excitation (see Table V), only Re{j ′ } or Im{j ′ } is used to determine j sol,hsol and in an analog fashion to determine j nsol . While this approach also removes the breakdown below 1 MHz, it introduces a breakdown when going above that frequency.  A similar behavior is observed for the electric ring current as shown in Fig. 6(d) and the magnetic dipole in Fig. 6(e). The errors in the current components correspond to the fields, as predicted in Table II, and only the adaptive normalization (in addition to a stabilized RHS) fully removes the issues for all fields and for arbitrarily low frequencies.
The same is true for the spherical mode excitations TE 11 and TM 11 , as shown in Fig. 6(f) and (g), respectively. Note that the ad hoc approach is not possible in this case, but only the adaptive normalization can overcome the problem.
As the analysis so far showed that the problem lies in an insufficient accuracy of the solution components, another investigation is conducted for the sphere excited by a magnetic dipole: we vary the relative residual ϵ res of the GMRES solver. In this numerical experiment, the ACA is not used and we employ C = 1 (to emphasize the effect of the residual). The results for the standard normalization in Fig. 7(a) reveal that indeed by lowering the residual, the breakdown can be shifted to lower frequencies. However, even when employing a direct solver (LU decomposition), the breakdown cannot be fully avoided. For large problems, iterative solvers become the only option; a typical choice for ϵ res is around 1×10 −3 to 1×10 −4 . This underlines the importance for using an adaptive normalization as shown in Fig. 7(b), where we see that the error in the solution is dominated by the geometric approximation and discretization and not by the stopping criterion. Furthermore, the number of iterations remains constant when decreasing the frequency, whereas the iteration count for the standard normalization decreases with decreasing frequency: with the standard normalization, only the solenoidal component is recovered, but not the nonsolenoidal component. Hence, the iterative solver only determines a part of the unknowns sufficiently accurately but stops too early for the remaining unknowns resulting in the erroneous electric NF. Also, note that separate residuals for solenoidal and nonsolenoidal components are not possible (as the nonsolenoidal component of the RHS depends  on all current components), and if so, would lead to longer convergence times. Similarly, the real and imaginary parts of the RHS depend each on the real and the imaginary part of the current, thus preventing a separate residual for real and imaginary parts of the current without an approximation.

B. Scattering From Multiply-Connected Geometries
To validate our formulation also for multiply-connected geometries, we consider the double torus in Fig. 1 for a plane-wave and a Hertzian dipole excitation. The results in Fig. 8 demonstrate the effectiveness of our approach.
The results for the PEC model of a Fokker Dr.I shown in Fig. 9 excited by a Hertzian dipole in Fig. 10 highlight again that the breakdown in accuracy can occur at relatively high frequencies but is fully overcome by the proposed scheme.

C. Inductive and Capacitive Gap Excitation
As a last example, we consider the voltage gap excitation of the capacitive and the inductive structure shown in Fig. 11. The capacitive structure consists of two plates of size 1 m × 1 m separated by 0.01 m, where the nonuniformity of the mesh originates from the small size of the strip connecting the two plates. The inductive structure is a 1 m × 1 m ring with 0.25 m width. A comparison to the radiated fields when using the standard normalization is shown in Fig. 12. Clearly, for both structures, only the adaptive scheme yields the physically correct scalings of the current components, as given in Table V. For the inductive structure, the nonsolenoidal current is incorrect resulting in a difference in the electric NF. Dually, for the capacitive structure, the solenoidal current is incorrect resulting in a difference in the magnetic NF, all   [42]. Analogously, the determined capacitance of 918.05 pF agrees well with 917.65 pF from an electrostatic simulation in CST MWS.

VI. CONCLUSION
We have shown that scaling coefficients in the quasi-Helmholtz preconditioners, which lead to a well-conditioned impedance matrix, do not necessarily lead to current solutions that allow an accurate determination of all fields. In fact, despite a well-conditioned matrix, the current components are in general solved for with different relative accuracies, which is for general excitations not sufficient. The numerical results demonstrated that the proposed adaptive frequency normalization method effectively resolves this effect for all excitations considered and, thus, allows for the accurate computation of NF and FF. In fact, depending on the structure and the specific excitation, standard approaches lead to inaccurately computed fields already in the MHz region, whereas with the proposed method, the scattered and radiated fields can be determined correctly with an overall reduced number of iterations down to the static limit. Given our theoretical considerations and the wide range of studied excitations, we conclude that our method can stabilize the EFIE in the presence of arbitrary excitations.

APPENDIX A RHS SCALINGS
In order to derive the real and imaginary parts of the solenoidal and nonsolenoidal components of the RHS for the excitations in Table III (i.e., for k → 0), the corresponding excitation fields e ex are expressed via their currents j ex and m ex and Green's function representation Γ t · e ex dS(r) using a Taylor series expansion of Green's function with the distance R = |r − r ′ | between the source point r ′ and the observation point r and a Taylor series expansion of the gradient of Green's function with R = r −r ′ . We obtain the tabulated values by considering that the scalar potential contribution vanishes when tested by a solenoidal function, i.e., Γ t · ∇ ′ · j ex ∇G dV (r ′ ) dS(r) = 0 if ∇ Γ · t = 0 (39) that a constant vector tested by a solenoidal function vanishes, i.e., Γ t · c dS(r ′ ) = 0 if ∇ Γ · t = 0 and c = const (40) that a solenoidal excitation current tested by a solenoidal function vanishes, i.e., Γ t · j ex dV (r ′ ) dS(r) = 0 if ∇ Γ · t = 0 and ∇ · j ex = 0 (41) that [43] Γ t · R .
In case N TE = N TM , this leads to w = O(k), and for all N TE = N TM − g with g ≥ 1, it leads to w = O(k 2 ). In consequence, we have for all possible combinations of N TE and N TM . For the special case of scattering from a sphere placed in the origin, we have Again, for N TE ≥ N TM + 1, this leads to w = O(1). On the other hand, for N TE < N TM + 1, we have the case N TE = N TM leading to w = O(k) and the cases where N TE = N TM − g with g ≥ 1 resulting in which is a case not satisfying (47).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.