On a Low-Frequency and Refinement Stable PMCHWT Integral Equation Leveraging the Quasi-Helmholtz Projectors

Classical Poggio–Miller–Chan–Harrington–Wu–Tsai (PMCHWT) formulations for modeling radiation and scattering from penetrable objects suffer from ill-conditioning when the frequency is low or when the mesh density is high. The most effective techniques to solve these problems, unfortunately, either require the explicit detection of the so-called global loops of the structure, or suffer from numerical cancellation at extremely low frequency. In this contribution, a novel regularization method for the PMCHWT equation is proposed, which is based on the quasi-Helmholtz projectors. This method not only solves both the low frequency and the dense mesh ill-conditioning problems of the PMCHWT, but it is immune from low-frequency numerical cancellations and it does not require the detection of global loops. This is obtained by projecting the range space of the PMCHWT operator onto a dual basis, by rescaling the resulting quasi-Helmholtz components, by replicating the strategy in the dual space, and finally, by combining the primal and the dual equations in a Calderón-like fashion. Implementation-related treatments and details alternate the theoretical developments in order to maximize impact and practical applicability of the approach. Finally, numerical results corroborate the theory and show the effectiveness of the new schemes in real case scenarios.

The EFIE's low frequency ill conditioning problem (lowfrequency breakdown) can be solved by performing a loopstar decomposition and rescaling the components of both the basis and the testing functions [19]- [23]. Loop-star decompositions, however, cannot solve the EFIE's ill conditioning problem associated with densely discretized geometries (mesh refinement breakdown). Both breakdowns can instead be solved by Calderón or by hierarchical preconditioning (see [12], [24] and references therein). Some drawbacks, however, are associated with these techniques in their standard incarnations. Hierarchical bases require the detection of global topological loops for non simply connected geometries, an operation which is computationally cumbersome. Moreover, standard Calderón methods, when they are not combined with quasi-Helmholtz decompositions, suffer from numerical cancellations in the solution current [12]. At the same time, standard quasi-Helmholtz decompositions such as loopstar/tree lead to ill-conditioning when combined with Calderón preconditioning [25]- [27]. Recently [12], to solve this, quasi-Helmholtz projectors have been introduced that can be used to obtain a LF stable, well-conditioned, cancellation and globalloops detection free, integral equation for scattering problems involving perfect conductors.
The origin of PMCHWT equation's conditioning and cancellation problems is similar to that of the EFIE's [13], [14]. In fact, on simply connected smooth geometries, the PMCHWT operator can be regarded as a compact perturbation of a linear combination of EFIE operators. As a consequence, loop-star decompositions can be employed to solve the low frequency breakdown of the PMCHWT equation [13], [28]- [31]. For multiply connected geometries, the picture becomes more complicated. In contrast to the EFIE operator, the MFIE operator featuring in the PMCHWT equation does differentiate between local and global loop currents. This translates in the fact that the EFIE rescaling method fails to solve the low frequency breakdown of the PMCHWT equation when applied to multiply connected geometries. This problem has been recently analyzed and addressed by a hierarchical scheme [18] that still requires, however, the expensive detection of global loops. Summarizing, currently there is not a formulation available for penetrable objects which is numerically stable and well conditioned till extremely low frequencies, immune from mesh related breakdowns and that does not require the detection of global loops. This lamentable gap will be filled by this work.
The contribution of this paper is three-fold: (i) a new regularization strategy for the PMCHWT based on the quasi-Helmholtz projectors. This scheme fixes the deterioration of the system condition number and the loss of accuracy due to cancellation of the current components resulting from the lowfrequency breakdown. As opposed to what is the case when using other quasi-Helmholtz decomposition based techniques, it does not require the explicit construction of the global loops in the boundary element space used to approximate the current. The application of the quasi-Helmholtz Projector to the PMCHWT turned out to be far from trivial. In fact, unlike the EFIE case [12], where the quasi-Helmholtz components of the testing space are rescaled, the range space of the PMCHWT operator is first projected onto a dual basis, and then the quasi-Helmholtz components of this projected current are rescaled. This new regularization strategy transforms the PMCHWT in an equation with the correct frequency scaling, even for multiply connected geometries. (ii) The extension of the new rescaling to the dual space, where suitable dual elements and dual projectors are adopted. (iii) A suitable combination of (i) and (ii) to obtain a Calderón -like PMCHWT equation which, to the frequency benefits of the formulation obtained in (i), adds the fact of being immune from mesh-refinement ill-conditioning. Very preliminary results were included in the conference proceedings [32]. This paper is organized as follows. Section II fixes the notation and introduces some background material. Section III presents a brief analysis of the standard PMCHWT conditioning and numerical cancelation problems. In Section IV, the novel frequency rescaling method is introduced. Section V extends this rescaling method to the dual space. This dual formulation is used in Section VI to construct a Calderón preconditioner that solves the PMCHWT equation's mesh refinement breakdown. In Section VII, a number of implementation-related details are discussed to facilitate the practical applicability of the new schemes presented here. Numerical results are presented in Section VII. Finally, our conclusions are summarized in Section IX.
II. NOTATION AND BACKGROUND Consider a dielectric body Ω characterized by a permittivity ǫ ′ and a permeability µ ′ , which is embedded in a medium characterized by material parameters ǫ, µ. Scattering by Ω can be modeled using the PMCHWT equation where the integral operators T k and K k are defined as Furthermore, Γ is the boundary of Ω,n r is its outward normal at point r, k = ω √ µǫ, η k = µ/ǫ and η r = η k /η ′ k . The quantities with primes (k ′ , η ′ ) are defined similarly with material parameters ǫ ′ and µ ′ . In (1), J (r) and M (r) are the equivalent electric and magnetic current densities on Γ, from which the scattered electromagnetic fields can be computed when an incident electromagnetic wave E i (r), H i (r) impinges on Ω.
In order to construct a numerical solution to (1), the surface Γ is triangulated by a mesh of average mesh size h. Moreover, on the mesh, the unknown currents M (r) and J (r) are most commonly expanded in Rao-Wilton-Glisson (RWG) functions f i (r) [5] normalized such that the integrated flux through their defining edges equals 1 and (1) is tested with the rotated RWG functionsn r × f i (r), i = 1, 2, ..., N . The following matrix equation is obtained: It is well known [19]- [23] that the RWG basis admits a quasi-Helmholtz decomposition in terms of Loop and Star basis functions, i.e. the RWG coefficient vector can be decomposed as where the matrices Λ, Σ, and H are the RWG to Loop, Star, and discrete harmonic function coefficient matrices, respectively. The Loop basis functions are solenoidal (their divergence is zero) and are associated to an internal node of the mesh. Their associated transformation matrix is defined as: where the symbols used are defined in Fig. 1. The Star basis functions are associated to the cells of the mesh. Their associated transformation matrix is defined as: Finally, H is the so-called global loop to RWG transformation matrix. Global loops are the discretized counterparts of the harmonic space of the Helmholtz decomposition. For the sake of brevity, we will omit a complete description of this family of functions and the interested reader is referred to [19] and [33].
Here it is just important to recall that the column dimension of H is 2N handles , where N handles is the number of handles of Γ.

III. A BRIEF ANALYSIS OF PMCHWT CONDITIONING PROBLEMS AND NUMERICAL CANCELLATIONS FOR MULTIPLY CONNECTED GEOMETRIES
To properly justify the need for the new formulations that will follow in the next sections, we will now briefly analyze the pathologic behavior of the PMCHWT as pertains frequency, discretization, and numerical stability. This section is included for completeness, but will be very synthetic for the sake of conciseness. The reader interested in more details could refer to the more extended analysis in [18] done under the prospective of hierarchical bases.

1) Ill-Conditioning in Frequency:
The following frequency scaling for the EFIE and MFIE operators holds [12], [28], [34] Note that (17) is not optimal since more stringent bounds can be derived if the global loops are separated into poloidal and toroidal loops [35]. This, however, is not a limitation for the analysis that follows and thus this distinction is omitted here for the sake of simplicity.
After defining the Loop-Star transformation matrix the following frequency scaling is obtained for the transformed PMCHWT operator From this, the fact that the condition number of the PMCHWT will grow for decreasing frequencies follows directly from a simple application of the Gershgorin disk theorem [36].
2) Ill-Conditioning in Discretization: When applied to simply connected geometries, the behavior of the PMCHWT condition number as a function of the average mesh size h follows directly from the behavior of the EFIE operator T . In fact, the PMCHWT operator can be written as where the latter is a compact operator (on simply connected smooth geometries) since it is off-diagonal with compact blocks. Otherwise said, the PMCHWT is a compact perturbation of a block operator containing T on the diagonal blocks. The overall h-related behavior will, as a consequence, be the one of the EFIE operator T . In other words, a h −2 conditioning growth for h → 0 is expected. For multiply connected geometries, the same conclusion holds. It must however be noted that the MFIE operator K behaves as an identity on the finite-dimensional harmonic space.

3) Numerical Cancellations in the Solutions:
The low frequency studies on the MFIE operator [34], [35], [37] have shown that the condition number of the equation is not the only factor that should be considered in evaluating the numerical stability of a method. In fact, at very low frequency, numerical cancellations could occur in the solution current vector for both ill-conditioned and well-conditioned operators. We will now show that numerical cancellations occur also for the PMCHWT and that the new formulation that will follow will need to be able to address this issue together with the conditioning problems analyzed previously.
For a plane wave excitation, the excitation vector of the PMCHWT scales as [18] This, after a Sherman-Morrison analysis [18] leads to the following scaling for the solution current When computing the charge density and the scattered field, the star part of these currents is multiplied with ω −1 . However, for low frequencies, this component (which scales as O(ω)) falls below machine precision and undergoes numerical cancellation. This results in a wrong radiated field. This phenomenon is the PMCHWT counterpart of what it is known to occur for MFIEs and both preconditioned and unpreconditioned EFIEs [12], [37].
IV. ON THE USE OF QUASI-HELMHOLTZ PROJECTORS: THE FREQUENCY REGULARIZED EQUATION All the PMCHWT equation problems, delineated in the previous section, have a counterpart for the EFIE. For this equation, a solution was proposed in [12] obtained by introducing the quasi-Helmholtz projectors. For the RWG basis functions, the projectors are defined as where + denotes the pseudoinverse matrix. For the EFIE, both the low-frequency ill-conditioning and the numerical cancellations can be solved by rescaling the operator as follows [12] T where and where k 0 is a constant introduced to make the equation dimensionally consistent. In the numerical results section, k 0 will be chosen equal to 1 per meter, in order to obtain consistency with [12]. Other choices are equally valid and will not influence the conclusions of this paper. Given (26) and the structure of the PMCHWT operator, containing the EFIE operator on the diagonal blocks, the first intuitive choice to solve the problems of the PMCHWT would be to adopt the operator M in a block-wise manner, i.e. proposing the following Helmholtz decomposed operator Such an immediate and intuitive choice, however, would lead to a catastrophic consequences as pertains the low frequency behavior. In fact the frequency analysis in this case would read The above expression shows that an order ω −1 term in the off-diagonal blocks (corresponding to the interaction of the harmonic subspaces, associated to the matrix H, of the electric and magnetic currents) would results in an unstable condition number at low frequency. Otherwise said, the strategy working for the EFIE will not work for the PMCHWT when multiply connected geometries are involved, and therefore a different approach is required. This will be investigated in the remaining part of this section. The reason why the frequency behavior of (31) is unstable is that the regularization employed here is tailored to the diagonal blocks of the PMCHWT operator (the T operators), but fails to regularize the K operators in the off-diagonal blocks. More specifically, after regularization, the scaling of the harmonic components of the off-diagonal blocks is inversely proportional to ω. In order to develop a regularization strategy for the complete PMCHWT operator, rather than only for its diagonal blocks, it is necessary to leverage on the dual elements. Define The currents M (r) and J (r) are expanded in RWG basis functions, as in the previous sections. The auxiliary currents M ′ (r) and J ′ (r) are expanded in Buffa-Christiansen (BC) dual basis functions g i (r), defined on the dual mesh by a barycentric refinement [38] By testing (32) with the rotated RWG functionsn r × f i (r), the following equation is obtained where the Gram matrix G is defined as Next, define the auxiliary unknowns where with With this, we find that In the right hand side of (43), the range of the PMCHWT operator is first projected onto the dual space (spanned by BC functions), and then rescaled in frequency. This is fundamentally different from the EFIE rescaling (26), where the frequency rescaling is performed in the testing space rather than the range space. Finally, the same rescaling is performed on the right hand side of the PMCHWT equation, leading to the new quasi-Helmholtz projected PMCHWT equation (qHP-PMCHWT) which we propose here where The beneficial properties of the new equation (44) with respect to (43) can again be appreciated by a frequency analysis. To this purpose, it should be preliminarily noted that where denotes a nonzero block with O(1) frequency scaling. With this, we find the following frequency scaling for the qHP-PMCHWT operator It is clear that the condition number of this matrix remains constant in the low frequency limit, so that the new qHP-PMCHWT is immune from the low-frequency ill-conditioning. It is now possible to verify the absence of numerical cancellations. For a plane wave excitation, we find that which leads to a solution current scaling as Now the Σ component has the same frequency scaling as the Λ component. As a result, no numerical cancellations occur between these two terms. The H component, on the other hand, does get cancelled at extremely low frequency. This is in itself not a problem, because its contribution to the scattered field is also O(ω) relative to the contribution of the Λ component. Moreover, this cancellation is physical. The previous derivation assumed a plane wave excitation, which does not excite global loop current in the static limit. Other types of incident field discributions (e.g. Although the qHP-PMCHWT in (44) is immune from lowfrequency problems and from numerical cancellations, it still suffers from discretization related ill-conditioning. This will be solved via a Calderon strategy. Before and in order to do this, however, it is necessary to extend, by duality, the qHP-PMCHWT in the dual space. This will be the focus of next section.

V. THE FORMULATION IN THE DUAL SPACE
A good strategy to obtain dual discretization is to proceed to the systematic replacement of RWG functions, with BC functions and of standard quasi-Helmholtz projectors with dual ones. Otherwise said, the dual formulation is obtained via the replacements This results in the following dual qHP-PMCHWT where Similarly to the primal case (analyzed in the previous section), the dual equation has a stable frequency behavior. This can be directly inferred by the quasi-Helmholtz decomposition analysis below The existence of the dual formulation, will allow to build up a Calderon preconditioned qHP-PMCHWT equation which, in addition to the immunity to low frequency and numerical cancellation problems of the new equation presented in Section IV, will have the additional feature of being immune from the discretization related ill-conditioning.

VI. THE FREQUENCY AND REFINEMENT REGULARIZED EQUATION
After stabilizing in frequency the PMCHWT in the standard and dual mesh, a Calderón strategy will deliver an equation that will be stable also for increasing discretization densities. A Calderón strategy requires the preconditioning of the primal operator with the dual one. In our specific case the new Calderón qHP-PMCHWT equation we propose reads We will now proceed to study the properties of this new equation. the frequency stability and the absence of numerical cancellation in the solution current directly follow from the fact that both properties holds for the formulations in (44) and in (57). In order to study the behavior of the operator of the new equation for increasing discretization density (h → 0), it is useful to explicitly compute the blocks in the matrix product in (64) and perform some simplifications. This provides for the block (1, 1), for the block (1, 2), for the block (2, 1), and for the block (2, 2). After performing the simplifications, all the four blocks in (65)-(68) are in the form M −1 XM where X is the corresponding block of the standard Calderónpreconditioned PMCHWT matrix (see [13]). Since the original Calderón PMCHWT equation does not suffer from refinement ill-conditioning problems, the same behavior is transferred to the new equation we present here. In fact, matrices M and M −1 are linear combinations of projectors. As such, their spectrum is piecewise flat and it does not introduce a variation as a function of the mesh parameter h.
From the above discussion, one might be tempted to use M and M −1 as right and left preconditioners for the CP-PMCHWT equation, respectively, and to eliminate the inner rescaling (MM −1 ). This inner rescaling is however required in order to perform a number of cancellations, as will be elucidated in the next section.

VII. IMPLEMENTATION RELATED DETAILS
This section focuses on the practical details to implements, in a numerically stable way till arbitrary low frequencies, the new frequency and frequency-refinement stable integral equations in (44) and (64). The three required building blocks for both equations are the stable discretizations of the electric operator, of the magnetic operator and of the right-hand-sides. These three elements will be detailed below.
The projected primal electric terms must be implemented using the following decomposition where the matrix relationships T h P ΛH = 0 and P Λ G −1 T h = 0 have been used. The projected primal magnetic terms must be implemented using the following decomposition where K 0 is the discretized static MFIE operator matrix and K e = K − K 0 . The matrix relationship P Λ G −1 K 0 P ΛH = 0 has been used. This matrix relationship, when finite precision integration rules are used, is only approximately satisfied. This is why the above splitting between static and dynamic contributions of the magnetic operator is necessary. The explicit cancellation of P Λ G −1 K 0 P ΛH in such a decomposition, in fact, ensure the numerical stability of our scheme up to arbitrarily low frequences. The projected dual electric terms must be implemented using the following decomposition where the matrix relationships T h P ΣH = 0 and P Σ G −1 T h = 0 have been used.
The projected dual magnetic terms must be implemented using the following decomposition where K 0 is the discretized static MFIE operator matrix and K e = K − K 0 . The matrix relationship P Σ G −1 K 0 P ΣH = 0 has been used.
As it is customary in schemes that must stabilize the extremely low-frequency regime, also the right-hand-side of the equations must be handled with care. In practice, in the right hand side of both (44) and (64), the following decomposition must be implemented where h ex and e ex represent the right hand side vectors computed using the extracted kernel method, i.e. where the exponential of the Green's function is replaced by (e ikk·r − 1) as standardly done in low-frequency stable methods [12], [21].

VIII. NUMERICAL RESULTS
A first set of tests have been done on a dielectric sphere of unitary radius and dielectric constant ǫ r = 3. Figure 2 shows the condition number of the new equation (64) as a function of the frequency. The performance of the standard PMCHWT, of a Loop-Star preconditioned PMCHWT and of the formulation presented in this work are compared. It is clear that the formulation proposed here is immune from the lowfrequency breakdown. Moreover, although both Loop-Star and the new formulation reach a constant condition number at lowfrequencies, the constant condition number reached by the new formulation (cond=1.8) is much lower than the value reached by Loop-Star (cond=63890) as predicted by the theory.
The fact that the new formulation is also immune from the mesh refinement breakdown is shown in Figure 3 where the condition number is depicted against 1/h. It is clear that differently from Loop-Star techniques the new formulation has a constant condition number for increasing mesh densities.
To check both immunities also on a non-simply connected geometry, this first set of tests has been repeated for a torus of external radius equal to 1.5m, internal radius equal to 0.5m and of dielectric constant ǫ r = 3. For this case the fact that the new equation is immune from the low frequency breakdown is confirmed by Figure 4, while the mesh refinement breakdown immunity is verified in Figure 5.
In order to show that the new scheme is also immune from the low-frequency current cancellation problem, the Loop and Star components of the current (plane wave exitation) at extremely low frequencies are benchmarked in Figures 6 and 7. It is clear that the new formulation provides the exact current. This is further verified in Figure 8 for which the correct far field is obtained differently from what can be achieved by a formulation that does not adopt an Helmholtz decomposed current. The error in the non-Helmholtz decomposed scheme, is due to the very low-frequency current cancellation analyzed in Section III-3. The fact that the new scheme is immune from this problem, which we have predicted theoretically, is clearly confirmed in practice by this numerical test. Finally, in order to show the accuracy of the new formulation also with near field sources, the current and the far field due to a dipole source are shown in Figures 9, 10, and 11.
To test the new formulation on a more complex geometry we have used a model of a quadccopter frame structure which has a non simply connected geometry comprising five global loops (refer to Figure 12). The electric and magnetic current for plane wave scattering at 10 −40 Hz are shown in Figure 12. The convergence behavior of our new formulation is shown in Figure 13. The new formulation solved the problem in only 5 iterations which compared very favorably with other schemes.

IX. CONCLUSIONS
This work has presented a new integral equation formulation applicable to scattering by penetrable media, based on the PMCHWT equation. Differently from all previously proposed approaches, this new formulation is concurrently immune from ill-conditioning when the frequency is low or when the mesh density is high; from numerical cancellations at very lowfrequencies; and from the necessity of detecting the global loops of the structure. This result has been achieved by a new projector-based strategy that proved effective for both simply and non-simply connected geometries. All theoretical