# Eddy Current Modeling in Multiply Connected Regions via a Full-Wave Solver Based on the Quasi-Helmholtz Projectors

Tiffany Chhim, Adrien Merlini, Lyes Rahmouni, John Erick Ortiz Guzman, Francesco Andriulli

<figure><img src="image_1.png" /></figure>

## To cite this version:

Tiffany Chhim, Adrien Merlini, Lyes Rahmouni, John Erick Ortiz Guzman, Francesco Andriulli. Eddy Current Modeling in Multiply Connected Regions via a Full-Wave Solver Based on the Quasi-Helmholtz Projectors. IEEE Open Journal of Antennas and Propagation, 2020, 1, pp.534-548. 10.1109/OJAP.2020.3027186. [hal-03401788](https://hal.science/hal-03401788)

HAL Id: hal-03401788

https://hal.science/hal-03401788

Submitted on 26 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Received 31 July 2020; revised 28 August 2020; accepted 8 September 2020. Date of publication 29 September 2020; date of current version 19 October 2020.

Digital Object Identifier 10.1109/OJAP.2020.3027186

# Eddy Current Modeling in Multiply Connected Regions via a Full-Wave Solver Based on the Quasi-Helmholtz Projectors

TIFFANY L. CHHIM<fnref n="1" />, ADRIEN MERLINI<fnref n="1" />,<fnref n="2" /> (Member, IEEE), LYES RAHMOUNI<fnref n="1" />, JOHN ERICK ORTIZ GUZMAN<fnref n="1" />,<fnref n="3" />, AND FRANCESCO P. ANDRIULLI<fnref n="1" />,⁰ (Senior Member, IEEE)

CORRESPONDING AUTHOR: F. P. ANDRIULLI (e-mail: francesco.andriulli@polito.it)

This work was supported in part by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under Grant 724846 (Project 321), and in part by the Italian Ministry of Education, University and Research within the Program under Grant PRIN2017, Grant 2017HZJXSZ, and Grant CUP:E6419002530001.

**ABSTRACT** Eddy currents are central to several industrial applications and there is a strong need for their efficient modeling. Existing eddy current solution strategies are based on a quasi-static approximation of Maxwell's equations for lossy conducting objects and thus their applicability is restricted to low frequencies. On the other hand, available full-wave solvers such as the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) equation become highly ill-conditioned and inaccurate in eddy current settings. This work presents a new well-conditioned and stable full-wave formulation which encompasses the simulation of eddy currents. Our method is built upon the PMCHWT equation and thus remains valid over the entire frequency range. Moreover, our scheme is also compatible with structures containing holes and handles (multiply connected geometries). The effectiveness of quasi-Helmholtz projectors is leveraged to obtain a versatile solver, which is computationally efficient and allows for a seamless transition between low and high frequencies. The stability and accuracy of the new method are demonstrated both theoretically and through numerical experiments on canonical and realistic structures.

**INDEX TERMS** Eddy currents, preconditioning, full-wave, multiply connected, quasi-Helmholtz decomposition.

## I. INTRODUCTION

THE BOUNDARY Element Method (BEM) is among the most widespread techniques to simulate scattering and radiation phenomena in electromagnetics. Its main benefit is the necessity to discretize only the surfaces separating the different media, as opposed to the entire volume, thus greatly reducing the size of the problem. Integral equation approaches have been used to create a wide family of formulations to accommodate different scenarios. For instance the Electric Field Integral Equation (EFIE) [1] is commonly used to simulate perfect electric conductors, while the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT)

equation [2]–[4] has been used to handle dielectric and conducting bodies.

The simulation of eddy currents, which are generated inside conductors in the vicinity of a time-varying magnetic field, is of major interest for many industrial applications [5]. In particular, induced eddy currents are widely employed for non destructive testing in various manufacturing areas to detect the presence of material defects [6]. As of now, the modeling of this type of scenarios requires ad hoc solvers built on the quasi-static approximation of Maxwell's equations after neglecting displacement currents [7], [8]. Unfortunately, this simplification greatly reduces their range

<footnote n="1">¹Department of Electronics and Telecommunications, Politecnico di Torino, 10129 Turin, Italy</footnote>

<footnote n="2">²Microwaves Department, IMT Atlantique, 29238 Brest, France</footnote>

<footnote n="3">³Electronic Engineering Department, Corporación Universitaria Autónoma de Nariño, Pasto 520002, Colombia</footnote>
of validity. Since the quasi-static approximation is only valid at low frequencies, these solvers are particularly ill-suited for electrically large conductors [9], [10]. In addition, the simulation of multiply connected geometries usually requires further preprocessing such as the introduction of cutting surfaces or structural loops detection to guarantee the uniqueness of the solution [11], [12].

On the other hand, the PMCHWT formulation is valid to model lossy conductors for moderate conductivities, but for high conductivities the discretization density must be increased to maintain the accuracy of the integrals [13] and numerical instabilities may occur. Lamentably, in fact, the PMCHWT formulation suffers from a severe breakdown at low frequencies including in the eddy current regime. More precisely, the finite electrical conductivity of the material simulated at low frequencies causes a critical ill-conditioning of the equation [14], which in turn impacts the performance of iterative solvers typically used in fast methods by slowing down or even preventing altogether their convergence. In addition, at very low frequencies, a detrimental loss of accuracy in the solution also occurs. This is the counterpart, for frequency dependent material parameters, of the well-known low frequency breakdown of perfect conductors or purely dielectric objects, which has been treated thoroughly in the literature using a variety of strategies [15]–[17]. Solutions to the low frequency breakdown in lossy materials have also been proposed, often leveraging loop-star/loop-tree decompositions [18] that however have a significant computation overhead caused by their conditioning and the detection of global cycles. In addition, the augmented EFIE [19] has successfully been extended to lossy conductors [20] and inhomogeneous media [21] but these extensions require additional matrices to be computed and stored, thus increasing their computational cost. Among the different approaches, the use of quasi-Helmholtz projectors [22]–[24] offers many benefits over previous techniques, including an improved stability, an implicit handling of multiply connected geometries, and a compatibility with fast solvers operating with quasi-linear computational complexity [25], [26]. However, the approaches developed in [22], [24] are only available for perfectly conducting objects. Similarly, the formulation introduced in [23] has been derived for purely dielectric (lossless) materials but the proposed low frequency regularizer is not applicable to the lossy case because it cannot rescale separately the upper and lower diagonal blocks of the PMCHWT which do not follow the same frequency behavior under eddy current conditions. This work extends the use of quasi-Helmholtz projectors to the eddy current regime, which requires an ad hoc analysis and preconditioning strategy due to the peculiar low frequency behavior associated to the complex permittivity of the material. Separate rescaling strategies need to be established leveraging both standard and dual operators to overcome the difference in scaling of the upper and lower diagonal blocks in such a way that: (i) both diagonal blocks are cured, (ii) the off-diagonal blocks are not corrupted by a new source of breakdown, (iii) no null

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Definition of relevant quantities and parameters. An exciting electromagnetic field impinges on the conductor immersed in the air medium.</figcaption></figure>

space is introduced by the rescaling, and (iv) the accuracy is preserved at very low frequencies.

The proposed full-wave solver is capable of handling eddy current modeling while effectively attaining (i) to (iv). In fact it shows a stable conditioning, is free from the loss of significant digits at very low frequencies, steadily shifts between low and high frequencies, and is compatible with both simply and multiply connected geometries, while maintaining a low computational complexity. Very preliminary results have been presented in the conference contribution [27].

The paper is organized as follows. The background and notation are set in Section II. A thorough analysis of the PMCHWT equation low frequency behavior in our setting is detailed in Section III, and is followed in Section IV by the presentation of the new strategy employed to obtain a robust full-wave formulation. Section V contains additional details related to the implementation. The development is then supported by several numerical results in Section VI, after which Section VII concludes the paper.

## II. BACKGROUND AND NOTATION

Let $\Omega_1 \subset \mathbb{R}^3$ be a lossy conductor with boundary $\Gamma = \partial\Omega_1$ and outward pointing normal $\hat{n}$, residing in the outside medium $\Omega_0 = \mathbb{R}^3 \setminus \bar{\Omega}_1$. Throughout this article, the indices $\{0, 1\}$ represent the exterior air medium of conductivity $\sigma_0 = 0$ and the interior medium, respectively, as illustrated in Figure 1. The boundary $\Gamma$ can be either simply or multiply connected, i.e., it can contain holes and handles. The conductor is characterized by its constant permeability $\mu_1 = \mu_0\mu_r$, constant conductivity $\sigma_1$, and complex permittivity $\epsilon_1 = \epsilon_0\epsilon'_r - j\sigma_1/\omega$, where $\epsilon'_r$ is the real-valued relative permittivity and $\omega$ is the angular frequency. The Electric Field Integral Operator (EFIO) is defined as

$$ \mathcal{T}_k = -jk\mathcal{T}_{A,k} + \frac{1}{jk}\mathcal{T}_{\phi,k} \quad (1) $$

$$ (\mathcal{T}_{A,k}f)(\mathbf{r}) = \hat{\mathbf{n}} \times \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}')f(\mathbf{r}') \, d\mathbf{r}' \quad (2) $$

$$ (\mathcal{T}_{\phi,k}f)(\mathbf{r}) = \hat{\mathbf{n}} \times \nabla \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}') \nabla' \cdot f(\mathbf{r}') \, d\mathbf{r}' \quad (3) $$

and the Magnetic Field Integral Operator (MFIO) is defined as

$$ (\mathcal{K}_k f)(\mathbf{r}) = \hat{\mathbf{n}} \times \int_{\Gamma} \nabla G_k(\mathbf{r}, \mathbf{r}') \times f(\mathbf{r}') \, d\mathbf{r}' \quad (4) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Notation used for the definition of RWG functions and Loop-Star transformation matrices.</figcaption></figure>

where $k = \omega\sqrt{\mu\epsilon}$ is the wave number, and

$$G_k(\mathbf{r}, \mathbf{r}') = \frac{e^{-j\kappa R}}{4\pi R} \qquad (5)$$

is the Green's function with $R = \|\mathbf{r} - \mathbf{r}'\|$. The eddy current regime is characterized by the following conditions [10], [28]

$$\left\{ \begin{array}{l} \omega\epsilon_0 \ll \sigma_1 \\ L\omega\sqrt{\mu_0\epsilon_0} \ll 1, \end{array} \right. \qquad (6)$$

where $L$ is the characteristic size of the conductor. In such conditions, the wave number inside the object is

$$k_1 \approx \sqrt{-j\omega\sigma_1\mu_1} \approx \frac{1-j}{\delta} \qquad (7)$$

where $\delta = \sqrt{2/(\omega\sigma_1\mu_1)}$ represents the skin depth and influences the concentration of current densities near the surface of the conductor [10].

The PMCHWT integral equation reads [2]–[4]

$$\left( \begin{array}{cc} \eta_0 \mathcal{T}_{k_0} + \eta_1 \mathcal{T}_{k_1} & -(\mathcal{K}_{k_0} + \mathcal{K}_{k_1}) \\ \mathcal{K}_{k_0} + \mathcal{K}_{k_1} & \frac{1}{\eta_0} \mathcal{T}_{k_0} + \frac{1}{\eta_1} \mathcal{T}_{k_1} \end{array} \right) \left( \begin{array}{c} j_s \\ \mathbf{m}_s \end{array} \right) = \left( \begin{array}{c} -\hat{\mathbf{n}} \times \mathbf{E}^i \\ -\hat{\mathbf{n}} \times \mathbf{H}^i \end{array} \right) \qquad (8)$$

where the unknowns are the electric and magnetic surface current densities $j_s = \hat{\mathbf{n}} \times \mathbf{H}$ and $\mathbf{m}_s = -\hat{\mathbf{n}} \times \mathbf{E}$, and where $\eta_{0,1} = \sqrt{\mu_{0,1}/\epsilon_{0,1}}$ is the characteristic impedance of the exterior or interior medium.

The next step towards obtaining a boundary element matrix system is to discretize (8) by approximating the surface $\Gamma$ with a mesh of planar triangular elements on which the divergence-conforming Rao-Wilton-Glisson (RWG) basis functions [1] are defined between pairs of adjacent triangles $c_n^\pm$ as

$$f_n(\mathbf{r}) = \begin{cases} \frac{\mathbf{r} - \mathbf{r}_n^+}{2A_n^+}, & \mathbf{r} \in c_n^+ \\ \frac{\mathbf{r}_n^- - \mathbf{r}}{2A_n^-}, & \mathbf{r} \in c_n^- \end{cases} \qquad (9)$$

where $A_n^\pm$ is the area of the triangle $c_n^\pm$ and $\mathbf{r}_n^\pm$ is the position vector of the vertex that does not belong to the common edge (Figure 2). The functions are normalized so that the flux integral through their defining edges equals one [22].

These functions are then used to expand the unknown current densities as

$$j_s \approx \sum_{n=1}^{N_e} j_n \mathbf{f}_n \qquad (10)$$

$$\mathbf{m}_s \approx \sum_{n=1}^{N_e} m_n \mathbf{f}_n \qquad (11)$$

where $N_e$ is the number of internal edges in the mesh. The equations are finally tested with curl-conforming rotated RWG functions $\{\hat{\mathbf{n}} \times \mathbf{f}_n\}$ to produce the linear system

$$\left( \mathcal{T}_{\text{upper}} \quad -\mathcal{K} \quad \mathcal{T}_{\text{lower}} \right) \left( \begin{array}{c} j \\ m \end{array} \right) = \left( \begin{array}{c} e \\ h \end{array} \right), \qquad (12)$$

where

$$\mathcal{K} = \mathcal{K}_{k_0} + \mathcal{K}_{k_1} \qquad (13)$$

$$\begin{aligned} \mathcal{T}_{\text{upper}} &= \eta_0 \mathcal{T}_{k_0} + \eta_1 \mathcal{T}_{k_1} \\ &= -j k_0 \eta_0 \mathcal{T}_{A,k_0} - j k_1 \eta_1 \mathcal{T}_{A,k_1} \\ &\quad + \frac{\eta_0}{j k_0} \mathcal{T}_{\Phi,k_0} + \frac{\eta_1}{j k_1} \mathcal{T}_{\Phi,k_1} \end{aligned} \qquad (14)$$

$$\begin{aligned} \mathcal{T}_{\text{lower}} &= \frac{1}{\eta_0} \mathcal{T}_{k_0} + \frac{1}{\eta_1} \mathcal{T}_{k_1} \\ &= -j \frac{k_0}{\eta_0} \mathcal{T}_{A,k_0} - j \frac{k_1}{\eta_1} \mathcal{T}_{A,k_1} \\ &\quad + \frac{1}{j k_0 \eta_0} \mathcal{T}_{\Phi,k_0} + \frac{1}{j k_1 \eta_1} \mathcal{T}_{\Phi,k_1} \end{aligned} \qquad (15)$$

and

$$(\mathcal{K}_k)_{mn} = \langle \hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{K}_k(\mathbf{f}_n) \rangle \qquad (16)$$

$$(\mathcal{T}_k)_{mn} = \langle \hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{T}_k(\mathbf{f}_n) \rangle \qquad (17)$$

$$(\mathcal{T}_{A,k})_{mn} = \langle \hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{T}_{A,k}(\mathbf{f}_n) \rangle \qquad (18)$$

$$(\mathcal{T}_{\Phi,k})_{mn} = \langle \hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{T}_{\Phi,k}(\mathbf{f}_n) \rangle \qquad (19)$$

$$(j)_n = j_n \qquad (20)$$

$$(m)_n = m_n \qquad (21)$$

$$(e)_m = \langle \hat{\mathbf{n}} \times \mathbf{f}_m, -\hat{\mathbf{n}} \times \mathbf{E}^i \rangle \qquad (22)$$

$$(h)_m = \langle \hat{\mathbf{n}} \times \mathbf{f}_m, -\hat{\mathbf{n}} \times \mathbf{H}^i \rangle \qquad (23)$$

in which $\langle u, v \rangle = \int_{\Gamma} u \cdot v \, d\Gamma$.

In the next sections, we will also have to use the Loop-Star transformation matrices. For the sake of completeness and to set the notation we define these matrices below very briefly, but the reader interested in more details should refer to [22] and references therein. We denote with $\Lambda$, $H$, and $\Sigma$ the transformation matrices from the Loop, global Loop, and Star subspaces, respectively, to the RWG space. We will also denote with $\mathcal{A} = (\Lambda \ H \ \Sigma)$ the matrix used to perform the complete decomposition. The matrices $\Lambda$ and $\Sigma$ are defined as

$$(\Lambda)_{ij} = \begin{cases} 1 & \text{if node } j \text{ equals } \mathbf{v}_i^+ \\ -1 & \text{if node } j \text{ equals } \mathbf{v}_i^- \\ 0 & \text{otherwise} \end{cases} \qquad (24)$$

$$(\Sigma)_{ij} = \begin{cases} 1 & \text{if cell } j \text{ equals } c_i^+ \\ -1 & \text{if cell } j \text{ equals } c_i^- \\ 0 & \text{otherwise} \end{cases} \qquad (25)$$

where we have used the notation from Figure 2. To ensure linear independence for both matrices, columns should be
appropriately eliminated: for each connected component of the geometry, one corresponding column must always be removed from $Σ$, and when the component is closed, another one must be eliminated from $Λ$.

The matrix $H$ corresponds to the harmonic subspace [22], of dimension $N_{holes}+2N_{handles}$, and it cannot be described as simply as $Λ$ and $Σ$; the reader is referred to [15], [22], [29] for further details on its construction. An explicit definition of $H$ will not be given here, since the retrieval of the global loops is a costly operation that would compromise the overall complexity of the solver, and therefore will be avoided in the scheme we present here. In fact in this work we will stabilize the PMCHWT in the eddy current regime by using quasi-Helmholtz projectors as defined in [22] which are capable of bypassing the explicit definition of $H$, making the expensive building of the harmonic subspace unnecessary. The quasi-Helmholtz projectors are defined as

$$ 
\begin{aligned}
P^{\Sigma} &= \Sigma (\Sigma^{\mathrm{T}} \Sigma)^{+} \Sigma^{\mathrm{T}} && (26) \\
P^{\Lambda H} &= I - P^{\Sigma} && (27)
\end{aligned}
$$

where $I$ is the identity matrix and the superscript $+$ denotes the Moore-Penrose pseudo inverse. The range of the $P^Σ$ projector corresponds to the RWG non-solenoidal subspace, so that by complementarity, $P^{ΛH}$ projects on the entire solenoidal subspace, which includes the harmonic functions. Additionally, we will use the following dual quasi-Helmholtz projectors, defined as [22]

$$ 
\begin{aligned}
P^{\Lambda} &= \Lambda (\Lambda^{\mathrm{T}} \Lambda)^{+} \Lambda^{\mathrm{T}} && (28) \\
P^{\Sigma H} &= I - P^{\Lambda} && (29)
\end{aligned}
$$

These projectors are the counterparts of the previous ones for the Buffa-Christiansen (BC) basis functions which are defined on the dual barycentric mesh [30]. In this case $P^Λ$ projects on the non-solenoidal subspace of the dual mesh, meaning that by complementarity, the harmonic subspace is implicitly contained in $P^{\Sigma H}$ which projects on the solenoidal subspace. These dual projectors also require the use of a mixed RWG/BC Gram matrix with elements defined as

$$ (G)_{mn} = \langle \hat{n} \times f_m^{\mathrm{RWG}}, f_n^{\mathrm{BC}} \rangle. \quad (30) $$

The inverse of this Gram matrix bridges the operators discretized with primal RWG functions to the projectors acting on the dual BC functions. The reader should note that the above Gram matrix is the only instance in this work where the barycentric refinement is required. For this reason, the method presented here will not suffer from the computational burden of handling dense matrices on the barycentric refinement since this Gram matrix can be obtained analytically without expensive computations [22]. Moreover the projectors can be built with a quasi-linear complexity [31], making them perfectly compatible with fast solvers.

### III. STUDY OF THE PMCHWT IN THE EDDY CURRENT REGIME

#### A. LOW FREQUENCY ASYMPTOTIC ANALYSIS

In this section we present an in depth analysis for the behavior of the PMCHWT matrix in the context of eddy current modeling at low frequencies, as characterized by (6). The reader, however, should note that the regime we consider, $\omega \to 0$ with $\sigma$ constant, is not the only possible regime for low frequency lossy materials. An alternative could be the one in which $\omega \to 0$ with $\sigma\omega$ equals to a constant (i.e., constant skin depth). The reader interested in this regime should refer to the detailed analysis of Bonnet and Demaldent in [14]. In this work we have chosen to limit our solution strategies to the above regime for the sake of brevity and its particular relevancy for applications [5]. The use of quasi-Helmholtz projectors we will propose here, however, can also be extended to other low frequency regimes.

To investigate the conditioning behavior of the system matrix, the latter is decomposed using the Loop-Star transformation matrix $A$ as

$$ Z_{\Lambda H \Sigma} = \begin{pmatrix} A^{\mathrm{T}} & A^{\mathrm{T}} \end{pmatrix} \begin{pmatrix} T_{\mathrm{upper}} & -K \\ K & T_{\mathrm{lower}} \end{pmatrix} \begin{pmatrix} A \\ A \end{pmatrix}. \quad (31) $$

This procedure is only done here for the sake of the analysis and it should not be employed in practice due to a discretization related ill-conditioning introduced by the Loop and Star matrices [31]. Applying this technique effectively highlights the asymptotic scalings for the Loop, harmonic and Star components of each operator block to better understand the origin of the instability of the formulation in eddy current conditions.

Let us start our analysis from the $K$ block (13), for which the key property to determine its asymptotic scalings is that the static part of the operator actually cancels out when solenoidal functions are used simultaneously as testing and source functions, unless both are harmonic (partial cancellations occur in this last case) [32], [33]. In other words,

$$ \Lambda^{\mathrm{T}} K_{0} \Lambda = H^{\mathrm{T}} K_{0} \Lambda = \Lambda^{\mathrm{T}} K_{0} H = 0. \quad (32) $$

The reader familiar with these relationships should keep in mind that in our definition the matrix $K_0$ results from a testing of the MFIO with curl-conforming functions (rotated RWGs, see definition (16)). These properties warrant a closer inspection of the MFIO kernel. The gradient of the Green's function (5) can be expanded using the Taylor series as

$$ 
\begin{aligned}
\nabla G_k(\mathbf{r}, \mathbf{r}') & \approx \frac{1}{4\pi} \nabla \left( \frac{1}{R} - jk - \frac{k^2 R}{2} + \dots \right) \\
& \approx \frac{1}{4\pi} \left( \nabla \left( \frac{1}{R} \right) - \frac{k^2}{2} \nabla R + \dots \right).
\end{aligned}
\quad (33) $$

In the general case, the first term of the expansion constitutes the static part of the operator and dominates as $O(1)$ (here and in the following we will always assume and omit “for $\omega \to 0$”). However, every time one of the cases in (32) occurs, the $O(1)$ term will be canceled and the leading term
becomes the one containing the $O(k^2)$ factor. For the exterior and interior wave numbers, we have

$$O(k_0^2) = O(\omega^2 \mu_0 \epsilon_0) = O(\omega^2) \quad (34)$$

$$O(k_1^2) = O(-j\omega\sigma_1\mu_1) = O(\omega). \quad (35)$$

Using the definition of $K$ in (13), when the static part is canceled, the remaining dynamic part scales as $O(k_0^2) + O(k_1^2) = O(\omega^2) + O(\omega) = O(\omega)$.

For the $T$ blocks, the expansions for $T_{\text{upper}}$ and $T_{\text{lower}}$ in (14) and (15) are used to emphasize the different behaviors of the vector and scalar potentials $T_{A,k}$ and $T_{\Phi,k}$. This becomes clear by rewriting the accompanying factors using the definitions of $k$ and $\eta$ to obtain

$$k_0\eta_0 = \omega\mu_0 = O(\omega) \quad (36)$$

$$k_1\eta_1 = \omega\mu_1 = O(\omega) \quad (37)$$

$$k_0/\eta_0 = \omega\epsilon_0 = O(\omega) \quad (38)$$

$$k_1/\eta_1 = \omega\epsilon_1 = \omega\epsilon_0\epsilon'_r - j\sigma_1 = O(1) \quad (39)$$

where the frequency term in (39) is dominated by the conductivity term, differing from the case of regular dielectrics [23]. From the expressions (14) and (15), the electric operators blocks can be shown to scale as

$$T_{A,\text{upper}} = O(\omega) + O(\omega) = O(\omega) \quad (40)$$

$$T_{\Phi,\text{upper}} = O(\omega^{-1}) + O(1) = O(\omega^{-1}) \quad (41)$$

$$T_{A,\text{lower}} = O(\omega) + O(1) = O(1) \quad (42)$$

$$T_{\Phi,\text{lower}} = O(\omega^{-1}) + O(\omega^{-1}) = O(\omega^{-1}). \quad (43)$$

The inverse scalings between $T_{A,k}$ and $T_{\Phi,k}$ in the upper and lower blocks become evident and in particular, both $T_{\Phi}$ blocks constitute the main reason behind the matrix conditioning breakdown, since they diverge at the limit. Recalling that [22]

$$\Lambda^T T_{\Phi,k} = 0, \quad T_{\Phi,k} \Lambda = 0, \quad H^T T_{\Phi,k} = 0, \quad T_{\Phi,k} H = 0 \quad (44)$$

the Loop-Star decomposed PMCHWT matrix scalings can be written in block form as

$$Z_{\Lambda H \Sigma} = \begin{array}{c} \Lambda \quad H \quad \Sigma \\ \Lambda \quad j\omega\mu_0 \quad j\omega\mu_0 \quad j\omega\mu_0 \\ H \quad j\omega\mu_0 \quad j\omega\mu_0 \quad j\omega\mu_0 \\ \Sigma \quad j\omega\mu_0 \quad j\omega\mu_0 \quad j(\omega\epsilon_0)^{-1} \\ \Lambda \quad j\omega \quad j\omega \quad 1 \\ H \quad j\omega \quad 1 \quad 1 \\ \Sigma \quad 1 \quad 1 \quad 1 \end{array} \begin{array}{c} \Lambda \quad H \quad \Sigma \\ j\omega \quad j\omega \quad 1 \\ j\omega \quad 1 \quad 1 \\ \sigma \quad \sigma \quad \sigma \\ \sigma \quad \sigma \quad \sigma \\ \sigma \quad \sigma \quad j(\omega\mu_0)^{-1} \end{array} \quad (45)$$

The pathological behavior of the matrix is caused by the presence of scalings of different orders in $\omega$. Specifically, the application of the Gershgorin circle theorem shows that

the singular value spectrum of the matrix has a branch which asymptotically diverges as $O(\omega^{-1})$ and at least another which is centered around zero, resulting in a clear degeneration of the matrix condition number.

## B. LOSS OF SOLUTION ACCURACY

While it is necessary to address the conditioning issue of the PMCHWT formulation, it is not sufficient to guarantee an accurate solution. The loss of significant digits that occurs at low frequencies needs to be taken into careful consideration to prevent the corruption of the desired quantities of interest. The rescaling process based on a quasi-Helmholtz decomposition, whether using transformation matrices or projectors, does not solely serve to cure conditioning breakdowns, but also to preserve the needed components that would normally be lost in finite precision arithmetic. Nevertheless, it is not always possible to preserve every single component of the solution, which is why an additional study is necessary to determine which components actually affect the accuracy of the quantities of interest. This stability is targeted for several classes of excitations, for instance the plane wave, inductive, and capacitive excitations [19], [22]. The scattered far field is usually the desired quantity in the case of the plane wave, and thus should be preserved and computed accurately, while circuit simulations with inductive and capacitive excitations require the volumic current inside the conductor, derived from the electric near field.

First, the scalings of the quasi-Helmholtz decomposition of the excitations need to be obtained. For the plane wave, both electric and magnetic fields behave identically [23], whereas they must be examined separately for the inductive and capacitive excitations. These are frequently modeled with a voltage delta gap, for which a magnetic frill around the feed point is a more realistic equivalent [34] that creates the following fields [35]

$$E^i(\mathbf{r}) = - \int_{\text{frill}} \nabla G_k(\mathbf{r}, \mathbf{r}') \times \mathbf{m}(\mathbf{r}') \, d\mathbf{r}' \quad (46)$$

$$H^i(\mathbf{r}) = \frac{1}{j\omega\eta_0^2} \nabla\nabla \cdot \int_{\text{frill}} G_k(\mathbf{r}, \mathbf{r}') \mathbf{m}(\mathbf{r}') \, d\mathbf{r}' + \frac{k_0^2}{j\omega\eta_0^2} \int_{\text{frill}} G_k(\mathbf{r}, \mathbf{r}') \mathbf{m}(\mathbf{r}') \, d\mathbf{r}' \quad (47)$$

where $\mathbf{m}$ is the given magnetic current forming the frill. These expressions show that the electric and magnetic fields are determined, respectively, by the MFIO and the EFIO, for which the reasoning to derive the scalings has already been described in Section III-A. On a multiply connected structure, the magnetic frill corresponds to a poloidal loop around the feed point. Depending on whether there exists or not a toroidal loop that passes through the frill, the problem becomes either inductive or capacitive [19], and the behavior for the electric field changes accordingly [33]. Thus, the asymptotic scalings for each type of right hand side are gathered in Table 1 (a), separated into real and imaginary parts.
the dominant terms of each block of x leads to the asymptotic scalings of the electric and magnetic surface current densities, presented in Table 1 (b). As an example, the scaling of the loop part of the electric current density induced by a plane wave excitation is obtained by multiplying the first line of (48) and (49) by the corresponding scalings in Table 1 (a):

$$ \begin{pmatrix} \omega^{-1/2} + j\omega^{-1} \\ 1 + j\omega^{1/2} \\ \omega^{3/2} + j\omega \\ 1 + j\omega^{1/2} \\ 1 + j\omega^{1/2} \\ 1 + j\omega^{1/2} \end{pmatrix}^T \cdot \begin{pmatrix} \omega^2 + j\omega \\ \omega^2 + j\omega \\ 1 + j\omega \\ \omega^2 + j\omega \\ \omega^2 + j\omega \\ 1 + j\omega \end{pmatrix} = O(1 + j\omega^{1/2}). \quad (50) $$

The distinction between real and imaginary parts of the different elements is made to determine precisely which are the dominant components, listed in the last column of Table 1, that should be preserved to correctly retrieve the corresponding quantities.

The electric far and near fields computed from the electric and magnetic surface current densities can now be derived. The far field is obtained from the following expression [37]

$$ E(\mathbf{r}) \approx -\frac{j\omega\mu_0}{4\pi} \frac{e^{-j\kappa_0 r}}{r} \int_{\Gamma} e^{j\kappa_0 r'} \hat{\mathbf{r}}_s(\mathbf{r}') d\mathbf{r}' + \hat{\mathbf{r}} \times \eta_0 \left( \frac{j\omega\epsilon_0}{4\pi} \frac{e^{-j\kappa_0 r}}{r} \int_{\Gamma} e^{j\kappa_0 r'} \hat{\mathbf{m}}_s(\mathbf{r}') d\mathbf{r}' \right). \quad (51) $$

Note that the solenoidal parts (loop and harmonic) of the currents cancel the static term of the exponential inside the integrals above, which is why the solenoidal scalings of the physical currents are multiplied by jω compared to the non-solenoidal part.

The electric near field inside the conductor is instead calculated with the EFIO and MFIO [37]

$$ E(\mathbf{r}) = \frac{\mu_1}{\sigma_1} \nabla \nabla \cdot \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}') j_s(\mathbf{r}') d\mathbf{r}' \\ - j\omega\mu_1^2 \int_{\Gamma} G_k(\mathbf{r}, \mathbf{r}') j_s(\mathbf{r}') d\mathbf{r}' \\ - \int_{\Gamma} \nabla G_k(\mathbf{r}, \mathbf{r}') \times \mathbf{m}_s(\mathbf{r}') d\mathbf{r}'. \quad (52) $$

The contribution from the solenoidal part of $j_s$ induces a complete cancellation of the first term and a partial cancellation in the second term similarly to the far field above. Performing a Taylor series expansion of the Green function in (51) and (52) when $\omega \to 0$ and keeping the dominant real and dominant imaginary terms will yield the low frequency behavior of the near and far field scattering operators. To obtain the correct behavior, however, the aforementioned cancellations should be explicitly enforced in the expansion when considering the solenoidal components (as is customarily done at very low frequency [24]). Once the six scalings per operator have been obtained, they can be multiplied with those of the physical solutions (Table 1 (b)) to derive the far and near fields asymptotic scalings that are presented in Table 1 (c) and (d).

To illustrate the relevancy of this analysis, we take the example of the plane wave excitation. Table 1 (c) indicates that $\Re(j_\Lambda)$ and $\Im(j_\Sigma)$ are required to compute the far field from the electric current density, while Table 1 (b) shows that the dominant component of the electric part of the solution is $\Re(j_\Lambda)$ only. Therefore, the electric scattered field cannot be computed accurately at low frequencies since $\Im(j_\Sigma)$ is lost from the electric current density when no specific cure is employed. Consequently, the preconditioned equation must produce a solution with leading components that include all dominant parts of the quantities of interest, so that they can be calculated with satisfying accuracy.

## IV. STABILIZATION SCHEME BASED ON QUASI-HELMHOLTZ PROJECTORS

Based on the analysis of the previous section, the requirements to achieve a stable formulation consist in preconditioning the system matrix as well as preventing undesired loss of significant digits in the solution. For this purpose, rescaling coefficients are selected through a Loop-Star decomposition study, by defining diagonal matrices $L = \text{diag}(a_L, b_L, c_L, d_L, e_L, f_L)$ and $R = \text{diag}(a_R, b_R, c_R, d_R, e_R, f_R)$, applied as follows

$$ LZ_{\Lambda H \Sigma} R = O \begin{pmatrix} a_L a_R \omega \mu_0 & a_L b_R \omega \mu_0 & \cdots & a_L f_R \\ b_L a_R \omega \mu_0 & \cdots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ f_L a_R & \cdots & \cdots & f_L f_R (\omega \mu_0)^{-1} \end{pmatrix}. \quad (53) $$

The first condition to fulfill is the prevention of a matrix breakdown. For that, the coefficients must be chosen so as to remove diverging elements and without introducing new ones nor creating null spaces caused by vanishing rows or columns.

As explained in Section III-B, the coefficients also need to account for the loss of accuracy in the solution. When using the Loop-Star method, the base equation $zx = b$ becomes $LZ_{\Lambda H \Sigma}Ry = L\text{diag}(A^T, A^T)b$, with $x = \text{diag}(A, A^T)Ry$. The asymptotic scalings of the rescaled solution $y$ obtained with the preconditioned equation are determined by the right rescaling matrix $R$ as

$$ y = R^{-1} \text{diag}(A^{-1}, A^{-1}) x = \begin{pmatrix} O(a_R^{-1} j_\Lambda) \\ O(b_R^{-1} j_H) \\ O(c_R^{-1} j_\Sigma) \\ O(d_R^{-1} m_\Lambda) \\ O(e_R^{-1} m_H) \\ O(f_R^{-1} m_\Sigma) \end{pmatrix}. \quad (54) $$

For each type of excitation, the dominant components of $y$ must comprise at least the elements necessary to accurately compute the far (for the plane wave) or near field (for inductive and capacitive problems), as indicated by the last
columns of Table 1 (c) and (d). In other words, the choice of rescaling coefficients must simultaneously cure the matrix breakdown shown in (45) and allow an adequate retrieval of the solution.

We first start with eliminating the diverging blocks of the matrix that scale as $O(\omega^{-1})$ by setting $c_L = c_R = (\omega\epsilon_0)^{1/2}$ and $f_L = f_R = (\omega\mu_0)^{1/2}$. Doing this, however, provokes the loss of the Loop and harmonic parts of the magnetic current density for the plane wave excitation ($m_\Lambda$ and $m_H$). Therefore, we choose to adjust them by selecting $d_R = e_R = (\omega/\sigma_1)^{1/2}$. Omitting the factors $\mu_0$, $\epsilon_0$, and $\sigma_1$ for increased readability, the system matrix now scales as

$$ O \begin{array}{c} \Lambda \\ H \\ \Sigma \\ \Lambda \\ H \\ \Sigma \end{array} \begin{array}{c} \begin{array}{ccc} \Lambda & H & \Sigma \end{array} \begin{array}{ccc} \omega & \omega & \omega^{3/2} \\ \omega & \omega & \omega^{3/2} \\ \omega^{3/2} & \omega^{3/2} & 1 \end{array} \begin{array}{ccc} \omega^{3/2} & \omega^{3/2} & \omega^{1/2} \\ \omega^{3/2} & \omega^{1/2} & \omega^{1/2} \end{array} \begin{array}{ccc} \omega & \omega & \omega \\ \omega & \omega & \omega \\ \omega & \omega & \omega \end{array} \end{array} \begin{array}{c} \begin{array}{ccc} \Lambda & H & \Sigma \end{array} \begin{array}{ccc} \omega & \omega & \omega^{1/2} \\ \omega & 1 & \omega^{1/2} \\ \omega^{1/2} & \omega^{1/2} & \omega \end{array} \end{array} . \quad (55) $$

Since several rows and columns vanish as the frequency decreases, the following step is to fix rows 1, 2 and 4, choosing $a_L = b_L = (\omega\mu_0)^{-1/2}$ and $d_L = (\omega\sigma_1)^{-1/2}$, which results in

$$ O \begin{array}{c} \Lambda \\ H \\ \Sigma \\ \Lambda \\ H \\ \Sigma \end{array} \begin{array}{c} \begin{array}{ccc} \Lambda & H & \Sigma \end{array} \begin{array}{ccc} \omega^{1/2} & \omega^{1/2} & \omega \\ \omega^{1/2} & \omega^{1/2} & \omega \\ \omega^{3/2} & \omega^{3/2} & 1 \end{array} \begin{array}{ccc} \omega & \omega & 1 \\ \omega & 1 & 1 \\ \omega & \omega & \omega \end{array} \end{array} \begin{array}{c} \begin{array}{ccc} \Lambda & H & \Sigma \end{array} \begin{array}{ccc} \omega^{1/2} & \omega^{1/2} & 1 \\ \omega & 1 & \omega^{1/2} \\ \omega^{1/2} & \omega^{1/2} & \omega \end{array} \end{array} . \quad (56) $$

We then select $a_R = (\omega\mu_0)^{-1/2}$ to adjust the vanishing first column. At this point, all requirements are fulfilled and two coefficients are left to choose. Both pertain to the harmonic functions, which in the context of quasi-Helmholtz projectors, should always share the same rescaling coefficient as either the adjacent Loop or Star component. Therefore, two possibilities arise that do not compromise the already stable matrix,

$$ \begin{cases} b_R = (\omega\mu_0)^{-1/2} \\ e_L = (\omega\mu_0)^{1/2} \end{cases} \begin{cases} b_R = (\omega\epsilon_0)^{1/2} \\ e_L = (\omega\sigma_1)^{-1/2} \end{cases} \quad (57) $$

Choosing the second set results in the second column of the matrix exhibiting a single non vanishing scaling which is unfortunately excessively small (coefficient $(5, 2) = O((\epsilon_0/\sigma_1)^{1/2})$, thus worsening the condition number. Additionally, the first set allows the retrieval of the surface current densities for the inductive excitation as well. Consequently, after retaining the first pair of coefficients, the complete set of left and right coefficients reads

$$ \begin{pmatrix} a_L \\ b_L \\ c_L \\ d_L \\ e_L \\ f_L \end{pmatrix} = \begin{pmatrix} (\omega\mu_0)^{-1/2} \\ (\omega\mu_0)^{-1/2} \\ (\omega\epsilon_0)^{1/2} \\ (\omega\sigma_1)^{-1/2} \\ (\omega\mu_0)^{1/2} \\ (\omega\mu_0)^{1/2} \end{pmatrix}, \quad \begin{pmatrix} a_R \\ b_R \\ c_R \\ d_R \\ e_R \\ f_R \end{pmatrix} = \begin{pmatrix} (\omega\mu_0)^{-1/2} \\ (\omega\mu_0)^{-1/2} \\ (\omega\epsilon_0)^{1/2} \\ (\omega\sigma_1)^{1/2} \\ (\omega\mu_0)^{1/2} \\ (\omega\mu_0)^{1/2} \end{pmatrix} \quad (58) $$

which gives the preconditioned system matrix in (59), shown at the bottom of the page. Using the Gershgorin circle theorem shows that it is effectively freed from any conditioning breakdown. Moreover, Table 1 (f) contains the scalings for the rescaled solution and, as expected, the last column of the table includes all dominant elements needed for the far and near fields for each type of excitation, confirming the validity of the selected rescaling coefficients. In addition, the current densities are also correctly retrieved for the plane wave and inductive excitations. The reader should notice that, as mentioned above in the coefficient selection process, each coefficient pertaining to the harmonic subspace in (58) has been chosen to always be shared by either the Loop or Star adjacent coefficient. Depending on whether the harmonic coefficient matches with the former or the latter, the corresponding projector will be $P^{\Lambda H}$ or $P^{\Sigma H}$, respectively. Therefore, the rescaling of the projectors, which is compliant with our analysis above, reads

$$ M_1 = (\omega\mu_0)^{-1/2} P^{\Lambda H} + (\omega\epsilon_0)^{1/2} P^{\Sigma} \quad (60) $$

$$ M_2 = (\omega\sigma_1)^{-1/2} P^{\Lambda} + (\omega\mu_0)^{1/2} P^{\Sigma H} \quad (61) $$

for the left side, and

$$ M_3 = (\omega\mu_0)^{-1/2} P^{\Lambda H} + (\omega\epsilon_0)^{1/2} P^{\Sigma} \quad (62) $$

$$ M_4 = (\omega/\sigma_1)^{1/2} P^{\Lambda H} + (\omega\mu_0)^{1/2} P^{\Sigma} \quad (63) $$

for the right side. As described in Section II, the inverse of the Gram matrix (30) is needed to apply the dual projector,

$$ LZ_{\Lambda H \Sigma} R = O \begin{array}{c} \Lambda \\ H \\ \Sigma \\ \Lambda \\ H \\ \Sigma \end{array} \begin{array}{c} \begin{array}{ccc} \Lambda & H & \Sigma \end{array} \begin{array}{ccc} 1 & 1 & \omega(\mu_0\epsilon_0)^{1/2} \\ 1 & 1 & \omega(\mu_0\epsilon_0)^{1/2} \\ \omega(\mu_0\epsilon_0)^{1/2} & \omega(\mu_0\epsilon_0)^{1/2} & 1 \end{array} \begin{array}{ccc} \omega(\mu_0\sigma_1)^{-1/2} & \omega(\mu_0\sigma_1)^{-1/2} & 1 \\ \omega(\mu_0\sigma_1)^{-1/2} & (\mu_0\sigma_1)^{-1/2} & 1 \\ (\mu_0\sigma_1)^{-1/2} & (\mu_0\sigma_1)^{-1/2} & (\epsilon_0/\sigma_1)^{1/2} \end{array} \end{array} \begin{array}{c} \begin{array}{ccc} \Lambda & H & \Sigma \end{array} \begin{array}{ccc} 1 & 1 & \omega(\mu_0\epsilon_0)^{1/2} \\ 1 & 1 & \omega(\mu_0\epsilon_0)^{1/2} \\ \omega(\mu_0\epsilon_0)^{1/2} & \omega(\mu_0\epsilon_0)^{1/2} & 1 \end{array} \end{array} \quad (59) $$
which results in the final formulation

$$ \begin{pmatrix} M_1 & \mathbf{0} \\ \mathbf{0} & \mathbb{M}_2 \mathbf{G}^{-1} \end{pmatrix} Z \begin{pmatrix} M_3 & \mathbf{0} \\ \mathbf{0} & M_4 \end{pmatrix} y = \begin{pmatrix} M_1 & \mathbf{0} \\ \mathbf{0} & \mathbb{M}_2 \mathbf{G}^{-1} \end{pmatrix} b \quad (64) $$

where y is the rescaled solution, from which we derive the current density solution

$$ x = \begin{pmatrix} M_3 & \mathbf{0} \\ \mathbf{0} & M_4 \end{pmatrix} y. \quad (65) $$

## V. IMPLEMENTATION DETAILS

When dealing with quasi-Helmholtz decompositions and related projectors, careful attention is also required on the implementation side to ensure an accurate computation of the operator matrices and right hand sides. Specifically, several terms fail to cancel out exactly in limited precision, leaving residual errors which can be amplified at very low frequencies. Such cancellations should therefore be implemented manually to guarantee the stability and correctness of the formulation.

In particular, the relations $P^{\Lambda H}T_{\Phi,k} = \mathbf{0}$ and $T_{\Phi,k}P^{\Lambda H} = \mathbf{0}$ must be explicitly enforced, as well as $P^{\Lambda}G^{-1}T_{\Phi,k} = \mathbf{0}$ [22], which results for the diagonal blocks in

$$ M_1 T_{\text{upper}} M_3 = M_1 T_{\text{A,upper}} M_3 + \omega \epsilon_0 P^{\Sigma} T_{\Phi,\text{upper}} P^{\Sigma} \quad (66) $$

$$ M_2 T_{\text{lower}} M_4 = M_2 T_{\text{A,lower}} M_4 + \omega \mu_0 P^{\Sigma H} G^{-1} T_{\Phi,\text{lower}} P^{\Sigma} \quad (67) $$

Likewise, the static component of $K_k$ needs to be canceled out since $P^{\Lambda}G^{-1}K_0P^{\Lambda H} = \mathbf{0}$ [23], which leads to

$$ \begin{aligned} M_2 K M_3 = M_2 K_d M_3 &+ \omega (\sigma_1 \epsilon_0)^{1/2} P^{\Lambda} G^{-1} K_0 P^{\Sigma} \\ &+ \omega (\mu_0 \epsilon_0)^{1/2} P^{\Sigma H} G^{-1} K_0 P^{\Sigma} \\ &+ P^{\Sigma H} G^{-1} K_0 P^{\Lambda H} \end{aligned} \quad (68) $$

for the lower off-diagonal block of the system matrix, where $K_d$ is the dynamic part of the discretized operator. In (68), $K_d$ should be computed by explicitly omitting the static term of the kernel (33), instead of calculating $K_k - K_0$, to avoid numerical cancellation. Note that in the case of the upper off-diagonal block, no cancellation is enforced since $H^T K_0 H \neq \mathbf{0}$, and so $P^{\Lambda H}K_0P^{\Lambda H} \neq \mathbf{0}$.

The plane wave right hand side must also be computed with the kernel extraction technique, by removing the static term of the exponential in the integral $(e^{-jkk \cdot r} - 1)$. In fact, the static part disappears with $P^{\Lambda H}$ and $P^{\Lambda}G^{-1}$ [22], which results in

$$ \begin{pmatrix} M_1 e & 0 \\ \mathbb{M}_2 G^{-1} h & 0 \end{pmatrix} = \begin{pmatrix} (\omega \mu_0)^{-1/2} P^{\Lambda H} & \mathbf{0} \\ \mathbf{0} & (\omega \sigma_1)^{-1/2} P^{\Lambda} G^{-1} \end{pmatrix} \begin{pmatrix} e_{\text{ext}} \\ h_{\text{ext}} \end{pmatrix} + \begin{pmatrix} (\omega \epsilon_0)^{1/2} P^{\Sigma} & \mathbf{0} \\ \mathbf{0} & (\omega \mu_0)^{1/2} P^{\Sigma H} G^{-1} \end{pmatrix} \begin{pmatrix} e \\ h \end{pmatrix} \quad (69) $$

where the subscript ext refers to the extracted kernel right hand side.

The new quasi-Helmholtz projector-based formulation, and its discretization, can be accelerated using fast

solvers [25], [26] to reach a quasi-linear complexity $O(N \log N)$. Indeed, the projector-based preconditioner is applicable in quasi-linear complexity by leveraging multigrid preconditioners [38]. The stabilization and the framework presented here are designed for the low-frequency limit. The scheme however is fully compatible with high frequency simulations for frequencies where a loop-star decomposition should not be used. In this case, it is sufficient to set to 1 all coefficients multiplying the projectors as is customary done in similar frameworks (see [24] and references therein). A smooth and automatic transition between the two regimes is possible and will be the topic of a future communication. At moderately low frequencies for which the PMCHWT can still be used, albeit with a high condition number, the overhead caused by the computation and application of the preconditioner is rapidly offset by the lower number of iterations required for the convergence of the employed iterative solver, as will be illustrated in the numerical results. At lower frequencies, the spread of the condition number would become even more severe and the standard PMCHWT would fail to provide the correct solution.

Finally, the complex wave number $k_1$ (7) present in the Green's function exponential gives rise to a real exponential which decays excessively fast as the skin depth becomes smaller, typically due to a high conductivity and/or high frequency. Therefore, the computation of the integral operators becomes increasingly inaccurate when using schemes such as the Gaussian integration. This is an issue characteristic of lossy conducting media which requires specific treatment such as the one found in [39]. After the fast decay of the exponential has been handled, the near-singularity of the right hand side and matrix self- or near-elements require special treatment as they cannot be obtained accurately using standard Gaussian quadrature. In this work, we have opted for singularity extraction approaches [40], [41] but other schemes such as singularity cancellation [42] can be employed.

## VI. NUMERICAL RESULTS

Several tests have been conducted to validate the stability and accuracy of our scheme. For all experiments of this section, we have used $\epsilon'_r = 1$. The first experiment involves a simply connected sphere of radius 1 m and conductivity $\sigma_1 = 10^3 \text{ Sm}^{-1}$. The geometry is discretized with 1048 triangular elements, which corresponds to 3144 unknowns. The condition numbers computed at different frequencies for the PMCHWT equation, the Loop-Star PMCHWT method and the new formulation are illustrated in Figure 3. It is clear that the original PMCHWT matrix displays a rapid degeneration of its conditioning, while both the Loop-Star PMCHWT and our new method remain stable until very low frequencies. However, the condition number achieved with the formulation we propose is significantly lower than that of the Loop-Star scheme, thanks to the well-conditioned nature of the quasi-Helmholtz projectors.
<figure><img src="image_5.png" /><figcaption>**FIGURE 3.** Sphere of radius 1 m with 1048 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Condition number as a function of the frequency.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 4.** Sphere of radius 1 m with 1048 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Electric current density amplitude on elements of constant longitude given an exciting plane wave of frequency $f = 10^{-40} \text{ Hz}$ along $-\hat{x}$ and with $B_0 = 1 \text{ T}$.</figcaption></figure>

We then compared the amplitude of the electric and magnetic current densities against an eddy current specific formulation [7], obtained with a plane wave excitation and computed at the centroids of a subset of elements, shown in Figure 4 and 5. Analytical solutions are available for the magnetic current density [43]–[45], out of which we have implemented the latter. The Loop-Star PMCHWT method and our formulation are in good agreement with both the analytical and eddy current solutions, whereas the PMCHWT equation fails to return correct results. This shows the capacity of the new strategy to maintain accuracy until arbitrarily low frequencies.

The scattered far field is then verified against the Mie solution in Figure 6. Like before, the PMCHWT method delivers a wrong result while the Loop-Star and new methods both result in a good match, as predicted by the theory.

<figure><img src="image_7.png" /><figcaption>**FIGURE 5.** Sphere of radius 1 m with 1048 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Magnetic current density amplitude on elements of constant longitude given an exciting plane wave of frequency $f = 10^{-40} \text{ Hz}$ along $-\hat{x}$ and with $B_0 = 1 \text{ T}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**FIGURE 6.** Sphere of radius 1 m with 1048 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Radar cross section given an exciting plane wave of frequency $f = 10^{-40} \text{ Hz}$ along $-\hat{z}$ and with $E_0 = 1 \text{ Vm}^{-1}$.</figcaption></figure>

To confirm the applicability of our new method to multiply connected structures, we simulated a torus with circular cross section, of major radius 1.5 m, minor radius 0.5 m and conductivity $\sigma_1 = 10^3 \text{ Sm}^{-1}$. The geometry is discretized with 1620 triangular elements, which leads to 4860 unknowns, and the harmonic subspace contains 2 global loops. The condition numbers are computed as before and are illustrated in Figure 7. The same behavior observed in the case of the sphere occurs for the torus, corroborating the theoretical development in the presence of global loops.

The electric and magnetic current densities (Figure 8 and 9) were validated with an exciting plane wave, demonstrating matching results between our new method, the Loop-Star PMCHWT strategy and an eddy current formulation, while the original PMCHWT gives erroneous results. The scattered field (Figure 10) was computed as well, with a similar outcome.
<figure><img src="image_9.png" /><figcaption>**FIGURE 7.** Torus of major and minor radii 1.5 m and 0.5 m with 1620 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Condition number as a function of the frequency.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 8.** Torus of major and minor radii 1.5 m and 0.5 m with 1620 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Electric current density amplitude on elements around the handle given an exciting plane wave of frequency $f = 10^{-5} \text{ Hz}$ along $-\hat{x}$ and with $B_0 = 1 \text{ T}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 9.** Torus of major and minor radii 1.5 m and 0.5 m with 1620 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Magnetic current density amplitude on elements around the handle given an exciting plane wave of frequency $f = 10^{-5} \text{ Hz}$ along $-\hat{x}$ and with $B_0 = 1 \text{ T}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**FIGURE 10.** Torus of major and minor radii 1.5 m and 0.5 m with 1620 elements and $\sigma_1 = 10^3 \text{ Sm}^{-1}$: Radar cross section given an exciting plane wave of frequency $f = 10^{-40} \text{ Hz}$ along $-\hat{z}$ and with $E_0 = 1 \text{ Vm}^{-1}$.</figcaption></figure>

To demonstrate the full-wave capabilities of the new formulation, i.e., that it is stable and accurate beyond the eddy current regime, we have compared the electric and magnetic (Figure 11 and 12, respectively) current densities obtained with the new formulation, the standard PMCHWT, and an eddy current solver to an analytical solution obtained from vector spherical harmonics (VSH) at 10 MHz. As expected the standard PMCHWT and the new solver provide the correct solution, while the eddy current solver does not. This illustrates the versatility of the new formulation that performs accurately at low and higher frequencies.

The conditioning effect of the different formulations is directly reflected in the convergence rate of the solution obtained via iterative solvers. Figure 13 illustrates the decaying speed of the residual error from the Generalized Minimal Residual (GMRES) iterative method, using a plane wave excitation. Our formulation successfully converged to

the correct solution in 255 iterations, and the Loop-Star PMCHWT required 3410 iterations, whereas the original PMCHWT converged to a completely incorrect solution.

We then injected a voltage of 1 mV in a circular ring of major radius 1.0 m, minor radius 0.2 m and conductivity $\sigma_1 = 10^3 \text{ Sm}^{-1}$ at 50 Hz. The object is discretized with 1750 triangular elements and 5250 unknowns. Using circuit theory, we first obtain a resistance of $R = 0.05 \Omega$ from the conductivity and dimensions of the conductor, which gives a predicted current of $I = V/R = 20 \text{ mA}$ after applying Ohm's law. This corresponds to an average eddy current density of $I/A = 0.159 \text{ Am}^{-2}$, where A is the area of the cross section. The eddy currents obtained with the new formulation (Figure 15) are constant along the axis of the ring and vary linearly on the cross section, with an average amplitude matching the above value, as shown in Figure 14.
<figure><img src="image_13.png" /><figcaption>**FIGURE 11.** Sphere of radius 1 m with 1048 elements and $σ_1 = 1 \text{ mSm}^{-1}$: Electric current density amplitude on elements of constant longitude given an exciting plane wave of frequency $f = 10 \text{ MHz}$ along $-\hat{z}$ and with $E_0 = 1 \text{ Vm}^{-1}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**FIGURE 12.** Sphere of radius 1 m with 1048 elements and $σ_1 = 1 \text{ mSm}^{-1}$: Magnetic current density amplitude on elements of constant longitude given an exciting plane wave of frequency $f = 10 \text{ MHz}$ along $-\hat{z}$ and with $E_0 = 1 \text{ Vm}^{-1}$.</figcaption></figure>

To illustrate that the computational overhead caused by the preconditioning is alleviated by the significantly lower number of iterations required for the new formulation, we have performed a comparative run-time study between the standard and preconditioned PMCHWT. In Table 2 we summarize the setup time, overall computation time, and number of iterations required for the simulation of a homogeneous sphere of conductivity $σ_1 = 10^{-3} \text{ Sm}^{-1}$, discretized with 2792 triangles and excited by a plane wave oscillating at 5 MHz for both the new formulation and the standard PMCHWT. The timings presented have been obtained without a fast matrix vector product algorithm.

Lastly, we simulated a plane wave illuminating the outer shell of a jet engine with an average aluminum conductivity of $2 \times 10^7 \text{ Sm}^{-1}$ to verify the applicability of our scheme to a

<figure><img src="image_15.png" /><figcaption>**FIGURE 13.** Torus of major and minor radii 1.5 m and 0.5 m with 1620 elements and $σ_1 = 10^3 \text{ Sm}^{-1}$: Convergence of the residual error for the Generalized Minimal Residual algorithm with relative tolerance $10^{-4}$ given an exciting plane wave of frequency $f = 10^{-40} \text{ Hz}$ along $-\hat{z}$ and with $E_0 = 1 \text{ Vm}^{-1}$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**FIGURE 14.** Torus of major and minor radii 1.0 m and 0.2 m with 1750 elements and $σ_1 = 10^3 \text{ Sm}^{-1}$: Eddy current along the x axis given a voltage excitation of value $V = 1 \text{ mV}$ at frequency $f = 50 \text{ Hz}$ (the torus is symmetrical about the z axis).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**FIGURE 15.** Torus of major and minor radii 1.0 m and 0.2 m with 1750 elements and $σ_1 = 10^3 \text{ Sm}^{-1}$: Eddy current density given a voltage excitation of value $V = 1 \text{ mV}$ at frequency $f = 50 \text{ Hz}$.</figcaption></figure>

realistic example. The electric current density was computed on the surface of the object, as shown in Figure 16. The same solution was obtained through the GMRES method in 3241 iterations with the new formulation, against 18317 with
<table><caption>TABLE 2. Setup and computation time comparison between the new formulation and the PMCHWT.</caption><thead><tr><td>Formulation</td><td>Setup time (s)</td><td># of iterations</td><td>Total time (s)</td></tr></thead><tbody><tr><td>PMCHWT</td><td>59</td><td>1062</td><td>150</td></tr><tr><td>New Formulation</td><td>68</td><td>276</td><td>98</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>FIGURE 16. Jet engine outer shell of length 0.9 m with 9196 elements and $\sigma_1 = 2 \times 10^7 \text{ Sm}^{-1}$; Electric current density norm given an exciting plane wave of frequency $f = 10^{-40} \text{ Hz}$ along $-\hat{z}$ and with $E_0 = 1 \text{ Vm}^{-1}$.</figcaption></figure>

the Loop-Star PMCHWT, confirming again the improved performance of our scheme.

## VII. CONCLUSION

In this work we presented a novel boundary element method for the simulation of eddy current scenarios which relies on the quasi-Helmholtz projectors. The scheme is free from any approximation of the Maxwell's equations, and contrary to existing eddy current models, it can handle frequencies in the eddy current regime (which standard full-wave solvers cannot do) but also beyond it (which standard eddy current solvers cannot do). Furthermore, the formulation is valid for both simply and multiply connected conductors and has a stable conditioning and accurate solution until arbitrarily low frequencies. Lastly, this scheme is compatible with fast solvers, making it a versatile and attractive scheme. The precision and reliability of the new technique have been confirmed through canonical and realistic examples.

## REFERENCES

[1] S. Rao, D. Wilton, and A. Glisson, "Electromagnetic scattering by surfaces of arbitrary shape," *IEEE Trans. Antennas Propag.*, vol. AP-30, no. 3, pp. 409–418, Mar. 1982.

[2] A. J. Poggio and E. K. Miller, *Integral Equation Solutions of Three-Dimensional Scattering Problems*, MB Assoc., New Delhi, India, 1970.

[3] Y. Chang and R. Harrington, "A surface formulation for characteristic modes of material bodies," *IEEE Trans. Antennas Propag.*, vol. AP-25, no. 6, pp. 789–795, Nov. 1977.

[4] T.-K. Wu and L. L. Tsai, "Scattering from arbitrarily-shaped lossy dielectric bodies of revolution," *Radio Sci.*, vol. 12, no. 5, pp. 709–718, 1977.

[5] E. Kriezis, T. D. Tsiboukis, S. M. Panas, and J. A. Tegopoulos, "Eddy currents: Theory and applications," *Proc. IEEE*, vol. 80, no. 10, pp. 1559–1589, Sep. 1992.

[6] J. García-Martín, J. Gómez-Gil, and E. Vázquez-Sánchez, "Non-destructive techniques based on eddy current testing," *Sensors*, vol. 11, no. 3, pp. 2525–2565, 2011.

[7] W. M. Rucker, R. Hoschek, and K. R. Richter, "Various BEM formulations for calculating eddy currents in terms of field variables," *IEEE Trans. Magn.*, vol. 31, no. 3, pp. 1336–1341, May 1995.

[8] D. Zheng, "Three-dimensional eddy current analysis by the boundary element method," *IEEE Trans. Magn.*, vol. 33, no. 2, pp. 1354–1357, Mar. 1997.

[9] A. Buffa, H. Ammari, and J.-C. Nédélec, "A justification of eddy currents model for the Maxwell equations," *SIAM J. Appl. Math.*, vol. 60, no. 5, pp. 1805–1823, 2000.

[10] H. K. Dirks, "Quasi-stationary fields for microelectronic applications," *Elect. Eng. Archiv fur Elektrotechnik*, vol. 79, no. 2, pp. 145–155, 1996.

[11] P. Dlotko, B. Kapidani, S. Pitassi, and R. Specogna, "Fake conductivity or cohomology: Which to use when solving eddy current problems with H-formulations?" *IEEE Trans. Magn.*, vol. 55, no. 6, pp. 1–4, May 2019.

[12] Z. Ren, "T-/spl OMEGA/formulation for eddy-current problems in multiply connected regions," *IEEE Trans. Magn.*, vol. 38, no. 2, pp. 557–560, Mar. 2002.

[13] Y.-H. Chu and W. C. Chew, "A robust surface-integral-equation formulation for conductive media," *Microw. Opt. Technol. Lett.*, vol. 46, no. 2, pp. 109–114, 2005.

[14] M. Bonnet and E. Demaldent, "The eddy current model as a low-frequency, high-conductivity asymptotic form of the Maxwell transmission problem," *Comput. Math. Appl.*, vol. 77, no. 8, pp. 2145–2161, 2019.

[15] D. Wilton and A. Glisson, "On improving the stability of the electric field integral equation at low frequencies," in *Proc. USNC/URSI Spring Meeting Dig.*, Los Angeles, CA, USA, Jun. 1981, p. 24.

[16] Z. G. Qian and W. C. Chew, "A quantitative study on the low frequency breakdown of EFIE," *Microw. Opt. Technol. Lett.*, vol. 50, no. 5, pp. 1159–1162, 2008.

[17] J. E. O. Guzman et al., "On the hierarchical preconditioning of the PMCHWT integral equation on simply and multiply connected geometries," *IEEE Antennas Wireless Propag. Lett.*, vol. 16, pp. 1044–1047, 2016.

[18] G. Vecchi, "Loop-star decomposition of basis functions in the discretization of the EFIE," *IEEE Trans. Antennas Propag.*, vol. 47, no. 2, pp. 339–346, Feb. 1999.

[19] Z.-G. Qian and W. C. Chew, "Enhanced A-EFIE with perturbation method," *IEEE Trans. Antennas Propag.*, vol. 58, no. 10, pp. 3256–3264, Oct. 2010.

[20] T. Xia et al., "An integral equation modeling of lossy conductors with the enhanced augmented electric field integral equation," *IEEE Trans. Antennas Propag.*, vol. 65, no. 8, pp. 4181–4190, Jun. 2017.

[21] H. H. Gan, T. Xia, Q. I. Dai, Y. Li, and W. C. Chew, "Augmented electric-field integral equation for inhomogeneous media," *IEEE Antennas Wireless Propag. Lett.*, vol. 16, pp. 2967–2970, 2017.

[22] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, "On a well-conditioned electric field integral operator for multiply connected geometries," *IEEE Trans. Antennas Propag.*, vol. 61, no. 4, pp. 2077–2087, Dec. 2013.

[23] Y. Beghein, R. Mitharwal, K. Cools, and F. P. Andriulli, "On a low-frequency and refinement stable PMCHWT integral equation leveraging the quasi-Helmholtz projectors," *IEEE Trans. Antennas Propag.*, vol. 65, no. 10, pp. 5365–5375, Aug. 2017.

[24] A. Merlini, Y. Beghein, K. Cools, E. Michielssen, and F. P. Andriulli, "Magnetic and combined field integral equations based on the quasi-Helmholtz projectors," *IEEE Trans. Antennas Propag.*, vol. 68, no. 5, pp. 3834–3846, May 2020.

[25] R. Coifman, V. Rokhlin, and S. Wandzura, "The fast multipole method for the wave equation: A pedestrian prescription," *IEEE Antennas Propag. Mag.*, vol. 35, no. 3, pp. 7–12, Jul. 1993.

[26] K. Zhao, M. N. Vouvakis, and J.-F. Lee, "The adaptive cross approximation algorithm for accelerated method of moments computations of EMC problems," *IEEE Trans. Electromagn. Compatibility*, vol. 47, no. 4, pp. 763–773, Dec. 2005.

[27] T. L. Chhim, J. E. Ortiz, L. Rahmouni, A. Merlini, and F. P. Andriulli, "A quasi-Helmholtz projector stabilized full wave solver encompassing the eddy current regime," in *Proc. IEEE Int. Conf. Electromagn. Adv. Appl. (ICEAA)*, 2019, pp. 1236–1238.

[28] R. Hiptmair, "Boundary element methods for eddy current computation," in *Boundary Element Analysis* (Lecture Notes in Applied and Computational Mechanics), vol. 29, M. Schanz and O. Steinbach, Eds.
<figure><img src="image_19.png" /></figure>

Berlin, Germany: Springer, 2007, pp. 213–248. [Online]. Available: https://doi.org/10.1007/978-3-540-47533-0_9

[29] K. Cools, F. P. Andriulli, F. Olyslager, and E. Michielssen, “Nullspaces of MFIE and Calderón preconditioned EFIE operators applied to toroidal surfaces,” *IEEE Trans. Antennas Propag.*, vol. 57, no. 10, pp. 3205–3215, May 2009.

[30] A. Buffa and S. Christiansen, “A dual finite element complex on the barycentric refinement,” *Math. Comput.*, vol. 76, no. 260, pp. 1743–1769, 2007.

[31] F. P. Andriulli, “Loop-star and loop-tree decompositions: Analysis and efficient algorithms,” *IEEE Trans. Antennas Propag.*, vol. 60, no. 5, pp. 2347–2356, Mar. 2012.

[32] S. Y. Chen, W. C. Chew, J. M. Song, and J.-S. Zhao, “Analysis of low frequency scattering from penetrable scatterers,” *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 4, pp. 726–735, Apr. 2001.

[33] I. Bogaert, K. Cools, F. Andriulli, and D. De Zutter, “Low frequency scaling of the mixed MFIE for scatterers with a non-simply connected surface,” in *Proc. IEEE Int. Conf. Electromagn. Adv. Appl.*, 2011, pp. 951–954.

[34] T. J. Cui and W. C. Chew, “Accurate model of arbitrary wire antennas in free space, above or inside ground,” *IEEE Trans. Antennas Propag.*, vol. 48, no. 4, pp. 482–493, Apr. 2000.

[35] L. Tsai, “A numerical solution for the near and far fields of an annular ring of magnetic current,” *IEEE Trans. Antennas Propag.*, vol. AP-20, no. 5, pp. 569–576, Sep. 1972.

[36] H. V. Henderson and S. R. Searle, “On deriving the inverse of a sum of matrices,” *SIAM Rev.*, vol. 23, no. 1, pp. 53–60, 1981.

[37] C. A. Balanis, *Advanced Engineering Electromagnetics*. Hoboken, NJ, USA: Wiley, 2012.

[38] A. Napov and Y. Notay, “An algebraic multigrid method with guaranteed convergence rate,” *SIAM J. Sci. Comput.*, vol. 34, no. 2, pp. A1079–A1109, 2012.

[39] Z. G. Qian, W. C. Chew, and R. Suaya, “Generalized impedance boundary condition for conductor modeling in surface integral equation,” *IEEE Trans. Microw. Theory Techn.*, vol. 55, no. 11, pp. 2354–2364, Nov. 2007.

[40] D. R. Wilton, J. Rivero, W. A. Johnson, and F. Vipiana, “Evaluation of static potential integrals on triangular domains,” *IEEE Access*, vol. 8, pp. 99806–99819, 2020.

[41] R. D. Graglia, “On the numerical integration of the linear shape functions times the 3-D green’s function or its gradient on a plane triangle,” *IEEE Trans. Antennas Propag.*, vol. 41, no. 10, pp. 1448–1455, Apr. 1993.

[42] F. Vipiana and D. R. Wilton, “Numerical evaluation via singularity cancellation schemes of near-singular integrals involving the gradient of helmholtz-type potentials,” *IEEE Trans. Antennas Propag.*, vol. 61, no. 3, pp. 1255–1265, Nov. 2012.

[43] W. R. Smythe, “Static and dynamic electricity,” in *Electromagnetics*, 3rd ed. New York, NY, USA: Hemisphere, 1989.

[44] C. Bidinosti, E. Chapple, and M. Hayden, “The sphere in a uniform rf field—Revisited,” *Concepts Magn. Reson. B Magn. Reson. Eng. Educ.*, vol. 31, no. 3, pp. 191–202, 2007.

[45] J. R. Nagel, “Induced eddy currents in simple conductive geometries: Mathematical formalism describes the excitation of electrical eddy currents in a time-varying magnetic field,” *IEEE Antennas Propag. Mag.*, vol. 60, no. 1, pp. 81–88, Dec. 2017.

<figure><img src="image_20.png" /></figure>

**TIFFANY L. CHHIM** received the M.Sc. degree in telecommunications engineering from IMT Atlantique, Brest, France, in 2017. She is currently pursuing the Ph.D. degree in Computational electromagnetics with Politecnico di Torino, Turin, Italy.

Her current research interests include preconditioning techniques and spectral analysis of boundary integral equations at high and low frequencies.

<figure><img src="image_21.png" /></figure>

**ADRIEN MERLINI** (Member, IEEE) received the M.Sc.Eng. degree from the École Nationale Supérieure des Télécommunications de Bretagne, France, in 2015, and the Ph.D. degree from the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique), Brest, France, in 2019.

From 2018 to 2019, he was a visiting Ph.D. student with the Politecnico di Torino, Italy, which he then joined as a Research Associate. Since 2019, he has been an Associate Professor with the Microwave Department, IMT Atlantique. His research interests include preconditioning and acceleration of integral equation solvers for electromagnetic simulations and their application in brain imaging.

Dr. Merlini received the Young Scientist Award at the URSI GASS 2020 meeting. In addition, he has coauthored a conference paper recipient of an honorable mention at the URSI/IEEE-APS 2020. He is a member of IEEE-HKN, the IEEE Antennas and Propagation Society, and the Lab-STICC Laboratory.

<figure><img src="image_22.png" /></figure>

**LYES RAHMOUNI** received the B.Sc. degree in aeronautical engineering and the M.Sc. degree in electronic engineering from the University of Science and Technology Houari Boumediene, Algeria, and the Ph.D. degree in computational electromagnetics from IMT Atlantique, France, in 2016, followed by two years of postdoctoral position at the same grande École. He is currently holder of a Research fellowship with Politecnico di Torino, Italy. His research interests focus on integral equation methods and their fast solutions.

He received the 2020 URSI Young Scientist Awards.

<figure><img src="image_23.png" /></figure>

**JOHN ERICK ORTIZ GUZMAN** received the Electronic Engineering degree from the Universidad Nacional de Colombia, Colombia, in 2013, and the Ph.D. degree from the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique), France, in 2018.

From 2018 to 2019, he was working with the Politecnico di Torino as a Research Associate. Since 2020, he has been a full time Professor with the Electronic Engineering Department, Corporación Universitaria Autonoma de Nariño (AUNAR), Pasto, Colombia. His research interests include modeling of electromagnetic problems in biological tissue with integral equations and their acceleration using preconditioning strategies.
<figure><img src="image_24.png" /></figure>

FRANCESCO P. ANDRIULLI (Senior Member, IEEE) received the Laurea degree in electrical engineering from the Politecnico di Torino, Italy, in 2004, the M.Sc. degree in electrical engineering and computer science from the University of Illinois at Chicago in 2004, and the Ph.D. degree in electrical engineering from the University of Michigan at Ann Arbor in 2008.

From 2008 to 2010, he was a Research Associate with the Politecnico di Torino. From 2010 to 2017, he was an Associate Professor from 2010 to 2014, and then a Full Professor with the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique, formerly ENST Bretagne), Brest, France. Since 2017, he has been a Full Professor with the Politecnico di Torino, Turin, Italy. His research interests are in computational electromagnetics with focus on frequency- and time-domain integral equation solvers, well-conditioned formulations, fast solvers, low-frequency electromagnetic analyses, and modeling techniques for antennas, wireless components, microwave circuits, and biomedical applications with a special focus on brain imaging.

Dr. Andriulli was the recipient of the Best Student Paper Award at the 2007 URSI North American Radio Science Meeting, the First Place Prize of the Student Paper Context of the 2008 IEEE Antennas and Propagation Society International Symposium, and the recipient of the 2009 RMTG Award for Junior Researchers and was awarded Two URSI Young Scientist Awards at the International Symposium on Electromagnetic Theory in 2010 and 2013, where he was also awarded the Second Prize in the Best Paper Contest. He also received the 2015 ICEAA IEEE-APWC Best Paper Award. In addition, he has coauthored with his students and collaborators other three first prize conference papers (eMTS 2016, URSI-DE Meeting 2014, and ICEAA 2009), the second prize conference paper (URSI GASS 2014), the third prize conference paper (IEEE-APS 2018), the three honorable mention conference papers (ICEAA 2011, URSI/IEEE-APS 2013, and URSI/IEEE-APS 2020), and other three finalist conference papers (URSI/IEEE-APS 2012, URSI/IEEE-APS 2007, and URSI/IEEE-APS 2006). Moreover, he received the 2014 IEEE AP-S Donald G. Dudley, Jr., Undergraduate Teaching Award, the Triennium 2014–2016 URSI Issac Koga Gold Medal, and the 2015 L. B. Felsen Award for Excellence in Electrodynamics. He is the Editor in Chief of the *IEEE Antennas and Propagation Magazine*, he serves as a Track Editor for the *IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION*, and as an Associate Editor for *IEEE ACCESS*, *URSI Radio Science Letters*, and *IET Microwaves, Antennas and Propagation*. He is the PI of the ERC Consolidator Grant "321." He is a member of Eta Kappa Nu, Tau Beta Pi, Phi Kappa Phi, and the International Union of Radio Science (URSI).