Dyadic multipole-based generalized source integral equations

The approach of multipole-based generalized source integral equation (GSIE) formulations for the scattering by essentially-convex impenetrable objects is extended to dyadic problems. This is demonstrated through the two-dimensional (2-D) problem of transverse-electric (TE) scattering. For this case, the principles of multipole-based generalized source design are utilized within the framework of dyadic TE-GSIEs introduced recently for reflective shield sources. The derivation of the dyadic auxiliary contribution to the modified Green's function is presented in detail and is shown to provide similar superior low-rank compressibility of off-diagonal method-of-moments' matrix blocks as its transverse-magnetic counterpart, while maintaining error controllability. The extension enables the treatment of impedance boundary scatterers in 2-D and is a crucial stepping stone toward a full three-dimensional vector formulation that is expected to exhibit enhanced low-rank compressibility for a broad range of scatterer geometries.

1 Introduction

The scattering of an electromagnetic wave by an impenetrable object can be modeled using a surface integral equation (IE). The discretization of the IE, for example by using the method of moments (MoM) (Harrington, 1993), leads to dense matrices. To enable the treatment of electrically large problems with many unknowns, fast solvers have been proposed (Song and Chew, 1995; Bleszynski et al., 1996; Phillips and White, 1997; Boag et al., 2002; Brick and Boag, 2010; Yang and Yılmaz, 2012; Wei and Yılmaz, 2014). Among these, kernel-independent solvers, e.g., Hackbusch (1999), Adams et al. (2005), Zhao et al. (2005), Tamayo et al. (2011), and Liu et al. (2021), and, in particular, fast direct solvers (Shaeffer, 2008; Heldring et al., 2011; Brick and Boag, 2011; Chai and Jiao, 2013; Guo et al., 2017) rely on the algebraic compression of certain MoM matrix blocks. Conventionally, the off-diagonal blocks are represented by their low-rank (LR) approximations, where the desired compression error threshold τ determines the block ranks ℛ(τ).

The performance of fast solvers can be determined by the compressibility of these blocks. A reduction of the asymptotic scaling of the complexity below that of a straightforward solution requires that the ranks of the largest blocks scale slower with the electrical length than the number of problem unknowns. For conventional IE kernels, this is only guaranteed for a reduced dimensionality interaction, e.g., between parts of elongated or quasiplanar scatterers (Michielssen et al., 1996; Martinsson and Rokhlin, 2007; Corona et al., 2015; Brick, 2021). For arbitrary surface geometries, this is usually not the case. There, for blocks that represent interactions that include broadside components, the ranks scale asymptotically linearly with the number of unknowns, resulting in a computational complexity reduction by merely a constant factor (Brick and Yılmaz, 2016; Brick et al., 2026). While more sophisticated structures (Michielssen and Boag, 1996; Guo et al., 2017; Kaplan and Brick, 2022; Sayed et al., 2022) enable a slower asymptotic scaling of the costs, they often exhibit late cross-over with LR-compression-based solvers (Brick, 2021) – particularly for direct solvers where specialized arithmetic (see Liu et al., 2021) is required. The need for slowing the scaling of the rank, which can serve both conventional LR and later generations of algebraic compression fast solvers, has led to the development of generalized IE formulations (Boag and Lomakin, 2012; Brick et al., 2014; Klinkenbusch et al., 2018; Sharshevsky et al., 2020; Zvulun et al., 2023; Dahan and Brick, 2024; Dahan et al., 2026; Kalhöfer et al., 2025).

For essentially-convex scatterers (see Fig. 1a), the unfavorable scaling of the rank in conventional IEs, which are derived by using Love's equivalence principle, stems from the line-of-sight (LoS) broadside interactions between opposite side scatterer subdomains (Boag and Lomakin, 2012) (see Fig. 1b–c). Fortunately, for such scatterers, the IE kernel can be modified such that it includes, for each source on the surface, an additional contribution that appears to emanate from within the scatterer. If designed properly, these contributions can approximately cancel the broadside component of the interactions (see Fig. 1d), leading to the reduction of their effective dimensionality and, therefore, of the rate in which the ranks scale with the frequency.

Figure 1 (a) TE scattering by an essentially-convex PEC surface S with inner convex hull S_c, (b) line-of-sight interactions of equivalent surface currents in conventional integral equations, and (c) broadside interaction between subdomains S^s and S^o. (d) GSIE MGD exhibits deep shadow in broad regions of the volume defined by S.

Following this idea, various types of generalized source integral equations (GSIEs) have been proposed (Brick et al., 2014; Klinkenbusch et al., 2018; Sharshevsky et al., 2020; Zvulun et al., 2023). Particularly, it was shown that additive auxiliary kernels based on dependent multipole contributions can provide a deep shadow and a related attenuation of the broadside interactions for very low thresholds (Kalhöfer et al., 2025). Unlike earlier GSIE designs, multipole-based kernels do not involve expensive field integrals over shielding domains. Thanks to the singular source point of a multipole expansion, even surface points in narrow boundary regions can be augmented with such auxiliary contributions, making the approach applicable to a broader range of geometries. As the expansion requires only little knowledge of the scatterer geometry, the multipole-based kernel can be easily integrated into kernel-independent solvers and arbitrarily seamlessly modified if needed. This concept was introduced and reported in detail in Kalhöfer et al. (2025), through the two-dimensional (2-D) scalar problem of transverse-magnetic (TM) scattering by perfect electric conducting (PEC) objects, which focused on the implications to the design of fast direct solvers and the study of their performance. Further generalization of the formulation, toward enabling the treatment of large impenetrable three-dimensional (3-D) objects, requires its extension to vector- and dyadic-kernel formulations.

In this work, a dyadic extension of the multipole-based GSIE in Kalhöfer et al. (2025) is presented. This is done for the representative example of transverse-electric (TE) scattering problems. Following the vectorial approach in Dahan et al. (2026), two vector contributions (i.e., two sets of multipole amplitudes) are computed, one for each of the orthogonal surface-current components, resulting in a modified dyadic Green's function. Each of the two auxiliary kernels is produced using magnetic multipoles, i.e., scalar entities that produce a z-directed magnetic field contribution and a (with respect to the z-axis) TE field contribution. This representation avoids duplication of the number of design variables and enables, via duality, a convenient repurposing of routines from Kalhöfer et al. (2025). This process is presented in detail, the resulting IE formulation and its moment solution are validated, and the influence of the construction parameters on the compressibility of the MoM matrix and on the accuracy of the related solution is studied.

The remainder of this manuscript is structured as follows: Sect. 2 formulates the problem of designing a dyadic GSIE kernel. Section 3 describes the design of such a multipole-based kernel for TE problems. Section 4 analyzes the properties of the proposed kernel and demonstrates its effectiveness through numerical examples. Section 5 concludes the work.

2 Problem formulation

Consider the problem of an incident wave, described by electric and magnetic fields E^inc and H^inc, respectively (time-harmonic dependence e^jωt is assumed and suppressed) illuminating a highly conducting structure modeled as an impenetrable surface S with unit normal vector $\widehat{\mathit{n}}$ on which an impedance boundary condition (IBC) is imposed (see Fig. 1a). Corresponding scattered fields, E^sca and H^sca arise such that the total fields satisfy the Leontovich boundary condition $\widehat{\mathit{n}}\times \mathit{E}={\mathit{\eta }}_{\mathrm{s}}\widehat{\mathit{n}}\times (\widehat{\mathit{n}}\times \mathit{H})$ (Senior and Volakis, 1995) with the surface impedance η_s on S⁺ (the exterior side of S). When η_s=0, this reduces to the conventional electric field integral equation (EFIE) for a PEC scatterer.

Using Love's equivalence principle (Harrington, 1993), the scattered fields are conventionally represented as produced by a surface current J radiating in free space (see Fig. 1b) such that

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{E}}^{\mathrm{sca}}\left(\mathit{r}\right)& {\displaystyle }={\int }_S{\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\cdot \mathit{J}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }\text{(2)}& {\displaystyle }{\mathit{H}}^{\mathrm{sca}}\left(\mathit{r}\right)& {\displaystyle }={\displaystyle \frac{j}{k\mathit{\eta }}}{\int }_S\mathbf{\nabla }\times {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\cdot \mathit{J}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }\end{array}$$

outside S, where η is the wave impedance and ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })$ is the dyadic Green's function of the free space

$$\begin{array}{}\text{(3)}& {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })=\left(\mathbf{I}+{\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}}}\mathbf{\nabla }\mathbf{\nabla }\right)g^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\end{array}$$

with the wave number $k=\mathrm{2}\mathit{\pi }/\mathit{\lambda }$, the identity dyadic I, and the scalar Green’s function $g^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })=-jk\mathit{\eta }\exp (-jk|\mathit{r}-{\mathit{r}}^{\prime }\left|\right)/\mathrm{4}\mathit{\pi }|\mathit{r}-{\mathit{r}}^{\prime }|$. From the boundary condition on the tangential electric field follows the IBC EFIE

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \widehat{\mathit{n}}\left(\mathit{r}\right)\times {\int }_S{\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\cdot \mathit{J}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }\\ & -{\displaystyle \frac{j{\mathit{\eta }}_{\mathrm{s}}}{k\mathit{\eta }}}\widehat{\mathit{n}}\times \left[\widehat{\mathit{n}}\times {\int }_S\mathbf{\nabla }\times {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\cdot \mathit{J}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }\right]\\ & ={\mathit{\eta }}_{\mathrm{s}}\widehat{\mathit{n}}\times \left[\widehat{\mathit{n}}\times {\mathit{H}}^{\mathrm{inc}}\left(\mathit{r}\right)\right]-\widehat{\mathit{n}}\left(\mathit{r}\right)\times {\mathit{E}}^{\mathrm{inc}}\left(\mathit{r}\right),\mathit{r}\in S^{+}\end{array}\end{array}$$

for the unknown current J.

Using the MoM, J is approximated by a linear combination of N basis functions and Eq. (4) is tested using N testing functions. This translates Eq. (4) into a system of linear equations Zi=v, with the vectors i containing the unknown coefficients of the basis functions and v containing the values of the tested right-hand side. The entries of the impedance matrix Z are tested field-integral contributions produced by individual basis functions, by means of 𝓖^f. The strength of the coupling through the impedance matrix between sources and observers within the supports of the nth and mth basis and testing functions, respectively, is indicated by the magnitude of the corresponding matrix element Z_mn.

For electrically large scatterers, the straightforward construction and storage of $\mathbf{Z}\in {\mathbb{C}}^{N\times N}$ are computationally cumbersome. To reduce these costs, algebraic compression techniques are often employed. Partitioning S into a hierarchy of subdomains, the basis and testing functions can be partitioned into hierarchies of clusters. From this follows the partitioning of Z into blocks. A block ${\mathbf{Z}}^{\mathrm{os}}\in {\mathbb{C}}^{N_{\mathrm{o}}\times N_{\mathrm{s}}}$ of Z that corresponds to basis and testing function clusters of sizes N_s and N_o, respectively, is expressed in a compressed form if the disjoint subdomains S^s⊆S and S^o⊆S that correspond to these clusters satisfy some size- and distance-based admissibility criteria. An admissible block can be expressed using its numerical rank-ℛ approximation ${\mathbf{Z}}^{\mathrm{os}}\approx {\mathbf{AB}}^{†}$, where $\mathbf{A}\in {\mathbb{C}}^{N_{\mathrm{o}}\times \mathcal{R}}$ and $\mathbf{B}\in {\mathbb{C}}^{N_{\mathrm{s}}\times \mathcal{R}}$. The rank ℛ depends on the desired accuracy, where the smallest ℛ(τ) that is guaranteed to provide a relative spectral norm error of τ is the number of singular values (SVs) σ_n (sorted in descending order) with $n\in \mathrm{1}\mathrm{\dots }min(N^{\mathrm{o}},N^{\mathrm{s}})$ that are greater than τσ₁, where ${\mathit{\sigma }}_{\mathrm{1}}=\Vert {\mathbf{Z}}^{\mathrm{os}}{\Vert }_{\mathrm{2}}$ is the matrix block's spectral norm. The efficiency of LR-compression-based methods depends on ℛ scaling slower with k than N_s and N_o. In the conventional IE (1), for various geometries, blocks that represent broadside interactions between regions are considered admissible for compression (see Fig. 1c). The ranks for such blocks grow, in general, linearly (in 2-D) with the domain sizes (Bucci and Franceschetti, 1989; Gustafsson, 2025; Gustafsson and Brick, 2026) and, therefore, also with N_s and N_o, preventing reduction of the asymptotic scaling of the memory and computation time through LR compression.

For essentially-convex objects, defined as surfaces with small (sub-wavelength) variation from a convex shape, the IE kernel can be modified to include auxiliary components that can be associated with sources internal to S, without violating the governing equations in the exterior. That is, an auxiliary dyadic 𝓖^aux augments the free-space Green's dyadic to form a modified Green's dyadic (MGD)

$$\begin{array}{}\text{(5)}& \mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })={\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })+{\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}^{\prime }).\end{array}$$

The scattered field is represented by a modified source distribution $\stackrel{\mathrm{̃}}{\mathit{J}}$ that radiates according to the MGD, such that

$$\begin{array}{}\text{(6)}& {\mathit{E}}^{\mathrm{sca}}\left(\mathit{r}\right)={\int }_S\mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })\cdot \stackrel{\mathrm{̃}}{\mathit{J}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }.\end{array}$$

The resulting GSIE is given by

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \widehat{\mathit{n}}\left(\mathit{r}\right)\times {\int }_S\mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })\cdot \stackrel{\mathrm{̃}}{\mathit{J}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }\\ & -{\displaystyle \frac{j{\mathit{\eta }}_{\mathrm{s}}}{k\mathit{\eta }}}\widehat{\mathit{n}}\times \left[\widehat{\mathit{n}}\times {\int }_S\mathbf{\nabla }\times \mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })\cdot \stackrel{\mathrm{̃}}{\mathit{J}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{S}^{\prime }\right]\\ & ={\mathit{\eta }}_{\mathrm{s}}\widehat{\mathit{n}}\times \left[\widehat{\mathit{n}}\times {\mathit{H}}^{\mathrm{inc}}\left(\mathit{r}\right)\right]-\widehat{\mathit{n}}\left(\mathit{r}\right)\times {\mathit{E}}^{\mathrm{inc}}\left(\mathit{r}\right),\mathit{r}\in S^{+}\end{array}\end{array}$$

and can be viewed as an extension of the combined source IE (Mautz and Harrington, 1979). The auxiliary component should be designed carefully to enable the accurate representation, using the MGD, of the scattered field while reducing the effective dimensionality and, by that, increase the rank deficiency. While a formalistic framework for the design of such MGDs is yet to be developed, several design approaches (Sharshevsky et al., 2020; Dahan and Brick, 2024; Kalhöfer et al., 2025; Dahan et al., 2026) were shown to be suitable for 2-D problems.

For 2-D z-invariant problems, the formulation simplifies. For example, for TM excitation and a PEC boundary, Eq. (7) reduces to

$$\begin{array}{}\text{(8)}& {\int }_{S}g(\mathit{r},{\mathit{r}}^{\prime }){\stackrel{\mathrm{̃}}{J}}_z\left({\mathit{r}}^{\prime }\right)\mathrm{d}{s}^{\prime }=-E_z^{\mathrm{inc}}\left(\mathit{r}\right),\mathit{r}\in S^{+}\end{array}$$

where $g(\mathit{r},{\mathit{r}}^{\prime })=g^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })+g^{\mathrm{aux}}(\mathit{r},{\mathit{r}}^{\prime })$ is a scalar modified Green's function, as in Kalhöfer et al. (2025), $g^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })=-\frac{k\mathit{\eta }}{\mathrm{4}}{H}_{\mathrm{0}}^{\left(\mathrm{2}\right)}\left(k\right|\mathit{r}-{\mathit{r}}^{\prime }\left|\right)$, and $H_m^{\left(\mathrm{2}\right)}(\cdot )$ is the mth-order Hankel function of the second kind. For an IBC, the modification applies also to the scattered magnetic field, such that

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& {\int }_S\left[g(\mathit{r},{\mathit{r}}^{\prime })+{\displaystyle \frac{j{\mathit{\eta }}_{\mathrm{s}}}{k\mathit{\eta }}}\widehat{n}\cdot {\mathbf{\nabla }}_{\mathrm{t}}g(\mathit{r},{\mathit{r}}^{\prime })\right]{\stackrel{\mathrm{̃}}{J}}_z\left({\mathit{r}}^{\prime }\right)\mathrm{d}{s}^{\prime }\\ & ={\mathit{\eta }}_{\mathrm{s}}{H}_t^{\mathrm{inc}}\left(\mathit{r}\right)-E_z^{\mathrm{inc}}\left(\mathit{r}\right),\mathit{r}\in S^{+}\end{array}\end{array}$$

where $H_t^{\mathrm{inc}}$ is the tangential component in the direction of $\widehat{t}=\widehat{z}\times \widehat{n}$ and ∇_t is the transverse gradient. For TE excitation, the kernel relating the transverse modified vector source distribution $\stackrel{\mathrm{̃}}{\mathit{J}}$ to the electric field is a 2×2 dyadic 𝓖. The IE reduces to

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\int }_S\left[\widehat{t}\cdot \mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })+{\displaystyle \frac{j{\mathit{\eta }}_{\mathrm{s}}}{k\mathit{\eta }}}\widehat{z}\cdot {\mathbf{\nabla }}_{\mathrm{t}}\times \mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })\right]\cdot \stackrel{\mathrm{̃}}{\mathit{J}}\left({\mathit{r}}^{\prime }\right)\mathrm{d}{s}^{\prime }\\ & =-{\mathit{\eta }}_{\mathrm{s}}{H}_z^{\mathrm{inc}}\left(\mathit{r}\right)-E_t^{\mathrm{inc}}\left(\mathit{r}\right),\mathit{r}\in S^{+}\end{array}\end{array}$$

where, now,

$$\begin{array}{}\text{(11)}& {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })=\left(\mathbf{I}+{\displaystyle \frac{\mathrm{1}}{k^{\mathrm{2}}}}{\mathbf{\nabla }}_{\mathrm{t}}{\mathbf{\nabla }}_{\mathrm{t}}\right)g^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime }).\end{array}$$

Focusing on the design of the multipole-based 𝓖, in the remainder of this work, η_s=0.

To increase the LR compressibility of admissible blocks, $\mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })$ should be designed so that the effective dimensionality of the interactions is reduced. This is achieved by constructing 𝓖^aux to approximately cancel ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })$ in large parts of S. As a result, the broadside components of the interactions between subdomains are greatly attenuated. For interactions between touching subdomains (weak admissibility, e.g. Martinsson and Rokhlin, 2005; Ho and Greengard, 2012; Brick et al., 2014; Sharshevsky et al., 2020), the end-fire interaction components become dominant, leading to a dimensionality reduction in a straightforward sense. This translates to a very slow scaling (almost constant) ℛ(τ), for τ values greater than the shadow level provided by the kernel. For strictly separated domains (strong admissibility, e.g. Zhao et al., 2005; Shaeffer, 2008; Chai and Jiao, 2013), the modified admissible interaction is composed entirely of the greatly diminished broadside component. As near-neighbor and self interactions are not significantly attenuated, this allows for a more aggressive compression of the admissible blocks, following the approach in Dahan and Brick (2024) (powered by the spectral norm estimation in Kelley et al., 2023) to maintain block ranks that are independent of k.

It is, therefore, desired to design 𝓖 to exhibit a deep shadow, that will enable significant compression for broad ranges of τ values. The auxiliary contribution 𝓖^aux should be computable fast for arbitrary $(\mathit{r},{\mathit{r}}^{\prime })$ pairs. In the GSIE body of work, this has been done by computing “stencil” auxiliary contributions that can be rotated and translated for each ${\mathit{r}}^{\prime }\in S$. For vector sources and MGD stencils, two vector auxiliary contributions should be computed. To this end, it is convenient to describe the essentially-convex surface by a coordinate ξ measured along an inner convex hull S_c from some reference point and a coordinate ζ(ξ) representing the distance from S_c (see Fig. 1a). An auxiliary contribution is then attributed to each of two orthogonal components of a vector surface distribution $\widehat{\mathit{\zeta }}\left({\mathit{r}}^{\prime }\right)$ and $\widehat{\mathit{\xi }}\left({\mathit{r}}^{\prime }\right)$, corresponding to the coordinates ξ and ζ, respectively (see Fig. 2). The multipole-based design of these stencils is discussed in the next section.

Figure 2 For some source point ${\mathit{r}}_{\mathrm{0}}^{\prime }$ on the scatterer surface S, the corresponding multipole location ${\mathit{r}}^{\prime \prime }\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ and polar coordinates of the observation point r relative to both ${\mathit{r}}_{\mathrm{0}}^{\prime }\in S$, and ${\mathit{r}}^{\prime \prime }\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ are shown. Note that r, ${\mathit{r}}_{\mathrm{0}}^{\prime }$, and ${\mathit{r}}^{\prime \prime }\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ are defined with respect to a common origin (not shown), typically located at the center of the scattering object. Additionally, in magenta, the local directions $\widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ and $\widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ of the coordinate system with respect to the inner convex hull S_c are illustrated. The design parameters of the modified kernel are highlighted in blue. Dashed black rays show the directions where the kernel at infinity is forced to zero and the same happens in the gray directions due to the symmetry of the construction method.

3 Kernel design

The multipole contribution ${\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}^{\prime })$ to the MGD is designed following similar geometric considerations to those from the TM case (Kalhöfer et al., 2025). For a point ${\mathit{r}}^{\prime }\in S$, the multipole contributions are defined to emanate from a point ${\mathit{r}}^{\prime \prime }\left({\mathit{r}}^{\prime }\right)={\mathit{r}}^{\prime }-\mathit{\upsilon }\widehat{\mathit{\zeta }}\left({\mathit{r}}^{\prime }\right)$ located a distance υ towards S_c from r^′. Let ${\mathit{r}}_{\mathrm{0}}^{\prime }\in S$ denote a source location for which a stencil MGD is designed. Following the approach in Dahan et al. (2026), it is constructed of two vector contributions attributed to $\widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$- and $\widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$-oriented surface current dipoles.

The multipole contributions are assumed to correspond to z-oriented magnetic currents. As such, they can be derived from z-oriented electric vector potentials. Let $\widehat{\mathit{\chi }}$ denote either $\widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ or $\widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$. The corresponding vector potentials are denoted ${\mathit{F}}_{\mathit{\chi }}\left(\mathit{r}\right)=\widehat{\mathit{z}}{F}_{\mathit{\chi }}\left(\mathit{r}\right)$, $\mathit{\chi }\in \mathit{\{}\mathit{\xi },\mathit{\zeta }\mathit{\}}$. They can be written as series of cylindrical wave functions in a polar system with a radial coordinate $\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })=|\mathit{r}-{\mathit{r}}^{\prime \prime }({\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right|$ and an azimuthal $\mathit{\phi }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ angle defined between $\widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ and $\mathit{r}-{\mathit{r}}^{\prime \prime }\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ (see Fig. 2). The expansions take the form

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& F_{\mathit{\chi }}\left(\mathit{r}\right)={\displaystyle \frac{\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{0}}^{\mathrm{\infty }}{H}_m^{\left(\mathrm{2}\right)}\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\\ & \left[a_{\mathit{\chi }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)+b_{\mathit{\chi }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\right]\end{array}\end{array}$$

where a_χm and b_χm are dimensionless complex amplitudes. The corresponding vectors forming the MGD are given by ${\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\chi }}={\mathbf{\nabla }}_{\mathrm{t}}\times {\mathit{F}}_{\mathit{\chi }}\left(\mathit{r}\right)$. It is convenient to express them in terms of their radial and azimuthal components $\widehat{\mathit{\rho }}$ and $\widehat{\mathit{\phi }}$, defined with respect to ${\mathit{r}}^{\prime \prime }\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$, such that

$$\begin{array}{}\text{(13)}& {\displaystyle }\begin{array}{rl}& \widehat{\mathit{\rho }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\chi }}={\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{1}}^{\mathrm{\infty }}{\displaystyle \frac{H_m^{\left(\mathrm{2}\right)}\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}{k\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}\\ & \left[a_{\mathit{\chi }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)-b_{\mathit{\chi }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\right]m\end{array}\text{(14)}& {\displaystyle }\begin{array}{rl}& \widehat{\mathit{\phi }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\chi }}={\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{0}}^{\mathrm{\infty }}{H}_m^{\left(\mathrm{2}\right)\prime }\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\\ & \left[a_{\mathit{\chi }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)+b_{\mathit{\chi }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\right]\end{array}\end{array}$$

In the design of the stencil MGD, a_χm and b_χm should be set so that the magnitudes of the components of $\mathcal{G}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ inside the scatterer and to faraway segments of S are significantly reduced compared to those of ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$. To this end, let us express ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ using similar vector components. With $\stackrel{\mathrm{̃}}{\mathit{\rho }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })=|\mathit{r}-{\mathit{r}}_{\mathrm{0}}^{\prime }|$ and $\stackrel{\mathrm{̃}}{\mathit{\phi }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ denoting the angle between $\widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)$ and $\mathit{r}-{\mathit{r}}_{\mathrm{0}}^{\prime }$ (see Fig. 2), 𝓖^f can be written as

$$\begin{array}{}\text{(15)}& {\displaystyle }\begin{array}{rl}& \widehat{\stackrel{\mathrm{̃}}{\mathit{\rho }}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\cos \left(\stackrel{\mathrm{̃}}{\mathit{\phi }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right){\displaystyle \frac{H_{\mathrm{1}}^{\left(\mathrm{2}\right)}\left(k\stackrel{\mathrm{̃}}{\mathit{\rho }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}{k\stackrel{\mathrm{̃}}{\mathit{\rho }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}\end{array}\text{(16)}& {\displaystyle }\begin{array}{rl}& \widehat{\stackrel{\mathrm{̃}}{\mathit{\phi }}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sin \left(\stackrel{\mathrm{̃}}{\mathit{\phi }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)H_{\mathrm{1}}^{\left(\mathrm{2}\right)\prime }\left(k\stackrel{\mathrm{̃}}{\mathit{\rho }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}& \widehat{\stackrel{\mathrm{̃}}{\mathit{\rho }}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sin \left(\stackrel{\mathrm{̃}}{\mathit{\phi }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right){\displaystyle \frac{H_{\mathrm{1}}^{\left(\mathrm{2}\right)}\left(k\stackrel{\mathrm{̃}}{\mathit{\rho }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}{k\stackrel{\mathrm{̃}}{\mathit{\rho }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}\end{array}\text{(18)}& {\displaystyle }\begin{array}{rl}& \widehat{\stackrel{\mathrm{̃}}{\mathit{\phi }}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\cos \left(\stackrel{\mathrm{̃}}{\mathit{\phi }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)H_{\mathrm{1}}^{\left(\mathrm{2}\right)\prime }\left(k\stackrel{\mathrm{̃}}{\mathit{\rho }}\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right).\end{array}\end{array}$$

As shown in Kalhöfer et al. (2025), the multipole amplitudes can be determined by trying to minimize the MGD magnitude in a far-field observation sector. For large ρ, the approximations $\stackrel{\mathrm{̃}}{\mathit{\phi }}\approx \mathit{\phi }$, $\stackrel{\mathrm{̃}}{\mathit{\rho }}\approx \mathit{\rho }-\mathit{\upsilon }\cos \left(\mathit{\phi }\right)\approx \mathit{\rho }$, and $H_m\left(k\mathit{\rho }\right)\approx j^m\sqrt{\mathrm{2}/\mathit{\pi }k\mathit{\rho }}{e}^{j(\mathit{\pi }/\mathrm{4}-k\mathit{\rho })}$ (DLMF, 2025, Eq. (10.2.6)) hold and Eqs. (13)–(18) can be combined and written as

$$\begin{array}{}\text{(19)}& {\displaystyle }\begin{array}{rl}& \mathcal{G}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\stackrel{\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\to \mathrm{\infty }}{⟶}{\displaystyle \frac{\mathit{\eta }{e}^{j(\mathit{\pi }/\mathrm{4}-k\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}}{\mathrm{2}\sqrt{\mathit{\lambda }\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}}\widehat{\mathit{\phi }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\\ & (\sum _{m=\mathrm{0}}^{\mathrm{\infty }}{j}^{m-\mathrm{1}}\left[a_{\mathit{\zeta }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)+b_{\mathit{\zeta }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\right]\\ & +e^{jk\mathit{\upsilon }\cos \left(\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}\sin \left(\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right))\end{array}\text{(20)}& {\displaystyle }\begin{array}{rl}& \mathcal{G}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\stackrel{\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\to \mathrm{\infty }}{⟶}{\displaystyle \frac{\mathit{\eta }{e}^{j(\mathit{\pi }/\mathrm{4}-k\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}}{\mathrm{2}\sqrt{\mathit{\lambda }\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}}\widehat{\mathit{\phi }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\\ & (\sum _{m=\mathrm{0}}^{\mathrm{\infty }}{j}^{m-\mathrm{1}}\left[a_{\mathit{\xi }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)+b_{\mathit{\xi }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\right]\\ & -e^{jk\mathit{\upsilon }\cos \left(\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}\cos \left(\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right))\end{array}\end{array}$$

To minimize $\mathcal{G}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ in a sector defined symmetrically with respect to φ, the components of ${\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ should inherit the symmetry or anti-symmetry of their counterparts in ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$. From this follows that $a_{\mathit{\zeta }m}=b_{\mathit{\xi }m}=\mathrm{0}$ for all m. The design problem is further simplified by setting $a_{\mathit{\xi }m}=b_{\mathit{\zeta }m}=\mathrm{0}$ for all m≥M, with the largest auxiliary multipole integer order parameter M. The remaining 2M−1 amplitudes are computed by solving a minimization problem. Here, this is done by trying to minimize the least-squares norm of the residual of $M^{\prime }\ge M$ equations, defined by setting the tangent far-field components of $\mathcal{G}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })$ to zero at M^′ angular sampling points ${\mathit{\phi }}_{m^{\prime }}$. In this work, these angles are uniformly spaced such that

$$\begin{array}{}\text{(21)}& {\mathit{\phi }}_{m^{\prime }}=\mathit{\pi }-{\displaystyle \frac{m^{\prime }-\mathrm{1}}{M^{\prime }-\mathrm{1}}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\phi }}{\mathrm{2}}},m^{\prime }\in \mathrm{1}\mathrm{\dots }{M}^{\prime }\end{array}$$

where Δφ is an opening angle design parameter (see Fig. 2). Two minimization problems are solved, one for each of the components, with Eqs. (19) and (20) translated to

$$\begin{array}{}\text{(22)}& {\displaystyle }\begin{array}{rl}& \sum _{m=\mathrm{1}}^{M-\mathrm{1}}{j}^{m+\mathrm{1}}\sin \left(m{\mathit{\phi }}_{m^{\prime }}\right)b_{\mathit{\zeta }m}\\ & =\phantom{-}{e}^{jk\mathit{\upsilon }\cos \left({\mathit{\phi }}_{m^{\prime }}\right)}\sin \left({\mathit{\phi }}_{m^{\prime }}\right),m^{\prime }\in \mathrm{2}\mathrm{\dots }{M}^{\prime }\end{array}\text{(23)}& {\displaystyle }\begin{array}{rl}& \sum _{m=\mathrm{0}}^{M-\mathrm{1}}{j}^{m+\mathrm{1}}\cos \left(m{\mathit{\phi }}_{m^{\prime }}\right)a_{\mathit{\xi }m}\\ & =-e^{jk\mathit{\upsilon }\cos \left({\mathit{\phi }}_{m^{\prime }}\right)}\cos \left({\mathit{\phi }}_{m^{\prime }}\right),m^{\prime }\in \mathrm{1}\mathrm{\dots }{M}^{\prime }\end{array}\end{array}$$

For $M^{\prime }=M$ and sufficiently low M values, the minimization problems can be solved exactly. Once solved, the stencil auxiliary contributions can be written as

$$\begin{array}{}\text{(24)}& {\displaystyle }\begin{array}{rl}& \widehat{\mathit{\rho }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{1}}^{M-\mathrm{1}}{b}_{\mathit{\zeta }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right){\displaystyle \frac{H_m^{\left(\mathrm{2}\right)}\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}{k\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}m\end{array}\text{(25)}& {\displaystyle }\begin{array}{rl}& \widehat{\mathit{\phi }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{1}}^{M-\mathrm{1}}{b}_{\mathit{\zeta }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)H_m^{\left(\mathrm{2}\right)\prime }\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\end{array}\text{(26)}& {\displaystyle }\begin{array}{rl}& \widehat{\mathit{\rho }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{1}}^{M-\mathrm{1}}{a}_{\mathit{\xi }m}\sin \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right){\displaystyle \frac{H_m^{\left(\mathrm{2}\right)}\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)}{k\mathit{\rho }(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })}}m\end{array}\text{(27)}& {\displaystyle }\begin{array}{rl}& \widehat{\mathit{\phi }}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot {\mathcal{G}}^{\mathrm{aux}}(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}_{\mathrm{0}}^{\prime }\right)\\ & =-{\displaystyle \frac{k\mathit{\eta }}{\mathrm{4}}}\sum _{m=\mathrm{0}}^{M-\mathrm{1}}{a}_{\mathit{\xi }m}\cos \left(m\mathit{\phi }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)H_m^{\left(\mathrm{2}\right)\prime }\left(k\mathit{\rho }\right(\mathit{r},{\mathit{r}}_{\mathrm{0}}^{\prime }\left)\right)\end{array}\end{array}$$

and employed, via translation and rotation, for arbitrary ${\mathit{r}}^{\prime }\in S$.

Figure 3 log ₁₀ of kernel component magnitudes normalized to $\mathit{\eta }/\mathit{\lambda }$ for the EFIE and the multipole-based GSIE with υ=1λ, Δφ=120°, and $M=M^{\prime }=\mathrm{8}$. (a) ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}^{\prime }\right)$, (b) $\mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })\cdot \widehat{\mathit{\zeta }}\left({\mathit{r}}^{\prime }\right)$, (c) ${\mathcal{G}}^{\mathrm{f}}(\mathit{r},{\mathit{r}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}^{\prime }\right)$, and (d) $\mathcal{G}(\mathit{r},{\mathit{r}}^{\prime })\cdot \widehat{\mathit{\xi }}\left({\mathit{r}}^{\prime }\right)$.

4 Examples

The numerical examples in this section are used to examine the key properties of the proposed 𝓖 and their dependence on its design parameters. First, 𝓖 is depicted for a representative set of construction parameters. Then, the compressibility of moment-matrix off-diagonal blocks and accuracy of the moment solution achieved using 𝓖 are studied as a function of the parameters. Finally, the error contractibility is studied through the solution of two representative scattering problems.

In the first example, 𝓖 was computed for the parameters υ=1λ, Δφ=120°, and $M=M^{\prime }=\mathrm{8}$. Henceforth, these construction parameters are used where not otherwise stated. For the two fundamental components, Fig. 3 compares 𝓖 to 𝓖^f. It can be seen that 𝓖 exhibits a deep shadow in the negative $\widehat{\mathit{\zeta }}$-direction and a strong field in the $\widehat{\mathit{\xi }}$-direction and positive $\widehat{\mathit{\zeta }}$-direction. The shadow depth is very similar to that achieved by the multipole-based TM-GSIE for the same construction parameters (Kalhöfer et al., 2025). This behavior is responsible for the reduced dimensionality and enhanced rank deficiency. This can also be seen in Fig. 4, showing the Green's dyadic and MGD components, tested in the tangent direction to the circle along its perimeter, for a source located on a circle of radius R=16λ. The depth of the MGD components shadow at the points across from the source suggests the strong accentuation of broadside interactions.

Figure 4 ξ-components for the kernels in Fig. 3, for a circular contour with R=16λ for ξ- and ζ-directions of a source located at φ^′ on the contour.

In the next set of examples, the influence of the design parameters on the performance of the MGD is studied. To this end, the equations are discretized using the MoM with triangular basis and testing functions defined on polygonal meshes. Where not otherwise specified, a mesh edge length of $h\approx \mathit{\lambda }/\mathrm{40}$ and a four-point Gaussian quadrature rule were used. First, this is applied to a circular scatterer with radius R=16λ. To analyze how different kernels influence the compressibility of matrix blocks, the SVs of the off-diagonal block that represents the interaction between a quarter of the circle and its remainder are computed. The accuracy of the method is measured by comparing the results for the scattered electric field on a circle of radius $R+\mathit{\lambda }/\mathrm{4}$ to the conventional solution.

Figure 5 shows how the construction parameters M (top), Δφ (middle), and υ (bottom) influence the compressibility (left) and accuracy (right). As in the TM case, the shadow gets deeper for greater M values, providing greater rank deficiency even at lower threshold values. Unlike in the TM case, the compressibility does not depend strongly on Δφ. This can be explained by the more broadside nature of the vector elemental sources, which is practically eliminated entirely by its attenuation in a narrow sector. The relative error levels depend only weakly on M and Δφ. The error is roughly 10⁻⁴ for $h=\mathit{\lambda }/\mathrm{40}$. It is hypothesized that increasing Δφ does not change the MGD significantly, that greater M are mostly relevant for deepening the shadow near the edges of the sector, and that both have little influence on key features of the near field behavior of the MGD, which determines the error level. For comparison, using $h=\mathit{\lambda }/\mathrm{80}$ (and R=8λ), error levels of the order of 10⁻⁵ are observed, suggesting already the controllability of the error with h. The parameter υ has little influence on the compressibility, which shows mostly for υ<λ for rather low singular values. These differences are likely to not show for most practical compression threshold values. The accuracy depends more strongly on υ, with a rapid increase in the error for $\mathit{\upsilon }<\mathit{\lambda }/\mathrm{2}$. For larger υ values, the error stagnates at the limit dictated by the discretization. Similar behavior was observed by Kalhöfer et al. (2025) for the TM case.

Figure 5 (a, c, e) SVs for various values of: (a) M, (c) Δφ, and (e) υ. (b, d, f) Relative difference as a function of the parameter (b) M, (d) Δφ, and (f) υ.

Lastly, the proposed GSIE was solved using the MoM for two structures illuminated by unit plane wave incident fields: a circular cylinder of radius R=16λ and a corrugated circular cylinder with minimum radius R=16λ, corrugation depth Δ_ζ=0.4λ and a cyclic corrugation with 84 cycles along the perimeter. The GSIE parameters were set to M=8, Δ_φ=120°, and υ=1λ. The results for the scattered field computed at a circle of radius $R+{\mathrm{\Delta }}_{\mathit{\zeta }}+\mathit{\lambda }/\mathrm{4}$ are shown in Fig. 6. The absolute errors in the computed fields are assessed through comparison of the GSIE solution, for various values of h, to its EFIE counterparts. For both examples, the error decreases steadily with h.

Figure 6 GSIE scattered field (red) and difference from the EFIE counterpart for different values (shades of blue) of h for a unit amplitude incident plane wave. The direction of incidence, the scatterer, a magnified section of its surface, and the observation circle (dotted) are shown in the insets. (a) Circular scatterer and (b) corrugated cylinder.

5 Conclusions

The multipole-based GSIE construction method was extended, from the previously proposed scalar (TM) IE formulation in Kalhöfer et al. (2025) to a dyadic (TE) one, using the principles in Dahan et al. (2026). The dyadic kernel design was presented in detail and its properties relevant for fast-solver design and performance were analyzed. The MGD was shown to provide the desired reduction of dimensionality of interactions and enhanced rank-deficiency, while maintaining error controllability of the MoM solution with the mesh density. The proposed formulation is the first of its kind that does not require kernel-tailored methods for removal of computational bottlenecks associated with the MGD computation. By avoiding the usage of large shield surfaces it overcomes many of the geometric restrictions of its predecessors and can be used effectively even for objects that include very narrow regions. The MGD components can be effortlessly repurposed for the analysis of TM and TE scattering by arbitrary impedance boundary impenetrable objects in 2-D. The 2-D dyadic derivation is also a critical stepping stone toward the development of full vector 3-D GSIEs. Further research on other extensions includes the development of GSIE-based formulations for penetrable and non-convex objects. Rigorous investigation of the spectral characteristics and various well-posedness aspects of GSIE formulations is expected to be more convenient for 2-D multipole formulations. The TE formulation is of particular relevance for the study of GSIE conditioning and preconditioning.