Accuracy Analysis of Div-Conforming Hierarchical Higher-Order Discretization Schemes for the Magnetic Field Integral Equation

The magnetic field surface integral equation for perfect electrically conducting scatterers suffers from accuracy problems when discretized with lowest-order Rao-Wilton-Glisson (RWG) functions. For high-frequency scattering scenarios, one of the various reported countermeasures are hierarchical higher-order (HO) functions. We demonstrate that the accuracy of these HO methods of up to 1.5th order may be further improved by employing a weak-form discretization scheme for the identity operator inside the magnetic field integral equation (MFIE), in particular for scatterers with sharp edges. As expected, the presented numerical results indicate that this approach becomes less effective for increasing order. Moreover, since the weak-form discretization overcomes only the anisotropy problems of the standard discretizations, parts of the accuracy problems of the MFIE persist for HO discretizations if the testing is performed with non dual-space conforming functions.


I. INTRODUCTION
E LECTROMAGNETIC radiation and scattering problems for perfect electrically conducting (PEC) objects are often solved by the method of moments (MoM) in conjunction with integral equations (IEs) [1]. Boundary IEs are especially popular due to their reduced complexity with surface discretization only as compared to volumetric formulations. The electric field IE (EFIE) which works with the tangential electric fields on the surface of a scatterer is suitable for open and closed structures, whereas the magnetic field IE (MFIE) for the tangential magnetic field on the surface is suitable for closed bodies only. When employed for closed bodies, the EFIE and MFIE are both only solvable with some kind of source constraint. The most common restriction is to choose specific equivalent sources, the so-colled Love-currents, and a superposition of EFIE and MFIE, the combined field IE (CFIE) [2]. However, both the EFIE and MFIE have their respective problems, and those problems also affect the CFIE, of course. In the classical discretization with div-conforming Rao-Wilton-Glisson (RWG) functions [3] or their higher-order (HO) counterparts, the issues of these two IEs are quite contrary to each other. While the EFIE is generally seen as the gold-standard in terms of accuracy, the conditioning of the matrix suffers at low frequencies and for dense discretizations. The MFIE is commonly praised for its excellent conditioning (due to the involved well-conditioned Gram matrix stemming from the discretization of the identity operator in the integral equation of the second kind), however, it suffers from inaccuracies which appear in particular in a pronounced manner for sharp-edged objects and low-frequency scenarios [4], [5], [6], [7]. These accuracy issues are still subject of ongoing research. The established explanation is that the evaluation of the magnetic field is not performed in the dual space of the magnetic field, when the same div-conforming basis and testing functions are employed. Unfortunately, this kind of testing is required for a well-conditioned Gram matrix when working with RWG functions only. In contrast, dual-space testing is easily achieved for the EFIE by using curl-conforming testing functions, which is one reason for the superior accuracy of the EFIE.
The first fully satisfying solution to the MFIE accuracy problems originated from the idea that the testing functions may not only be dual in the sense of curl-conformity but also defined on the dual (i.e., barycentrically refined) mesh [8], [9], [10], [11], [12], [13]. These so-called Buffa-Christiansen (BC) functions are, however, only the counterparts for the lowest-order (LO) RWG functions and entail increased computational expenditure due to the mesh handling. HO extensions to the BC functions exist, but have never been employed for the PEC MFIE [14], [15].
There exist numerous approaches to fix the MFIE accuracy issues only for a subset of its problematic scenarios. For lowfrequency scenarios, the razor-blade testing functions are worth mentioning [16]. For high-frequency scenarios, there are probably too many to list them all. Most of them may be categorized into approaches with div-conforming (mostly RWG) functions and into ones with curl-conforming functions. The latter include monopolar RWGs as well as rotated curl-conforming RWGs, which are inherently inappropriate to represent div-conforming surface current densities [17], [18], [19], [20]. This causes problems when combined with the EFIE to the CFIE since the occurring hypersingular line charges are at least rather difficult to handle, if not impossible [21]. One recent approach comprises interior volumetric testing functions in order to mitigate the implications [22] -nevertheless, the surface fields diverge with such a basis.
The approaches with div-conforming basis functions often choose curl-conforming testing functions or a mixture of curland div-conforming ones [23], [24], [25], [26]. Furthermore, it was demonstrated that div-conforming HO expansion approaches may improve the accuracy of both the EFIE and MFIE [27], [28], [29], [30], [31], [32], [33], [34], where it remains unclear whether the MFIE inaccuracy issues are really completely resolved or only partially. More recently, it was proposed to improve the accuracy of the high-frequency RWG MFIE by applying a weak-form (WF) discretization scheme for the discretized identity operator while working with classical div-conforming RWG testing functions [26]. Compared to the theoretically more sound approach with BC testing functions, the computational cost of this scheme is much lower since no barycentric refinement of the mesh is needed, and the additional effort related to the RWG Gram matrix is negligible.
In this article, the WF identity operator discretization scheme is introduced and investigated together with HO expansion functions up to order 1.5. In particalur, we demonstrate that a hierarchical HO basis [28], [29], [30] is only capable of curing the MFIE inaccuracies in part and additional measures seem to be required to boost the MFIE accuracy to really good levels. Preliminary results regarding the accurate full first-order discretization of the MFIE only are found in [34].
This article is structured as follows. Section II discusses the classical discretization of the MFIE with a div-conforming (HO) basis. In Section III, the WF hierarchical HO discretization scheme for the MFIE is introduced and analyzed. Numerical results in Section IV demonstrate the improved accuracy.

A. Electromagnetic Scattering From a PEC Object
The MFIE is related to the Love current condition on the surface s of a PEC scatterer, with the electric surface current density j, the outward unit normal vector n on s, and the total magnetic field h. Decomposing the magnetic field into incident (or excitation) and scattered parts as h = h inc + h sca , and replacing the scattered field by radiation operators acting on the current density, yields the MFIE including the identity operator where the integral is evaluated in the Cauchy principal value sense. Here, r denotes the observation coordinate, r the source coordinate, δ(r − r ) the Dirac delta distribution, and g(r, r ) the free-space Green's function which is employed for an equivalent free-space scenario with zero field inside the scattering object. This is possible due to the Love current description and the application of the equivalence theorem, see Fig. 1.

B. Classical RWG Discretization
We consider a div-conforming hierarchical HO basis β defined on a triangular mesh -the RWG functions constitute the LO subset of these functions and are referred to as β LO , while the complementary hierarchical HO part is referred to as β HO . The rotated versions of these functions are denoted by α = β × n.
The fields produced by the electric surface current densities, which are expanded by the aforementioned HO functions as are tested in a standard MoM scheme with the very same functions. In this way, the linear system of equations of the classical is obtained. The unknown coefficients of the basis functions are collected in the vector i. Both the right-hand side and the matrix entries are calculated similar to inner products as Subsets of the div-conforming hierarchical HO functions complete up to pth polynomial order of the expansion functions are described as follows [28], [29], [30]. The LO part, also called first-order "quasi"-gradient space, contains just the RWG functions and exhibits p = 0.5 order. The full-first order functions with p = 1 supplement the first-order "quasi"-gradient space by a first-order rotational space. Furthermore, we consider the case p = 1.5, which includes the second-order "quasi"-gradient space in addition [28], [29], [30]. Examples of the surface current distributions on (pairs of, where appropriate) triangular facets are depicted in Fig. 2.

C. Comparison Solutions to the PEC Scattering Problem
Utilizing the rotated curl-conforming and quasi-divconforming BC functionsα defined on the barycentrically refined mesh, we employ the BC-tested MFIE (10) with RWG ansatz functions as a comparison solution for an accurate discretization of the LO MFIE, i.e., p = 0.5. The Julia implementation of the boundary element analysis and simulation toolkit (BEAST) has been utilized to generate all results for the BC-tested MFIE shown during the course of this article [35].
Furthermore, the classical EFIE tested with rotated RWGs α is employed as an accurate comparison method, and as for obtaining a reference solution either on a finer mesh or with ansatz functions with p > 1.5.

III. WEAK-FORM IDENTITY OPERATOR DISCRETIZATION FOR THE HIGHER-ORDER MFIE
The following investigations are tailored to the used hierarchical HO functions. Other HO approaches, especially interpolatory ones, may not directly contain an RWG-alike subset of the functions. However, identifying and isolating the anisotropic influence of the RWG functions is crucially important to improve the accuracy of the high-frequency MFIE as considered in this work. By the RWG anisotropy, we mean the anisotropic spatial filtering of the identity operator due to the geometric restrictions of the RWG discretization, that is, low-order basis functions defined on triangles of different shapes and sizes.

A. The Lowest-Order Case
The foundation of the proposed formulation is a WF identity operator discretization scheme, which was originally proposed for RWG discretization. Instead of taking the pure Gram matrix G β LO ,β LO in (7), an added basis transformation matrix [26] (11) is introduced, where G α LO ,β LO is the Gram matrix of α LO and β LO functions and I is the identity matrix of size equal to the number of RWGs. The weighting factor γ determines how much the original, classical-MFIE identity operator discretization, i.e., the RWG Gram matrix, is present in W LO γ . With γ = 1, the pure classical MFIE is attained, whereas lower values increase the weighting of the WF part, which consists of a strong form rotation by 180 • (i.e., a minus sign) and two subsequent WF rotations by 90 • (each realized by the matrix multiplication . Too small weightings with γ → 0 should not be chosen since the matrix G α LO ,β LO suffers from a non-trivial kernel, which might negatively affect the condition number of the resulting system matrix.
The mentioned previous work has considered RWGs only. The therein improved discretization scheme for the MFIE reads [26] An appropriate weighting factor for this case p = 0.5 is γ 0.5 = 0.5, where the basis function order p = 0.5 was included as a subscript. Furthermore, the weighting may be adjusted to be slightly weaker for basis functions located near edges (of the object, not the mesh, of course) [36].

B. The Higher-Order Case
Since we assume that the anisotropy of the RWG functions is the major cause of the problems of the identity operator discretization, we tackle only the RWG part of the Gram matrix. 1 For doing so, the whole Gram matrix G β,β is split into LO and HO interaction blocks according to The respective LO expansion coefficients are represented by the vector i LO , the HO part by i HO . Then, the WF basis transformation scheme including two WF rotations is employed only for the LO subset of the basis functions as described in (11). This alters the Gram matrix into where the weighting factor γ p depends on the polynomial order p of the expansion functions but is fixed for the whole matrix otherwise. Note that the basis transformation (11) is in fact not even employed for the LO functions as a whole set of basis functions but just for the self-interaction Gram matrix block, i.e., for the inner products of LO with LO functions themselves. In consequence, the hierarchical HO MFIE including a WF identity operator representation follows as In the hierarchical HO discretization, a part of the RWG anisotropy is already taken care of by the expansion functions themselves; this effect should increase with increasing order of the expansion functions. Hence, we expect that γ p has to come closer to one, the higher the order p of the expansion functions is chosen. This means that the discretization comes closer to the classical one overall. In the following, we investigate two sets of HO functions with p = 1 and p = 1.5 and the respective choices for γ p .

IV. NUMERICAL RESULTS
In the following, we analyze four PEC scattering scenarios: two cubes with 0.5λ and 1λ side length, a sphere with a diameter of 9λ, and the stealth object Flamme. We analyze the matrix condition numbers, the iterative solver convergence and the RCS FF error of the solutions obtained with RWG basis functions and HO functions up to order 1.5. Other aspects, such as the accuracy improvements of the WF-MFIE in the NF (not just in the FF), or the convergence properties of the charges and currents have already been studied for the RWG case in [26]. We do not expect substantial changes with respect to the previous findings. In particular, the convergence properties and low-frequency convergence issues of charge and current of the WF scheme are inherited from the RWG MFIE.

A. Choosing a Suitable Weighting Factor for Accurate MFIE Scattering From a Cube
The weighting factor γ p exhibits some dependence on the considered scenario. For instance, the optimal choice was around 0.4 for a sharp-edged pyramid and around 0.55 for a smooth sphere in [26]. Here, we carry out an analysis over a set of scattering scenarios in order to gain sufficient insight toward the choice of the parameter in a problem-adapted manner.
For the first numerical investigation, we consider a 0.5λ square cuboid with 14 different triangular meshes ranging from 0.28λ to 0.035λ average edge lengths h. The results are summarized in Fig. 3. The number of unknowns ranges from 72 (coarsest mesh, RWGs) to 12 420 (finest mesh, p = 1.5). The radar cross section (RCS) of the PEC cube scattering for plane-wave incidence is evaluated in the far-field (FF) for the MFIE, the EFIE, the BC-tested MFIE, as well as the proposed WF-MFIE. The BC-tested MFIE is only considered for the case with RWG basis functions, i.e., p = 0.5. Then, the relative (arithmetic) average FF error with respect to a reference solution (a 2.5th order EFIE solution on the finest among the considered meshes) is computed for each choice of the mesh and also for each choice of γ p in the WF-MFIE. For all these mesh variations, the geometric average 2 of the individual arithmetically averaged FF errors is then calculated. Besides the RCS FF error levels of the considered IEs, the (min/max) spread of the WF-MFIE error as compared to the EFIE error for each individual simulation is indicated by error bars in Fig. 3. Obviously, the case γ = 1 of the WF-MFIE covers the case of the classical MFIE.
For attaining these RCS FFs, the linear systems of equations are solved with GMRES down to a residual of 10 −6 ; only some of the WF-MFIE versions with low γ values are limited to a maximum of 1000 iterations in addition since they do not converge at all due to the non-trivial kernel of G α LO ,β LO , which results in large residuals and, subsequently, errors. The Gram matrix inversion inside the matrix-vector-product of the WF-MFIE are carried out on the fly with a diagonally-preconditioned conjugate-gradient solver with a residual threshold of 10 −10 -this threshold should be chosen at least one or two orders of magnitude lower than the solver threshold at which the solution error (e.g., evaluated in the NF or FF) stagnates.
As expected, the best value of γ increases with increasing p. For the considered scenario, we find γ 0.5 ≈ 0.4 for p = 0.5, γ 1.0 ≈ 0.7 for p = 1.0, and γ 1.5 ≈ 0.9 for p = 1.5 to be optimal. However, as already mentioned, we know that the best weighting factor for the case p = 0.5 is scenario-dependent in the range γ 0.5 ∈ [0.4, 0.55]. In order to be on the safe side regarding the possibly negative influence of ker G α LO   mesh refinement analyses for these proposed values of γ p , p ∈ {0.5, 1.0}, have been presented in [34].
Going into a more in-depth analysis of Fig. 3, diminishing gains for HO expansions are observed. Compared to the LO discretization, the gains are observed only on a reduced level and the required value of γ 1.5 comes close to 1, i.e., the choice for the pure classical MFIE. In Fig. 3, the maximum accuracy gains as compared to the classical MFIE are −10 dB for p = 0.5, −4.6 dB for p = 1.0, and −2.3 dB for p = 1.5. For other scattering scenarios, the accuracy improvements for p = 1.0 and p = 1.5 are also in similar regions as for this 0.5λ cube (sometimes also better). Further, the error spread also increases with increasing p; both for the MFIE and the WF-MFIE. The meaning of this is best investigated by showing the underlying errors in more detail, e.g., for the best choices of γ.
Another interesting observation is that the RWG WF-MFIE, the BC-MFIE, and the classical full first-order MFIE are similarly accurate, and all three are better than the classical RWG MFIE, see Fig. 3(a) and (b). The HO MFIE is only able to achieve the same accuracy as accurate testing schemes with the RWG basis for a higher discretization density (i.e., on the same mesh but with HO functions). In Fig. 4, the number of electric current unknowns N instead of the edge length is plotted on the abscissa. The classical full first-order MFIE uses the additionally introduced unknowns inefficiently; it is barely more accurate than the RWG MFIE for the same number of unknowns. The 1.5th order HO MFIE achieves the same accuracy as the RWG WF-MFIE and the BC-MFIE 3 at the same number of unknowns. The 1.5th order WF-MFIE still has an accuracy advantage, though. In the past, such observations have led to the conclusion that the HO MFIE fixes the RWG MFIE accuracy problems [29], [31], [32]. Nonetheless, we have already demonstrated that further improvements for the accuracy of the HO MFIE are 3 Keep in mind that the compuational effort for the BC-MFIE is higher than for other discretizations with RWG basis for the same number of unknowns in Fig. 4. indeed possible, and also how to achieve these improvements for high-frequency scenarios. For order two and above, we do, however, expect that the weighting factor γ of the proposed WF scheme converges to one.

B. Detailed Accuracy Analysis and Scattering Results for a Cube
A similar mesh refinement analysis -this time just for the chosen values of γ p -is repeated for a 1λ cube with 20 meshes ranging from 0.37λ to 0.027λ edge lengths in Fig. 5 (with the number of unknowns between 162 and 100 800).
The accuracy improvements are similar to the first considered scenario. Again, the benefits of the WF discretization diminish with increasing p. Just as in the first example, there are still significant improvements for the full-first order discretizations, partially coming close to the EFIE accuracy; while these are not really observed anymore to the same extent for the 1.5th order solution. Another possibly interesting fact is that a full-first order EFIE exhibits decreased accuracy for the 0.5λ cube but unaffected accuracy for the 1λ cube, i.e., this is a frequencydependent effect.
Finally, we show some more of the underlying data for the respective meshes with h ≈ λ/10 of the presented refinement  studies. The considered square cuboid with 1λ edge length is aligned with the Cartesian coordinate axes. In Fig. 6, the copolarized bi-static RCS caused by a plane wave incidence with k = −e z and e x polarization is shown for the MFIE and the reference solution in the ϕ = 0 • cut. Since these two RCS curves are almost indistinguishable, the further IEs are only included for their respective relative FF RCS error in the considered ϑpolarization The error levels for the mesh considered in Fig. 6 are summarized in Table I for all considered IEs and orders, regarding the maximum error max , the arithmetic average avg , and the error of the monostatic RCS (0 • , 0 • ). The error of the WF-MFIE is almost 10 dB lower for the LO case and about 5 dB lower for the HO cases. While the full first-order MFIE is already more accurate than the BC-MFIE, the WF-MFIE shows additional accuracy gains. This observation demonstrates that the HO MFIE with classical div-conforming discretization is not able to solve the MFIE accuracy problem fully since further accuracy improvements are feasible. For the mesh refinement study, whose accuracy was analyzed in Fig. 5, we show the iterations of a GMRES solver for convergence to a threshold of 10 −6 in Fig. 7. This gives insights into the conditioning behavior of the different IEs. The EFIE shows the typical dense-mesh breakdown with an increasing number of iterations, while the three MFIEs show a stable behavior for all mesh densities. The BC-MFIE shows slightly better conditioning than the RWG MFIE, while the WF-MFIE (for all three sets of expansion functions) shows basically the same convergence behavior as the standard MFIE. This means that the WF scheme does not change the good conditioning of the MFIE.

C. Investigation of a Sphere-HO Benefits and WF Identity Limitations
The so-far analyzed objects were cubes, i.e., with edges. It is known from literature that the accuracy problems of the RWG MFIE are particularly severe for such objects. For the RWG MFIE, accuracy issues are also observed for smoother surfaces. Now, we consider a sphere with three different meshes with roughly the same number of unknowns. The finest mesh with λ/10 average triangle edge length comprises 88224 RWG unknowns, the coarsest mesh is solved for 88080 1.5th order unknowns, and the mesh with a discretization density in between the aforementioned ones exhibits 88056 first-order unknowns. Since the surfaces of the three meshes deviate slightly (leading to different solutions), the accuracies are judged by an individually calculated 1.5th order reference solution on three respective  Table II; and also at γ p values closer to one than before. Hence, we conclude that the proposed scheme is mostly suited for sharp-edged objects. Additionally, the limitations of the HO MFIE also become clear again. The full-first order CFIE is just as accurate as the RWG WF-CFIE, i.e., doubling the number of unknowns and degrees of freedom just increased the accuracy to levels which are also achieveable with an accurate LO discretization scheme. Hence, we conclude that the div-conforming hierarchical HO functions do not solve the accuracy issues of the (HO) MFIE in this scenario.

D. A Larger Scattering Scenario
As a scattering scenario closer to real-world applications, we consider scattering from the PEC stealth object Flamme for a plane wave incident from z direction and with e x polarization [37], [38]. First, we consider one mesh with 52782 triangular facets, i.e., 79173, 158346, and 263910 electric current unknowns for p = {0.5, 1, 1.5}. The average edge length of the triangles is λ/10. The reference solution is calculated with a 0.9-CFIE on a refined mesh exhibiting 1055640 unknowns. The error levels are summarized in Table III. Overall, accuracy improvements in the range of 5 dB for the WF-CFIE are observed. It has to be noted that -unlike as for the LO MFIE -the choice of the optimal γ p for γ ≥ 1.0 seems to be a bit more delicate. In Table III, the previous choices have been employed. For the Flamme, the optimal γ 1.0 is a bit larger and the optimal γ 1.5 is a bit lower, though.
In contrast to the previous scenarios, the full-first order CFIE is more accurate as than the accuracy improved WF-CFIE. However, the RWG WF-CFIE improvements are not as pronounced in this case. As before, the p = 1.0 CFIE just achieves about RWG EFIE accuracy levels, which might be also possible with an accurate dual-space-tested MFIE. Furthermore, the step towards p = 1.5 does not give the same improvements for the CFIE as for the EFIE. Overall, we conclude also in this case that the MFIE inaccuracy persists in the HO case, when the testing is done with div-conforming functions.

V. CONCLUSION
The classical testing for the LO MFIE with div-conforming functions is well known to be prone to accuracy issues. We have demonstrated that these issues persist with hierarchical HO discretizations of order up to 1.5 -in contrast to first-sight impressions based on the observation that the HO MFIE, in particular already the full first-order MFIE, achieves similar accuracy levels as accurate MFIE discretizations schemes with RWG basis (BC-MFIE or WF-MFIE). To relieve these accuracy issues, a WF discretization scheme was proposed, which takes care of the anisotropy of the RWG part of the HO functions and is able to mitigate the accuracy issues in part -this was demonstrated in particular for sharp-edged objects. The resulting HO discretization of order 1 and 1.5 can be more accurate than the classical counterparts, which only achieve the same level as the BC-MFIE or the WF-MFIE with RWG basis functions. Still, the accuracy issues of the HO MFIE due to the lack of dual-space testing remain of course even with our WF testing scheme and the achievable accuracy does not reach the level of the HO EFIE. As observed in this article, the effect of the WF scheme on the identity operator discretization is quite strong for the lowest order (0.5), since the Dirac-delta kernel leads to the well known RWG Gram matrices with only few interactions between neighboring and nearby triangles. The WF scheme introduces a smoothing effect-similar to the enlarged testing area of BC functions. With increasing basis function order, there are more interactions among the different basis functions. This has a similar smoothing effect as enlarging the testing area or the WF scheme and renders the WF scheme thus superfluous for increasing order. Hence, a reduced effect of the WF scheme is observed for full first-order and even more so for 1.5th order basis functions.