A Fast Converging Resonance-Free Global Multi-Trace Method for Scattering by Partially Coated Composite Structures

Global and local multitrace formulations (MTFs) provide a flexible and efficient method for the modeling of scattering and transmission of time harmonic electromagnetic (EM) waves by composite structures. This contribution extends the domain of applicability of global MTFs to cases where penetrable domains, perfectly conducting domain, and perfectly conducting thin sheets are all part of the design. A novel and easy to implement resonance-free equation to model scattering by perfect conductors is introduced. A Calderón preconditioner designed to limit the number of iterations required for the solution of the discrete system is designed and studied. The accuracy, flexibility, and efficiency of the method are demonstrated on a representative range of examples.


I. INTRODUCTION
T HE behavior of electromagnetic (EM) fields and waves in and around objects of arbitrary shape and varying material properties is a problem engineers are often tasked with. These objects can be fully metallic or dielectric, but are most commonly composite featuring multiple materials with dielectric and metallic regions. For arbitrary objects with piecewise homogeneous regions, the boundary element method (BEM) provides an attractive approach for the solution of the EM scattering and transmission problem that arises from these structures [1], [2].
The Poggio-Miller-Chan-Harrington-Wu-Tsai (PMCHWT) method is a boundary integral equation method well-known in literature for the solution of transmission problems with two or more subdomains [3]- [5]. It has been extended to deal with general metallic and dielectric structures including closed metallic objects [6]. The formulation presented in these papers remains susceptible to interior resonances. Moreover, PMCHWT is known to suffer from dense discretization breakdown [7], [8] which impairs the convergence of iterative solvers in the absence of effective preconditioning. For a single scatterer, or a collection of nontouching scatterers [9], Calderón preconditioning provides a highly efficient technique to suppress the growth in the number of iterations required for the solution of the systems resulting upon discretization of these integral equations. Unfortunately, these techniques cannot be applied in the presence of junctions, lines along which three or more domains meet. PMCHWT in the presence of junctions does not scale well with the number of degrees of freedom (DoFs). The problem becomes even more complex if both penetrable regions, thin conducting sheets, and perfectly conducting regions are present.
In search of a solution, a number of avenues have been explored. Nonoverlapping domain decomposition methods such as the boundary element tearing and interconnecting (BETI) method result in systems that remain amenable to preconditioning techniques and which can be solved in parallel to speed up the global solution. However, without further development, the BETI method for EM scattering does not scale well with the number of DoFs.
Single-trace solvers that can deal with junctions can be built if one is willing to part with conforming finite element spaces [10]. Unfortunately, no low-frequency or dense grid regularization schemes are known for this class of discretizations.
Global multitrace formulations (MTFs) for the solution of EM scattering and transmission problems have been independently introduced by Chu et al. [11], Ylä-Oijala et al. [12], and Claeys and Hiptmair [13]. The MTF starts from the intuitively appealing idea of introducing a small gap, occupied by the background material, between adjacent domains. Next, two pairs of Cauchy data are introduced at the interfaces between the domains that make up the structure. This is in contrast to single-trace formulations (STFs) such as PMCHWT which feature a single pair of traces at each interface. Even though the number of DoFs that needs to be solved for doubles, MTFs lead to more efficient solution strategies. The reason is that they remain amenable to operator or Calderón preconditioning schemes that are formally very similar to 0018-926X © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
those applied in the single-domain scenario. The introduction of a virtual gap separating adjacent domains conceptually removes any junctions, allowing for the formal application of preconditioning techniques originally developed for singledomain problems. Ylä-Oijala et al. [12] study a global MTF that can be used in the presence of both penetrable and perfectly conducting domains. The method is indirect, which means that the solutions of the integral equation cannot immediately be interpreted as the tangential traces of the electric field and magnetic field. It remains susceptible to interior resonances in the presence of perfectly conducting domains. The inclusion of thin sheets is not treated. It requires both a primal and dual finite element space for the discretization of the system without preconditioner. This precludes much needed acceleration methods such as the biparametric matrix assembly described in [14].
The analysis for the direct global MTF in the presence of both penetrable regions and perfect conducting regions has been presented in [13]. For acoustics, a combinedfield-type approach was suggested to allow incorporation of perfectly conducting regions without the introduction of resonance-related problems [15].
A good overview of both the local and global MTF can be found in [16] and [17]. Efficient implementations of the local MTF for composite structures comprising penetrable media can be found in [18]. The transmission conditions used in the local MTF have been further optimized in [19]. The algorithm was tailored to the modeling of close-to-resonance cavities in [20] and [21]. A generalized combined field integral equation (G-CFIE) MTF suitable for large multiscale problems was studied in [22]- [24]. Local MTFs are based on the use of field values of neighboring domains in the representation formulas of a given domain. Despite their conceptual simplicity, their analysis remains incomplete.
This article introduces three innovations to widen the domain of applicability to include realistic devices as encountered in engineering practice. First, it builds on preliminary results presented in [25] to construct a global MTF that can be applied to designs that include thin perfectly conducting sheets with arbitrary normal orientations that can cover completely or partially shared boundaries of neighboring domains, in addition to penetrable regions and perfectly conducting regions. This approach makes no distinction between sheets applied at coatings at material interfaces and free floating sheets internal to a domain. Second, a novel method is introduced to include PEC regions while suppressing internal resonances. In [15], a CFIE-type approach is presented to free the acoustic global MTF from internal resonances. The solution put forward here is based on a formulation that solves for both the electric and magnetic currents. In other words, the vanishing of the magnetic current on the boundary of perfect conductors is no longer an essential boundary condition (i.e., it is not a priori enforced by the choice of the spaces for trial and testing functions). Instead, it emerges as a natural condition (i.e., a property held by the solution of the variational formulation, by virtue of being the solution). The method is direct: its solution represents the tangential traces of the electric and magnetic fields. Third, this article introduces a Calderón preconditioner that is designed to be effective in the presence of combinations of penetrable domains, PEC regions, and thin sheets, in general stackings.
Thin sheets are considered to be embedded in one of the domains. In many cases, this is straightforward, but in the common case that the sheet appears as a coating at the interface between two regions, multiple choices can be made in assigning the sheet to a domain. Once a decision is made, generalized representation theorems are written down for these regions. This generalization can represent the tangential electric field in terms of the fields on the boundary of the domain, and the jumps of the magnetic field across the thin sheets embedded in that domain. A global MTF is then constructed by subtracting the representations in the exterior background region, and the representations in the interior domains, mimicking the construction of the singledomain PMCHWT.
The systems arising from boundary element discretizations of global MTFs can grow large, beyond the point where lower-upper (LU) factorization or other methods of cubic complexity cannot be used. The only viable alternative is the solution of the method by iterative methods of Krylov type. The complexity of these methods is pn 2 , where N is the number of DoFs used to describe the unknown field, and P is the number of iterations required for a Krylov solver to converge to a solution (up to an acceptable error). The number of iterations required by generalized minimal residual method (GMRES) depends on the field or values (a.k.a. the numerical range) of the system. Rigorous estimates can be found in [26], but unfortunately the bounds are not at all tight for wave problems, because of the nonnormality of the operators involved (see [27], [28]). Nevertheless, it is suggested in, e.g., [27], that a spectrum accumulating along a half axis in the complex plane and whose accumulation points are bounded away from both zero and infinity results in a small number of iterations. This situation can be achieved by multiplying both the leftand right-hand sides of the system by a suitable preconditioner. For conforming and stable discretizations of boundary integral equations with unknowns living on a manifold, this can be achieved by multiplying the system with the discretization of an integral operator of opposite derivative strength, discretized using a dual finite element space. This technique is referred to as Calderón preconditioning.
In this contribution, it is demonstrated that a Calderón preconditioner [29] can be constructed for the proposed global MTF. It is based on the general principles detailed in [30] and can be seen as an extension of the strategy proposed in [31]. We show that even in the presence of nontrivial geometric stackings including penetrable regions, perfectly conducting regions, and perfectly conducting sheets this approach leads to a significantly reduced iteration count and thus solution time. Following [9], only single-layer diagonal entries are retained to keep the computational cost of introducing the preconditioner minimal. This article is not focused on the low-frequency breakdown. The preconditioner demonstrated here will not stop the growth of the condition number as the frequency tends to zero. For a full solution Top: Ideas underlying the global multitrace method presented here will be explained by focusing on a geometry built from a penetrable slab (green-online), a perfectly conducting box (red-online, left), and a perfectly conducting sheet (red-online). These components are in fact touching but to stress the intuitive idea that gives rise to MTFs, they are shifted apart. Bottom: Schematic of the same situation. Relative orientation of the normals to geometrically coinciding surfaces affects how traces of double-layer potentials are computed.
to this problem, and the related problem of errors that creep into the solution and become manifest in secondary quantities such as the far-field, an approach based on quasi-Helmholtz splitting and rescaling is required [32].

II. GLOBAL MULTITRACE FORMULATION
Throughout this article, the frequency of all fields is assumed to be constant at some fixed value ω. To avoid having the main idea obfuscated by notation, we focus first on the geometry in Fig. 1. A penetrable domain 1 and a perfect electrical conductor domain 2 reside in an unbounded background domain 0 . The boundaries of 1 and 2 are denoted 1 and 2 , respectively. The normals n 1 and n 2 on 1 and 2 are chosen to point outward with respect to their defining domains. A thin sheet 3 represents a perfectly conducting coating on 2 but is conceptually treated as free floating in the background 0 . There is no standard choice for the normal field on 3 ; here, it is chosen to point toward 1 . The domain j is occupied by a material characterized by permittivity and permeability ( j , μ j ) or equivalently by wavenumber and impedance (κ j , η j ).
The single-and double-layer boundary integral operators T and K [1] are central to the construction of the global MTF with R = |x − y| and x ∈ p , y ∈ q . Although conceptually elegant, the implementation of global MTFs is delicate and requires modifications to the existing finite element and boundary element codes. The first point that requires attention is the correct computation of traces of double-layer potentials. Consider the situation in Fig. 1. When the observation point approaches the source region, the limiting value can be computed as the sum of the Cauchy principal value and a contribution ± j/2, where the sign depends on the direction the surface is approached from. Let Let σ i j be the relative orientation of the normals on i and j at a point in their intersection, and let α i j be +1 if the observation point x ∈ j is in front of i , and −1 if it is behind i . It holds that α 12 = −σ 12 α 21 . Then the tangential trace on 2 of the double-layer potential radiated by a source on 1 can be computed by reducing the general situation here to that of taking the standard trace of a double-layer operator (see [33]) For boundaries of regions, the normal is always chosen to point outward. As a consequence, whenever two boundaries meet, α i = +1. For thin sheets, there is no canonical way of choosing the orientation. In this scenario, the sign of the contribution j 1 /2 above needs to be carefully considered.
Another point to consider is that contributions like (1/2) j that appear in the interaction of boundaries and sheets that (partially) coincide should be considered geometric identities. This is best captured by the more verbose representation From this, it is clear that there are only nonzero contributions from the shared portions of 1 and 2 . In particular, for boundaries separated by a nonzero distance, the term disappears all together. Assembly routines encountered in most finite element codes for the treatment of local operators cannot be used here. A suitable algorithm: 1) iterates over the geometric elements of 2 ; 2) finds for each such element the geometrically coinciding element from 1 (if any exists); and 3) computes contributions and stores them in the correct matrix entries. One method to do this without increasing the asymptotic computational complexity is to store the geometric cells in an octree for fast lookup and collision detection. Where a conducting sheet is positioned exactly on the interface between two regions, a choice is required as to which domain the sheet belongs to. Here, we will always consider thin conducting sheets on the interface between dielectrics to be embedded in the background medium 0 . This is in contrast to the method presented in [34], where coatings are considered to be exactly at the boundary between regions. This results in a different treatment of free floating thin sheets and coatings and stands in the way of applying Calderón preconditioning. In [35], a local multitrace solver is introduced that allows the inclusion of perfectly conducting coatings. Also, this solver does not allow a Calderón preconditioner.
The starting point for building a global MTF is the representation theorem in the comprising regions. Structurally, the representation in terms of the traces is that of the same geometry, where the penetrable regions and the thin sheets are separated by a small layer of the background medium. The rigorous treatment of this, in case of a system comprising solely penetrable regions, relies upon extinction properties and can be found in [13]. The result of this careful analysis agrees with the intuition of taking the limit of the separation distance to zero. Equation (7), as shown at the bottom of the page.
A number of remarks is in place. The order and sign of the unknowns on the input and output sides of the block operator in (7) are chosen to bring out the coercivity of the operator. Because the number of iterations required by an iterative solver depends on the distribution of eigenvalues, which for discretized bilinear forms depends on the choice of basis, these choices do matter.
Even though these expressions are based on the representation formulas of Stratton-Chu, the operator in the left-hand side is not a Calderón projector. In particular, it does not map from a function space into itself; it uses the jump of the magnetic trace to compute the average of the electric trace. Nevertheless, these expressions allow for the construction of a global MTF that allows modeling of structures comprising penetrable media, perfectly conducting regions, and perfectly conducting thin sheets.
Under the assumptions that the interior region 1 contains neither sheets nor sources, the representation formula takes on the following form: Subtracting the exterior and interior representations gives the following condition (9), as shown at the bottom of the page, on the traces and jumps.
This system is not uniquely solvable. In particular, all the traces of homogeneous Maxwell solutions in the interior of the perfectly conducting domain 2 are in the nullspace. The treatment of such regions is the subject of Section III.

III. RESONANCE-FREE ELECTRIC-MAGNETIC FIELD INTEGRAL EQUATION
BEMs to treat scattering problems involving closed perfect conductors can be susceptible to internal resonances, i.e., unicity and existence of the solution can breakdown at a countable set of the so-called resonant frequencies. For the practitioner, these occurrences manifest as a sharp increase in the condition number of the system matrix, accompanied by a detrimental impact on the accuracy of the solution and the efficiency by which it can be calculated. In [6], a methodology for including both the dielectric and perfectly conducting regions is introduced that precludes the existence of internal resonances. The method is based on CFIE. It is not clear how a Calderón preconditioner can be built. Because of this, a novel methodology is developed in this contribution.
Our new method is inspired by the role that boundary layer operators at imaginary wavenumbers play in the construction of resonance-free CFIE-type integral equations [36]- [38]. In this section, a uniquely solvable integral equation for scattering by perfectly conducting regions is introduced. Given an incident EM field which is a solution of Maxwell's equations, the magnetic trace j on a perfect conductor is solution to both the MFIE and EFIE [33] The electric trace, on the other hand, vanishes: m = 0.
The system of equations for j can be extended with terms acting on m as long as the resulting system in (m, j ) allows a unique solution. After all, we know that the pair (0, j ) with j the shared solution of the EFIE and MFIE is a solution. That, in combination with uniqueness, makes it the solution. The freedom in choosing the terms involving m can be used to realize beneficial numerical properties. Consider the electric/magnetic field integral equation (EMFIE) The sign in front of the off-diagonal operators 1/2 is crucial. With opposing signs, the block operator would have been exactly the interior or exterior Calderón projector. These projectors have a very large nullspace and cannot be considered in searching for a system that can be uniquely solved for ( j, m). Unfortunately, numerical experiments show that this system remains susceptible to resonances. It is known that the single-layer operator at imaginary wave numbers does not allow resonances (see [39]) and thus provides a promising candidate for a modification of EMFIE that is free of resonances.
Replacing the top left operator with a Yukawa-type singlelayer operatorT retains the coercivity of the system and does lead to a unique solution at all frequencies withT In this treatment of perfectly conducting regions, the condition m = 0, tacitly assumed by most solution methods for the hard scattering problem appears as a natural condition, emerging upon solving the equations. The practical advantage is that when assigning electric and magnetic DoFs to the region boundaries, no distinction is necessary between penetrable and nonpenetrable objects and no manipulations as encountered in the construction of CFIEs are required.
Including the approach introduced above in the MTF amounts to changing the sign in front of one of the terms 1/2 and replacing the correct single-layer potential by its Yukawa equivalent. This yields (15), as shown at the bottom of the next page.
The manipulations carried through to arrive at this stage straightforwardly extend to other geometries. The three diagonal block operators in this example give a clear indication about how to treat the three types of region: penetrable domains, conducting domains, and conducting sheets.

IV. DISCRETIZATION AND PRECONDITIONING
The integral equation (15) is discretized in a straightforward manner. All the current densities are approximated as expansions in Rao-Wilton-Glisson (RWG) basis functions. No RWG functions are assigned to boundary edges of thin sheets, commensurate with the idea that the tangential magnetic field there is the same regardless the direction from which the boundary is approached (or equivalent, the idea that current cannot flow off the sheet). The equations are tested by the set of rotated (curl-conforming) RWG basis functions [40]. Note that the structure of the system does not depend on the adjacency information of the domains making up the study object. This is an advantage that has been recognized since the introduction of the first global multitrace integral equations. The resulting system is solved by feeding it to an iterative solver of Krylov type. In the numerical experiments, a GMRES solver has been used. Singular integrands can be computed using, e.g., the quadrature methods recorded in [41].
For meshes containing many edges, iterative methods become the only viable technique for the solution of (15) as the time for solution through the computation of an LU decomposition becomes prohibitively large. Unfortunately, the number of iterations needed by a Krylov solver grows as the mesh density parameter h tends to zero, resulting in an algorithm that has suboptimal asymptotic complexity.
This phenomenon can be linked to the spectral properties of the interaction matrix resulting from (15). It is affected to a large extent by the single-layer potential operators T on the diagonal-which are known to have a singular value distribution that stretches along the imaginary axis with one branch accumulating at zero and another at infinity. To mitigate this problem, a Calderón multiplicative preconditioner (CMP) that improves the eigenvalue distribution and consequently the condition number, and speed of convergence of the iterative solution is presented.
To demonstrate the flexibility of the approach introduced here, the preconditioning strategy is explained using a different geometry [see Fig. 2 (top)]: a thin PEC sheet is inserted between two penetrable regions. The MTF for this structure is given by (16), as shown at the bottom of the next page.
The system is built along the principles described in Sections II and III. The difference with the scenario treated there is that now there are two penetrable regions, and a single thin sheet. Another important difference are the signs in front of the contributions 1/2 stemming from interactions between the sheet and the penetrable regions. This is because the chosen normal to the sheet is in this configuration always pointing toward one of the penetrable regions and away from the other.
The boundary integral equation (16) is discretized and tested using rotated RWG functions, resulting in the linear system A Calderón preconditioner can be built starting from a dual discretization, using Buffa-Christiansen (BC) functions. The impact of the construction and use of the preconditioner can be minimized by only retaining the diagonal contributions. This is supported by both theoretical [30] and practical [14] considerations The extra factor "2" in the last diagonal element is included to compensate for the single occurrence of a single-layer operator on the system diagonal and further optimizes the final iteration count. Note that our choice does not attempt to mitigate ill-conditioning related to singular behavior near the boundary of open sheets, and its associated slow, logarithmic growth of the condition number. Strategies to deal with this are discussed in [42]. To map from the space of (rotated) RWG test coefficients to the space of BC expansion coefficients, multiplication by the inverse Gram matrix G −1 is required, with and m BC test functions and f (i) n RWG basis functions defined on the triangulation of i . The Calderón preconditioning of MTF (16) leads to the linear system V. NUMERICAL RESULTS

A. Resonance-Free Electric/Magnetic Field Integral Equation
To demonstrate the resonance-free approach introduced in Section III, Fig. 3 shows the plot of the condition number against frequency for a sphere with diameter d = 1.0 m. The structure is illuminated by a plane wave given by with wavelength ranging from λ = 0.62 m to λ = 3.14 m and κ = 2π/λ. The spikes in condition number at resonant frequencies are removed by use of a Yukawa-type single-layer operator.

B. Verification MTF Versus EFIE
As a first demonstration of the correctness of the method introduced here, the geometry from Fig. 1 is considered with the material in the penetrable region chosen equal to the material in the background medium (κ 0 , η 0 ) = (κ 1 , η 1 ) [see Fig. 2 (bottom)]. The resulting scattering problem is then solved using both the EFIE and our novel MTF.
The wavenumber of both the interior and exterior domains is set at 3 per meter. The impedance is set at η = 1 . The structure is illuminated by the plane wave (21). The results for the amplitude of the induced electric current in Fig. 4 show good agreement. Fig. 2. Top: Composite structure featuring two dielectric layers (blue-online) and a thin PEC sheet (yellow-online) covering the interface between the layers and part of the exterior boundary. These components are in fact touching but to stress the intuitive idea that gives rise to MTFs, they are shifted apart in this visual representation. Bottom: Perfectly conducting object comprising two connected components: a perfectly conducting domain and a perfectly conducting sheet. In a second verification, the same structure is treated, but now the penetrable region contrasts against the background: κ 1 = 2.4κ 0 , η 1 = η 0 . The problem is solved using both the classic STF and our novel MTF. Fig. 5 shows good agreement between the two approaches.
In Fig. 6 (left), the relative rms error in the far-field is compared against STF solution on the finest mesh. Good convergence is observed.

C. Impact of Preconditioning on the Iteration Count
To demonstrate the efficacy of the preconditioner proposed above, the number of iterations required for the GMRES solution to converge within 10 −8 is recorded in two situations.
The first example uses the geometry from Fig. 1 with κ = 3 per meter, κ 1 = 2.4κ 0 , and η 1 = η 0 = 1 . The structure is illuminated by the plane wave (21). The number of DoFs in the surface fields is 11 524. The number of iterations required for the GMRES solver to converge without preconditioner is 1368, and with the block diagonal Calderón preconditioner it is reduced to 168, corresponding to a speedup in solution time of about a factor 4.
For the second example, we reconsider the geometry where two dielectric boxes are juxtaposed [see Fig. 2 (top)]. One of the boxes is partially wrapped by a PEC sheet. The PEC covers both parts of the boundary of that box that are exposed to the background material and the parts of the boundary that are   6. Left: Far-field error (two-norm of the difference between the far-field computed from MTF and STF) calculated for the structure in Fig. 1 (top) converges to zero as the mesh density (number of edges per meter) is increased. Right: Number of iterations versus mesh discretization (number of edges per meter) for the vanilla and preconditioned system arising from a thin PEC sheet wrapped around one of two adjacent dielectric slabs. δ = 1/ h. shared with the other box. Note that this requires no special geometric preprocessing or modifications to the formulation compared with a situation where coating is entirely on the shared interface or entirely on the exterior boundary. The dielectric box on the left 1 has material properties κ 1 = 2.4κ 0 , and η 1 = η 0 , and the dielectric 2 box on the right κ 2 = 1.4κ 0 , and η 2 = η 0 . A thin PEC screen 3 covers part of the boundary of 2 . The background medium 0 is unbounded and is defined by material parameters κ 0 = 3.0 and η 0 = 1. The surface of the structure is approximated by a triangular mesh with edge size h. When the plane wave (21) traveling in theẑ-direction illuminates the structure, the total electric and magnetic fields can be computed from the solution of the system (16). Fig. 6 (right) shows the impact of Calderón preconditioning: although extra work is performed to assemble the preconditioner, the number of iterations required for convergence appears almost constant (at ≈100 iterations) with increasing mesh density. On the other hand, the iteration count for the unconditioned system grows at a significantly higher rate, thus taking an increasing number of iterations to reach a solution as mesh size decreases. Fig. 7 shows the z-component of the electric field. Comparison with the field plot obtained from the STF shows good agreement.

VI. CONCLUSION
This article has introduced: 1) a global MTF for the inclusion of thin perfectly conducting (PEC) sheets and coatings in composite structures; 2) a Calderón multiplicative preconditioner for the global MTF (MTF); and 3) a novel resonance-free method based on the standard EMFIEs to model scattering by perfect conductors. A range of examples featuring various combinations of thin PEC sheets, penetrable materials, and bulk PECs were demonstrated. Calderón preconditioning has been shown to become increasingly effective in speeding up the solution as the mesh size tends to zero. The resonance-free method developed was applied to the canonical example of scattering by a PEC sphere across a range of frequencies, without suffering from resonances. The methods introduced in this article enable the modeling of EM scattering by composite structures, including thin sheets and coatings, in an elegant, scalable (with mesh density), and resonancefree way.