# Low-Frequency Stabilizations of the PMCHWT Equation for Dielectric and Conductive Media: On a Full-Wave Alternative to Eddy-Current Solvers

Viviana Giunzioni, Alberto Scazzola, Adrien Merlini, Francesco P Andriulli

► To cite this version:

Viviana Giunzioni, Alberto Scazzola, Adrien Merlini, Francesco P Andriulli. Low-Frequency Stabilizations of the PMCHWT Equation for Dielectric and Conductive Media: On a Full-Wave Alternative to Eddy-Current Solvers. IEEE Transactions on Antennas and Propagation, In press, 10.1109/tap.2025.3558603. [hal-05087023](https://hal.science/hal-05087023v1)

HAL Id: hal-05087023

https://hal.science/hal-05087023v1

Submitted on 27 May 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Low-Frequency Stabilizations of the PMCHWT
Equation for Dielectric and Conductive Media: On
a Full-Wave Alternative to Eddy-Current Solvers

Viviana Giunzioni, Member, IEEE, Alberto Sazzola, Member, IEEE, Adrien Merlini, Senior Member, IEEE, and
Francesco P. Andriulli, Fellow, IEEE

*Abstract*—We propose here a novel stabilization strategy for the PMCHWT equation that cures its frequency and conductivity related instabilities and is obtained by leveraging quasi-Helmholtz projectors. The resulting formulation is well-conditioned in the entire low-frequency regime, including the eddy-current one, and can be applied to arbitrarily penetrable materials, ranging from dielectric to conductive ones. In addition, by choosing the rescal- ing coefficients of the quasi-Helmholtz components appropriately, we prevent the typical loss of accuracy occurring at low frequency in the presence of inductive and capacitive type magnetic frill excitations, commonly used in circuit modeling to impose a potential difference. Finally, because it relies on quasi-Helmholtz projectors instead of the standard Loop-Star decomposition, our formulation can be used for modeling multiply connected geometries, without incurring the computational overhead caused by the search for the global loops of the structure, while also being compatible with most fast solvers. The efficacy of the proposed preconditioning scheme when applied to both simply and multiply connected geometries is corroborated by numerical examples.

*Index Terms*—Eddy currents, preconditioning, full-wave, multiply connected, quasi-Helmholtz decomposition.

I. INTRODUCTION

THE electromagnetic modeling of dielectric penetrable objects is crucial in a variety of applications, ranging from antenna design to micro- and nano-electronic applications [1], [2]. The boundary element method (BEM) allows to recast this problem in a set of integral equations [3]. Differently from other popular techniques, its use automatically enforces radiation conditions and only requires the discretization of surfaces separating different materials [3], [4]. These characteristics make the BEM one of the most favorable and widespread tools to numerically solve the electromagnetic scattering problem.

The PMCHWT (Poggio-Miller-Chang-Harrington-Wu-Tsai) equation [5] and the Müller equation [6] are among the most common integral formulations for penetrable bodies, and are

respectively of the first- and second-kind. In spite of its favor-
able stability properties, the Müller equation in its standard,
non-conforming, discretization has been shown to yield less
accurate results than the PMCHWT formulation, especially
in high-contrast scenarios [7]–[9]. Unfortunately, because of
its first-kind nature, the PMCHWT formulation suffers from a
dense-discretization ill-conditioning [10], [11]. In addition, the
PMCHWT equation is also plagued by a severe low-frequency
breakdown [10], [12], coupled to further instabilities related
to the conductivity characterizing the obstacle medium.

To circumvent these issues, *ad-hoc* solvers relying on quasi-static approximations of the Maxwell’s system are usually employed to address the low-frequency regime. In particular, the simulation of eddy currents (i.e., conductive currents generated by a time-varying magnetic field [13], [14]), of great interest for many industrial applications [15], has historically required the use of quasi-static solvers that neglect displacement currents [16]. One of the most severe drawback related to these approximations is the challenge of handling the simulation of multi-scale scenarios which restricts their use in real case, complex apparatuses.

Several preconditioning techniques based on the Calderón
identities have been proposed for the PMCHWT formulation,
relying on the fact that the PMCHWT operator, when properly
discretized, is a valid preconditioner for itself [11], [17], [18].
While immune from the dense-discretization breakdown, these
formulations still suffer from contrast dependent instabilities
[11] and undesired current cancellations at very low frequency,
which dramatically affect the accuracy of the simulation
outcomes [12], [19]. An appropriate rescaling of the quasi-
Helmholtz decomposed solenoidal and quasi-irrotational com-
ponents of the system has been proven effective in addressing
this issue [10]. However, when the quasi-Helmholtz decom-
position is performed by means of the Loop-Star change of
basis the condition number behavior with respect to the mesh
refinement is even worsened [10]. Also, the computationally
costly explicit detection of the global-loop quasi-harmonic
subspace is required. This typically requires the search for the
holes and handles of the geometries or the evaluation of the
nullspace of the Loop-Star transformation matrix [10], which
can degrade the computational complexity of the overall prob-
lem that ideally should equal that of solving well-conditioned
integral equations accelerated with fast solvers [20].

Alternatively, the quasi-Helmholtz components of the sys-
tem can be retrieved by projection, that is, by multiplication

This work was supported by the European Innovation Council (EIC) through the European Union's Horizon Europe research Programme under Grant 101046748 (Project CEREBRO) and by the European Union – Next Generation EU within the PNRR project "Multiscale modeling and Engineering Applications" of the Italian National Center for HPC, Big Data and Quantum Computing (Spoke 6) – PNRR M4C2, Investimento 1.4 - Avviso n. 3138 del 16/12/2021 - CN00000013 National Centre for HPC, Big Data and Quantum Computing (HPC) - CUP E13C22000990001.

Viviana Giunzioni and Francesco P. Andriulli are with the Department of
Electronics and Telecommunications, Politecnico di Torino, 10129 Turin, Italy
(e-mail: francesco.andriulli@polito.it).

Adrien Merlini is with the Microwave Department, IMT Atlantique, 29238 Brest, France (e-mail: adrien.merlini@imt-atlantique.fr).
of the operator matrix with the quasi-Helmholtz projectors, introduced in [20], [21] and proved effective at curing the PMCHWT equation's low frequency breakdown in [19]. Their use brings along two significant advantages with respect to the standard Loop-Star decomposition, which are (I) the quasi-harmonic subspace is included in the solenoidal subspace spanned by the solenoidal projector, so that the global-loops detection is not required, and (II) the use of the quasi-Helmholtz projector does not affect the behavior of the condition number with respect to discretization [10].

Using this tool, we develop here a new preconditioning strategy capable of stabilizing the PMCHWT equation with respect to the physical parameters of the problem in the whole low-frequency regime, including the eddy-current one [13], and, contextually, preventing detrimental current cancellations in presence of both inductive and capacitive type magnetic frill excitations. Built upon the standard PMCHWT equation, our full-wave formulation allows for seamless transitions between low and high frequencies, as well as between dielectric and conductive scatterers, resulting in a versatile and computationally efficient solver, compatible with state-of-the-art acceleration strategies [22], [23].

The paper is organized as follows: after setting the necessary background and notation in Section II, Sections III and IV aim at analyzing the shortcomings of the PMCHWT formulation. In particular, Section III, presents an analysis of the asymptotic scalings of the quasi-Helmholtz components of the system as a function of frequency and of the conductivity of the scattering body. Section IV, is focused on the dominant component analysis of the equivalent currents induced by a voltage generator modelled as a magnetic frill, which is common type of excitation employed in circuit simulation. Then, in Section V, we define a quasi-Helmholtz projectors based preconditioning strategy and determine the multiplicative coefficients used in the rescaling of the quasi-Helmholtz components of the system, allowing for the stabilization of the PMCHWT equation with respect to both the electrical length and the conductivity of the scatterer and for the correct recovery of the currents dominant components. The numerical results in Section VI aim at illustrating the good conditioning properties of the proposed formulation applied to both simply and multiply connected geometries, with respect to both frequency and conductivity and the accuracy of the results. Finally, concluding remarks will be given in Section VII.

## II. BACKGROUND AND NOTATION

Let $\Omega_1 \subset \mathbb{R}^3$ be a closed, bounded conductive region with boundary $\Gamma = \partial\Omega_1$ characterized by the outward pointing normal $\hat{n}$, surrounded by the complementary vacuum region $\Omega_0 = \mathbb{R}^3 \setminus \bar{\Omega}_1$, with permittivity $\epsilon_0$ and permeability $\mu_0$. The domain $\Omega_1$, of characteristic length $L$, is filled with a linear, homogeneous, and isotropic conductive medium of permeability $\mu_1 = \mu_r \mu_0$ and complex permittivity $\epsilon_1 = \epsilon'_r \epsilon_0 - j\sigma/\omega$, where $\sigma$ is the conductivity, and $\omega$ is the angular frequency.

Given a wave number $k$ and the three-dimensional free-space Green's function associated to it [4], $G_k(r, r') := e^{-jkR}/(4\pi R)$, with $R := \|r - r'\|$, the electric and magnetic

field integral operators (EFIO and MFIO), $\mathcal{T}_k$ and $\mathcal{K}_k$, are defined as [3], [4]

$$ (\mathcal{T}_k \mathbf{f})(\mathbf{r}) := -jk(\mathcal{T}_{A,k} \mathbf{f})(\mathbf{r}) + \frac{1}{jk}(\mathcal{T}_{\Phi,k} \mathbf{f})(\mathbf{r}), \quad (1) $$

$$ (\mathcal{K}_k \mathbf{f})(\mathbf{r}) := \hat{\mathbf{n}} \times \text{p.v.} \int_{\Gamma} \nabla G_k(\mathbf{r}, \mathbf{r}') \times \mathbf{f}(\mathbf{r}') d\mathbf{r}', \quad (2) $$

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

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

where p.v. indicates Cauchy principal value. When an incident electromagnetic field ($\mathbf{E}^i, \mathbf{H}^i$) impinges on $\Omega_1$, the resulting electric and magnetic surface current densities $j_s := \hat{\mathbf{n}} \times \mathbf{H}$ and $m_s := -\hat{\mathbf{n}} \times \mathbf{E}$ satisfy the PMCHWT equation [5],

$$ \begin{pmatrix} \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{pmatrix} \begin{pmatrix} j_s \\ m_s \end{pmatrix} = \begin{pmatrix} -\hat{\mathbf{n}} \times \mathbf{E}^i \\ -\hat{\mathbf{n}} \times \mathbf{H}^i \end{pmatrix}, \quad (5) $$

where $k_{0,1} := \omega\sqrt{\epsilon_{0,1}\mu_{0,1}}$ is the wave number in $\Omega_{0,1}$ and $\eta_{0,1} := \sqrt{\mu_{0,1}/\epsilon_{0,1}}$ is the characteristic impedance of the exterior or interior medium. In the evaluation of the above square roots, the branch is chosen in accordance with the definition of propagative and lossy Green's functions, compatibly with the convention $e^{-jkR}$.

To discretize and numerically solve (5), we approximate the surface $\Gamma$ with a mesh of planar triangular elements with average edge length $h$, over which we can define different sets of basis functions. In the following, we will consider Rao-Wilton-Glisson (RWG) basis functions. For their definition and visual representation, we refer to [10] and references therein. By expanding the unknown current densities as linear combinations of RWG functions, as

$$ j_s = \sum_{n=1}^{N_e} j_n f_n, \quad m_s = \sum_{n=1}^{N_e} m_n f_n, \quad (6) $$

where $N_e$ is the number of edges in the mesh, and by testing the resulting equations with curl-conforming rotated RWG functions, following a conforming Petrov-Galerkin discretization procedure, we obtain the linear system of equations

$$ \begin{pmatrix} \mathcal{T}_u & -\mathcal{K} \\ \mathcal{K} & \mathcal{T}_l \end{pmatrix} \begin{pmatrix} j \\ m \end{pmatrix} = \begin{pmatrix} e \\ h \end{pmatrix}, \quad \begin{pmatrix} \mathcal{T}_u & -\mathcal{K} \\ \mathcal{K} & \mathcal{T}_l \end{pmatrix} =: Z \quad (7) $$

where

$$ \mathcal{K} := \mathcal{K}_{k_0} + \mathcal{K}_{k_1}, \quad (8) $$

$$ \mathcal{T}_u := \eta_0 \mathcal{T}_{k_0} + \eta_1 \mathcal{T}_{k_1}, \quad (9) $$

$$ \mathcal{T}_l := \frac{1}{\eta_0} \mathcal{T}_{k_0} + \frac{1}{\eta_1} \mathcal{T}_{k_1}. \quad (10) $$

The elements of the above matrices discretizing the EFIO and MFIO read

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

$$ (\mathcal{K}_k)_{mn} := (\hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{K}_k(\mathbf{f}_n))_{L^2(\Gamma)}, \quad (12) $$

$$ (\mathcal{T}_{A,k})_{mn} := (\hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{T}_{A,k}(\mathbf{f}_n))_{L^2(\Gamma)}, \quad (13) $$

$$ (\mathcal{T}_{\Phi,k})_{mn} := (\hat{\mathbf{n}} \times \mathbf{f}_m, \mathcal{T}_{\Phi,k}(\mathbf{f}_n))_{L^2(\Gamma)}. \quad (14) $$
The right-hand-side (RHS) resulting from the testing procedure is

$$ (e)_m := (\hat{\mathbf{n}} \times \mathbf{f}_m, -\hat{\mathbf{n}} \times \mathbf{E}^i)_{L^2(\Gamma)}, \quad (15) $$

$$ (h)_m := (\hat{\mathbf{n}} \times \mathbf{f}_m, -\hat{\mathbf{n}} \times \mathbf{H}^i)_{L^2(\Gamma)}, \quad (16) $$

while the elements of the unknown vector are simply given by $(j)_n := j_n$ and $(m)_n := m_n$. In the following we will focus on a rescaled version of the equation, that is, on matrix $\bar{Z}$, which reads

$$ \frac{\left( \begin{array}{cc} \mathbf{T}_u / \eta_0 & -\mathbf{K} \\ \mathbf{K} & \eta_0 \mathbf{T}_l \end{array} \right)}{\bar{Z}} \left( \begin{array}{c} \sqrt{\eta_0 j} \\ m / \sqrt{\eta_0} \end{array} \right) = \left( \begin{array}{c} e / \sqrt{\eta_0} \\ \sqrt{\eta_0} h \end{array} \right), \quad (17) $$

where

$$ \mathbf{T}_u / \eta_0 = \mathbf{T}_{k_0} - j k_0 \mu_r \mathbf{T}_{A, k_1} + \frac{1}{j k_0 \epsilon_r} \mathbf{T}_{\Phi, k_1}, \quad (18) $$

$$ \eta_0 \mathbf{T}_l = \mathbf{T}_{k_0} - j k_0 \epsilon_r \mathbf{T}_{A, k_1} + \frac{1}{j k_0 \mu_r} \mathbf{T}_{\Phi, k_1}. \quad (19) $$

This rescaling leads to an effective conductivity stabilization at mid-frequencies, as will be shown in the following, but cannot prevent the dramatic increase of the conditioning towards lower frequencies.

We introduce at this point the mixed Gram matrices $\mathbf{G}$ and $\mathbf{G} = -\mathbf{G}^T$, defined as $(\mathbf{G})_{mn} := (\hat{\mathbf{n}} \times \mathbf{f}_m, g_n)_{L^2(\Gamma)}$ and $(\mathbf{G})_{mn} := (\hat{\mathbf{n}} \times g_m, \mathbf{f}_n)_{L^2(\Gamma)}$, where $g_n$ is the dual with respect to the RWG function defined on the n-th edge. In this work, we employed Buffa-Christianes (BC) functions but other possibilities are feasible (see [10] and references therein for the definition of BC functions and for further details).

In Section III we will apply the Loop-Star decomposition technique to the PMCHWT operator for analyzing its spectral properties, and we denote with $\mathbf{A}$, $\mathbf{H}$, and $\mathbf{\Sigma}$ the transformation matrices from the loop, global loop, and star subspaces to the RWG subspace. For their explicit definition and interpretation, we refer the reader to [10], [21]. We also define the block matrix $\mathbf{A}$ as $\mathbf{A} := (\mathbf{\Lambda} \mathbf{H} \mathbf{\Sigma})$, that can be used to perform the required decomposition. It is worth noting that we are only introducing matrix $\mathbf{H}$ for the purpose of the analysis, but its explicit (potentially costly) evaluation will not be required for the implementation of the proposed formulation. The main building blocks of the actual preconditioning scheme presented in this work are the quasi-Helmholtz projectors, for both primal and dual spaces, defined as

$$ \mathbf{P}^{\mathbf{\Sigma}} := \mathbf{\Sigma} (\mathbf{\Sigma}^{\mathrm{T}} \mathbf{\Sigma})^{+} + \mathbf{\Sigma}^{\mathrm{T}}, \quad \mathbf{P}^{\mathbf{\Lambda H}} := \mathbf{I} - \mathbf{P}^{\mathbf{\Sigma}}, \quad (20) $$

$$ \mathbf{P}^{\mathbf{\Lambda}} := \mathbf{\Lambda}(\mathbf{\Lambda}^{\mathrm{T}}\mathbf{\Lambda})^{+} + \mathbf{\Lambda}^{\mathrm{T}}, \quad \mathbf{P}^{\mathbf{\Sigma H}} := \mathbf{I} - \mathbf{P}^{\mathbf{\Lambda}}, \quad (21) $$

where the subscript + indicates the Moore-Penrose pseudoinverse. In particular, $\mathbf{P}^{\mathbf{\Sigma}}$ and $\mathbf{P}^{\mathbf{\Lambda}}$ are the non-solenoidal projectors for primal and dual functions, built respectively from the star and the loop transformation matrices, given that the dual BC functions represent a rotation of the primal RWG functions; by complementarity, $\mathbf{P}^{\mathbf{\Lambda H}}$ and $\mathbf{P}^{\mathbf{\Sigma H}}$ are the solenoidal projectors for primal and dual functions, that means that they are projectors onto the solenoidal subspace—which includes the quasi-harmonic subspace. Taking advantage of algebraic multigrid preconditioning techniques, these projectors can be applied to a vector in quasi-linear complexity [10],

characteristic which makes them compatible with standard fast solvers.

### III. CONDITIONING ANALYSIS OF THE PMCHWT EQUATION

This section aims at studying some of the spectral properties of the PMCHWT system matrix $\bar{Z}$ in the low-frequency regime and for different conductivity levels $\sigma$ of the body $\Omega$. The analysis is based on the determination of the asymptotic scalings towards low frequency of the quasi-Helmholtz components of the decomposed system, where the decomposition is performed by means of loop and star transformation matrices, for analytical purposes only.

First, after recalling that $L = \mathrm{diam}(\Omega_1)$, we introduce the dimensionless scalar parameters

$$ \chi := k_0 L, \quad (22) $$

$$ \gamma := \sqrt{\omega \epsilon_0 / \sigma}, \quad (23) $$

$$ \xi := \sqrt{2 / \mu_r L / \delta}, \quad (24) $$

where $\delta$ denotes the skin depth parameter, or penetration length of the field inside the scatterer, which is function of the frequency and the conductivity and defined as $\delta := \sqrt{2/(\omega\sigma\mu_1)}$. Given the definitions above, we have $\chi = \gamma\xi$. Before moving to the study of the interior Green's function kernel, strictly dependent on the frequency-conductivity regime considered (the nomenclature of which has been partially borrowed from [2]), we start by analyzing the asymptotic behavior of matrices $\mathbf{T}_{A,k_0}$, $\mathbf{T}_{\Phi,k_0}$, and $\mathbf{K}_{k_0}$ for vanishing frequency, as well as their Loop-Star components. In the low-frequency limit, for $\chi \to 0$, the exterior Green's function and its gradient can be expanded as

$$ G_{k_0}(\mathbf{r}, \mathbf{r}') = \frac{1}{4\pi R} \left[ 1 - j \frac{R}{L} \chi - \frac{R^2}{2L^2} \chi^2 + \mathcal{O}(\chi^3) \right], \quad (25) $$

$$ \nabla G_{k_0}(\mathbf{r}, \mathbf{r}') = \frac{1}{4\pi} \left[ \nabla \left( \frac{1}{R} \right) - \frac{\chi^2}{2L^2} \nabla R + \mathcal{O}(\chi^3) \right]. \quad (26) $$

From these expansions, it can be recognized that, for $\chi \to 0$,

$$ \| \mathbf{T}_{A, k_0} / L \| = \mathcal{O}(1), \quad \| \mathbf{K}_{k_0} \| = \mathcal{O}(1), \quad (27) $$

$$ \| \mathbf{T}_{\Phi, k_0} / L \| = \mathcal{O}(1), \quad \| \mathbf{K}_{\mathrm{ext}, k_0} \| = \mathcal{O}(\chi^2), \quad (28) $$

where $\mathbf{K}_{\mathrm{ext},k} = \mathbf{K}_k - \mathbf{K}_0$. Here and in the following, the notation $f = \mathcal{O}(g)$ indicates that there exists $c > 0$ such that $|f| \le c|g|$ in the indicated limit. Note that in this contribution we consider the low frequency limit and thus the constant $c$ may be function of the mesh refinement parameter $h$, that is, $c = c(h)$. However, this characterization of $c$ in function of $h$ is outside the scope of this work. Equations (27), (28) result in the Loop-Star decompositions

$$ \mathbf{A}^{\mathrm{T}} (\mathbf{T}_{k_0}) \mathbf{A} = \mathcal{O} \left( \begin{pmatrix} \chi & \chi & \chi \\ \chi & \chi & \chi \\ \chi & \chi & \chi^{-1} \end{pmatrix} \right), \quad (29) $$

$$ \mathbf{A}^{\mathrm{T}} (\mathbf{K}_{k_0}) \mathbf{A} = \mathcal{O} \left( \begin{pmatrix} \chi^2 & \chi^2 & 1 \\ \chi^2 & 1 & 1 \\ 1 & 1 & 1 \end{pmatrix} \right), \quad (30) $$
where we have enforced the cancellation of the scalar potential
contribution when tested against or applied to a solenoidal
function [20]

$$
\Lambda^{\mathrm{T}} \mathbf{T}_{\Phi,k} = \mathbf{0}, \quad \mathbf{H}^{\mathrm{T}} \mathbf{T}_{\Phi,k} = \mathbf{0}, \quad \mathbf{T}_{\Phi,k} \Lambda = \mathbf{0}, \quad \mathbf{T}_{\Phi,k} \mathbf{H} = \mathbf{0} \tag{31}
$$

and the cancellation of the static magnetic field operator $\mathbf{K}_0$ [24]

$$
\Lambda^{\mathrm{T}} \mathbf{K}_0 \Lambda = 0, \quad \Lambda^{\mathrm{T}} \mathbf{K}_0 \mathbf{H} = 0, \quad \mathbf{H}^{\mathrm{T}} \mathbf{K}_0 \Lambda = 0. \quad (32)
$$

The decompositions in (30) will be useful to determine the
overall Loop-Star decomposed PMCHWT matrix $Z_{\Lambda H \Sigma}$,

$$
Z_{\Lambda H \Sigma} := \begin{pmatrix} \mathbf{A}^{\mathrm{T}} & \mathbf{0} \\ \mathbf{0} & \mathbf{A}^{\mathrm{T}} \end{pmatrix} \bar{Z} \begin{pmatrix} \mathbf{A} & \mathbf{0} \\ \mathbf{0} & \mathbf{A} \end{pmatrix}. \quad (33)
$$

We proceed now to the definition of the three low-frequency
regimes considered in this work and to the interior kernel
analysis.

A. Analysis in the Quasi-Static Regime

The quasi-static region is characterized by slowly varying
fields, for which the quasi-stationary condition $\chi \ll 1$ holds
true. In addition, the ratio between the angular frequency
and the conductivity should be high enough to make any
conductive current inside the object negligible with respect
to the displacement current induced by the impinging field. In
formulae, the quasi-static region is characterized by

$$
\begin{cases}
\chi \ll 1 \\
\gamma \gg 1
\end{cases}
. \tag{34}
$$

The asymptotic behavior of the blocks of the Loop-Star
decomposed PMCHWT equation, allowing the study of the
stability properties of the formulation, can be deduced from
the low-frequency asymptotic expansion of the interior Green's
function and of its gradient, under the assumption $\gamma \gg 1$.
This corresponds to analyzing the low-frequency limit along
a constant-$\gamma$ line (blue arrow in Figure 1, corresponding
to $\sigma \propto \omega$) with $\gamma \gg 1$. The regime characterized by
$\gamma = \text{const} \gg 1$ while $\chi \to 0$ will be referred to as quasi-
static regime, abbreviated as QSR.

Along this limit both the real and the imaginary part of
the interior wave number go to zero. In particular, writing the
expression of the interior wavenumber in the form

$$
k_1 = k_0 \sqrt{\epsilon'_r \mu_r} \sqrt{1 - j \frac{\sigma}{\omega \epsilon_0 \epsilon'_r}}, \quad (35)
$$

shows that, for $\chi \to 0$, $\Re\{k_1 L\} = \mathcal{O}(\chi)$ and $\Im\{k_1 L\} = \mathcal{O}(\chi)/\gamma^2$. This fact allows us to exploit the small argument Taylor expansions of the Green's function and of its gradient to determine the low-frequency scalings

$$
\| \mathbf{T}_{A,k_1} / L \| = \mathcal{O}(1), \quad \| \mathbf{K}_{k_1} \| = \mathcal{O}(1), \quad (36)
$$

$$
\| \mathbf{T}_{\Phi, k_1} L \| = \mathcal{O}(1), \quad \| \mathbf{K}_{\mathrm{ext}, k_1} \| = \mathcal{O}(\chi^2), \quad (37)
$$

leading to the quasi-Helmholtz decomposition

$$
Z_{\Lambda H \Sigma} = \mathcal{O} \begin{pmatrix} \chi & \chi & \chi & \chi^2 & \chi^2 & 1 \\ \chi & \chi & \chi & \chi^2 & 1 & 1 \\ \chi & \chi & \chi^{-1} & 1 & 1 & 1 \\ \chi^2 & \chi^2 & 1 & \chi & \chi & \chi \\ \chi^2 & 1 & 1 & \chi & \chi & \chi \\ 1 & 1 & 1 & \chi & \chi & \chi^{-1} \end{pmatrix} . \quad (38)
$$

By applying the Gershgorin circle theorem, we recognize
that the condition number of the PMCHWT equation grows
when decreasing the frequency, resulting in the so-called low-
frequency breakdown of the formulation. Moreover, we infer
that, as in the EFIE case [20], the pathological frequency
behavior of the condition number is in part due to the
opposed frequency scaling of the vector and scalar potential
components of the electric operator $T_k$. Numerical evidence
indicates that the condition number of the PMCHWT formu-
lation increases as $\omega^{-2}$ in the low-frequency limit, as for the
EFIE.

B. Analysis in the Non-Skin-Effect-Dominated Eddy-Current Regime

The eddy-current region is a quasi-stationary region ($\chi \ll 1$) in which the typical time-scale is longer than the electric relaxation time constant $\tau = \epsilon_0/\sigma$, i.e. $\omega\tau \ll 1$, equivalent to $\gamma \ll 1$. In addition, the non-skin-effect-dominated subset of the eddy-current region is characterized by a further condition on the skin depth parameter, which should be much larger with respect to the characteristic size of the object, resulting in the condition $\xi \ll 1$ in the regime under study. Therefore, in the non-skin-effect-dominated eddy-current region, the frequency, geometrical and electrical parameters of the scatterer are such that

$$
\begin{cases}
\chi \ll 1 \\
\gamma \ll 1 \\
\xi \ll 1
\end{cases}
. \tag{39}
$$

In the eddy-current region ($\chi \ll 1$ and $\gamma \ll 1$) the interior
wave number, which after some manipulations can be written
as

$$
k_1 = \frac{1-j}{\delta} \sqrt{1 + j \frac{\omega \epsilon_0 \epsilon'_r}{\sigma}}, \quad (40)
$$

can be approximated as $k_1 \approx (1-j)/\delta$. Hence, both the real
and the imaginary parts of $k_1 L$ decrease in magnitude as $\xi$
when decreasing the frequency.

In this region, we want to analyze the low-frequency limit
of the PMCHWT equation by keeping the ratio $\xi/\gamma$ constant
(i.e., constant conductivity), corresponding to the black arrow
in Figure 1. This regime, characterized by $\xi/\gamma = \text{const}$ while
$\chi \to 0$, will be referred to as non-skin-effect-dominated eddy-
current regime, abbreviated as non-SED ECR. Along this limit,
$\Re\{k_1 L\} = \mathcal{O}(\chi/\gamma)$ and $\Im\{k_1 L\} = \mathcal{O}(\chi/\gamma)$ for $\chi \to 0$. Once
again we can take advantage of the small argument expansions
of the Green's function and of its gradient to determine the
scalings

$$
\| \mathbf{T}_{A,k_1} / L \| = \mathcal{O}(1), \quad \| \mathbf{K}_{k_1} \| = \mathcal{O}(1), \quad (41)
$$

$$
\| \mathbf{T}_{\Phi, k_1} L \| = \mathcal{O}(1), \quad \| \mathbf{K}_{\text{ext}, k_1} \| = \mathcal{O}(\chi/\gamma)^2 . \quad (42)
$$
The Loop-Star decomposed PMCHWT equation in the non-SED ECR limit follows as

$$Z_{\Lambda H \Sigma} = \mathcal{O} \left( \begin{pmatrix} \chi & \chi & \chi & \chi^2 \gamma^{-2} & \chi^2 \gamma^{-2} & 1 \\ \chi & \chi & \chi & \chi^2 \gamma^{-2} & 1 & 1 \\ \chi & \chi & \chi^{-1} & 1 & 1 & 1 \\ \chi^2 \gamma^{-2} & \chi^2 \gamma^{-2} & 1 & \chi \gamma^{-2} & \chi \gamma^{-2} & \chi \gamma^{-2} \\ \chi^2 \gamma^{-2} & 1 & 1 & \chi \gamma^{-2} & \chi \gamma^{-2} & \chi \gamma^{-2} \\ 1 & 1 & 1 & \chi \gamma^{-2} & \chi \gamma^{-2} & \chi^{-1} \end{pmatrix} \right) \quad (43)$$

As in the previous case, the application of the Gershgorin circle theorem shows a low-frequency breakdown of the equation, as the condition number grows at least as $\omega^{-1}$ in the low-frequency limit.

Indeed, the matrix eigenvalues lie in three Gershgorin circles with radius $\mathcal{O}(1)$ and centered in values which are $\mathcal{O}(\chi)$, $\mathcal{O}(\chi^{-1})$, and $\mathcal{O}(\chi\gamma^{-2})$. Also notice then, when fixing the frequency, $\chi\gamma^{-2} = \sigma\eta_0 L$ ranges from $\chi$ to $\chi^{-1}$ when the conductivity varies from $\sigma_{\text{low}} := \omega\epsilon_0$ (corresponding to $\gamma = 1$) to $\sigma_{\text{high}} := (\omega L^2\mu_1)^{-1}$ (corresponding to $\xi = 1$).

### C. Analysis in the Skin-Effect-Dominated Eddy-Current Regime

The skin-effect-dominated eddy-current region is characterized by

$$\begin{cases} \chi \ll 1 \\ \gamma \ll 1 \\ \xi \gg 1 \end{cases} . \quad (44)$$

In this case, we aim at studying the asymptotic behaviour of the equation in the low-frequency limit along a constant-$\xi$ line, with $\xi \gg 1$ (red arrow in Figure 1, corresponding to $\sigma \propto \omega^{-1}$). This regime, characterized by $\xi = \text{const} \gg 1$ while $\chi \to 0$, will be referred to as skin-effect-dominated eddy-current regime, abbreviated as SED ECR.

The eddy-current conditions ($\chi \ll 1$ and $\gamma \ll 1$) guarantee the validity of the approximation of the interior wave number as $k_1 \simeq (1-j)/\delta$ mentioned before. The fundamental difference with respect to the other low-frequency regimes analyzed above is that in this case the product $k_1 L$ does not vanish in the low-frequency limit, preventing the use of the Maclaurin series of the Green's function and its gradient for the evaluation of the asymptotic behavior of the operators. However, by noting that the Green's function behavior is dominated by the exponential decay $e^{-R/\delta}$, through some technical steps it can be proved that

$$\| \mathcal{T}_{A,k_1}/L \| = \xi^{-1} \mathcal{O}(1), \quad \| \mathcal{K}_{k_1} \| = \xi^{-1} \mathcal{O}(1), \quad (45)$$

$$\| \mathcal{T}_{\Phi, k_1} L \| = \xi^{-1} \mathcal{O}(1), \quad \| \mathcal{K}_{\text{ext}, k_1} \| = \mathcal{O}(1) \quad (46)$$

for $\chi \to 0$, where $\xi^{-1} = \gamma/\chi$ is constant.

Indeed, as far as the vector and scalar electric potentials are concerned, it can be seen that, in the low-skin-depth limit, only the self and near interactions give rise to a non negligible contribution as the skin depth becomes smaller than the characteristic size of the mesh triangles $h$. Hence, we can grasp an intuition of the above scalings by noting that the Green's function kernel in modulo can be approximated as $|G_k(\mathbf{r}, \mathbf{r}')| \simeq e^{-R/\delta}/(4\pi R)$ and that the integral of such

function over an infinite plane, which is a good estimate of the integral over the triangle $T$ due to the exponential decay of the function for source points inside $T$ and testing points outside $T$, is

$$\int_{\mathbb{R}^2} \frac{e^{-|\mathbf{r}|/\delta}}{4\pi|\mathbf{r}|} d\mathbf{r} = \delta/2, \quad (47)$$

determining the scalings of the norms $\|\mathcal{T}_{A,k_1}\|$ and $\|\mathcal{T}_{\Phi, k_1}\|$ with the skin depth. Also, in the same low-skin-depth-limit, the magnetic field integral operator matrix becomes negligible, as the kernel exactly cancels for testing points lying on the source triangle. Finally, each of $\mathcal{T}_{A,k_1}$, $\mathcal{T}_{\Phi, k_1}$, and $\mathcal{K}_{k_1}$, denoted here with the placeholder $M$, can be written as the sum of a matrix accounting for the interactions between basis functions defined on non overlapping domains, $M_{\text{far}}$, and the remaining $M_{\text{near}}$, that is $M = M_{\text{near}} + M_{\text{far}}$. It can be shown that, in the low-frequency limit, $\|M_{\text{far}}\| = e^{-p\xi}\mathcal{O}(1)$, where $pL > 0$ represents the minimum distance between two non-overlapping domains of basis functions. Hence, $\|M\| \le \|M_{\text{near}}\| + \|M_{\text{far}}\| = \xi^{-1}\mathcal{O}(1) + e^{-p\xi}\mathcal{O}(1)$.

The accurate evaluation of the inner integrals for the definition of the operator matrices is not trivial and may require the implementation of different integration techniques specifically tailored to highly lossy media. In our work, we employed the integration scheme presented in [25] to produce accurate numerical results.

Due to the condition $\gamma \ll 1$ the diagonal blocks of $\bar{Z}$ can be approximated as

$$\mathcal{T}_u / \eta_0 \approx \mathcal{T}_{k_0} - j k_0 \mu_r \mathcal{T}_{A, k_1} + \frac{1}{\sigma \eta_0 L} (\mathcal{T}_{\Phi, k_1} L) \quad (48)$$

$$\mathcal{T}_l \cdot \eta_0 \approx \mathcal{T}_{k_0} - \sigma \eta_0 L (\mathcal{T}_{A, k_1} / L) + \frac{1}{jk_0 \mu_r} \mathcal{T}_{\Phi, k_1}, \quad (49)$$

It holds $\sigma\eta_0 L = \chi\gamma^{-2} \gg \chi^{-1}$ (following from the condition $\xi \gg 1$). As a result, in the SED ECR, the asymptotic behavior of the upper diagonal block is determined by $\mathcal{T}_{k_0}$, while the behavior of the lower diagonal block is determined by the dominating term $\sigma\eta_0 L(\mathcal{T}_{A, k_1}/L)$, for which it holds

$$\| \sigma \eta_0 L (\mathcal{T}_{A, k_1} / L) \| = \mathcal{O}(\chi \gamma^{-2} \xi^{-1}) = \mathcal{O}(\gamma^{-1}). \quad (50)$$

Under these conditions, the Loop-Star decomposition of the PMCHWT equation in the SED ECR is

$$Z_{\Lambda H \Sigma} = \mathcal{O} \left( \begin{pmatrix} \chi & \chi & \chi & 1 & 1 & 1 \\ \chi & \chi & \chi & 1 & 1 & 1 \\ \chi & \chi & \chi^{-1} & 1 & 1 & 1 \\ 1 & 1 & 1 & \gamma^{-1} & \gamma^{-1} & \gamma^{-1} \\ 1 & 1 & 1 & \gamma^{-1} & \gamma^{-1} & \gamma^{-1} \\ 1 & 1 & 1 & \gamma^{-1} & \gamma^{-1} & \gamma^{-1} \end{pmatrix} \right) . \quad (51)$$

Also in this regime the formulation is plagued by the low-frequency breakdown, that is, the condition number of $\bar{Z}$ increases at least as $\omega^{-1}$ in the low-frequency limit. Moreover, a further instability with respect to $\sigma$ can be recognized in the high-conductivity limit, related to the branch of singular values accumulating at $\mathcal{O}(\gamma^{-1})$, with $\gamma^{-1}$ always higher that $\chi^{-1}$ in the SED ECR, resulting in a condition number increase when fixing the frequency and increasing the conductivity of the scatterer.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Graphical representation of the low-frequency limits in the three regions: the blue arrow represents the limit direction in the QSR, the black arrow in the non-SED ECR, and the red arrow in the SED ECR.</figcaption></figure>

Before moving to the definition of a quasi-Helmholtz pro-
jectors based preconditioning strategy capable of curing the
instabilities identified so far, we first propose an asymptotic
analysis of the quasi-Helmholtz components of the solution
of the system, allowing to determine the current contribu-
tions required for the accurate evaluation of the electric and
magnetic fields inside and outside the scatterer in the low-
frequency limit. This preliminary analysis will indeed result
in further requirements to be satisfied in the design of the
preconditioning strategy.

IV. ASYMPTOTIC ANALYSIS OF THE CURRENT COMPONENTS

While a stable and low matrix condition number is required
for the convergence of an iterative solver [26], it is not
sufficient to guarantee the accuracy of the solution [10], [19],
[20]. Indeed, numerical loss of significance at low frequencies
can provoke the corruption of some non-dominant quasi-
Helmholtz current components, which, in turn, can lower the
accuracy of the scattered fields. In this section, we present
a low-frequency asymptotic analysis aimed at determining
which are the current components required for the correct
reconstruction of the electric and magnetic fields inside and
outside the body Ω. This information will constitute a further
constraint in the design of the preconditioning strategy, which
will have to allow for the correct recovery of these required
current components.

The analysis will take advantage of the low-frequency
asymptotic expansion of the PMCHWT equation along the
low-frequency limits previously identified, which are, along
a constant − γ line with γ ≫ 1 for the QSR, along a
constant − ( ξ/γ ) line with ξ ≫ 1 and γ ≫ 1 for the non-
SED ECR, and along a constant − ξ line with ξ ≫ 1 for the
SED ECR. These three limits are graphically represented as
arrows in Figure 1.

We aim at determining the asymptotic scalings of the loop,
star and quasi-harmonic components of the system, both in
their real and imaginary parts, along these limits. For the sake
of simplicity, the ω dependency will be explicitly extracted,
that is, the limits will be expressed in terms of ω → 0.

By following the same reasoning as above, after evaluating
the low-frequency asymptotic behavior of the interior and
exterior wavenumber, resulting in the operator scalings in real
and imaginary parts, we can estimate the asymptotic scalings
of the quasi-Helmholtz decomposed blocks of the PMCHWT
equation by enforcing the proper cancellations (equations (31),
(32)).

In addition, notice that the imaginary part of T$_{A,k}$ may
present a different scaling depending on whether the source
or test basis function is solenoidal or non-solenoidal. Indeed,
in the case of solenoidal basis functions, the second term in
the Green’s function Maclaurin expansion vanishes.

Then, from the above mentioned decomposition, the asymp-
totic scalings of the quasi-Helmholtz decomposed inverse ma-
trix Z − 1 AH Σ can be retrieved, for example by applying the Wood-
bury formula [27] or by leveraging symbolic calculations.
Numerical validation can further ensure the correctness of
the low-frequency asymptotic behavior of the quasi-Helmholtz
components of the system retrieved from this analysis. These
results are summarized in the following.

A. The Quasi-Static Regime

As already shown in Section III-A, the interior wavenumber asymptotic behavior in the quasi-static regime is ℜ { k₁L } = O (ω), ℑ { k₁L } = O (ω). This results in the Loop-Star decomposition of matrix Z in the low-frequency limit as

$$
\mathfrak{R}\{\mathbf{Z}_{\Lambda H \Sigma}\} = \mathcal{O}
\begin{pmatrix}
\omega^3 & \omega^3 & \omega^3 & \omega^2 & \omega^2 & 1 \\
\omega^3 & \omega^3 & \omega^3 & \omega^2 & 1 & 1 \\
\omega^3 & \omega^3 & \omega^{-1} & 1 & 1 & 1 \\
\omega^2 & \omega^2 & 1 & \omega & \omega & \omega \\
\omega^2 & 1 & 1 & \omega & \omega & \omega \\
1 & 1 & 1 & \omega & \omega & 1
\end{pmatrix}
\quad (52)
$$

$$
\mathfrak{I}\{\mathbf{Z}_{\Lambda H \Sigma}\} = \mathcal{O}
\begin{pmatrix}
\omega & \omega & \omega & \omega^2 & \omega^2 & \omega^2 \\
\omega & \omega & \omega & \omega^2 & \omega^2 & \omega^2 \\
\omega & \omega & \omega^{-1} & \omega^2 & \omega^2 & \omega^2 \\
\omega^2 & \omega^2 & \omega^2 & \omega & \omega & \omega \\
\omega^2 & \omega^2 & \omega^2 & \omega & \omega & \omega \\
\omega^2 & \omega^2 & \omega^2 & \omega & \omega & \omega^{-1}
\end{pmatrix} . \quad (53)
$$

The Loop-Star decomposition of the matrix inverse behaves
as

$$
\mathfrak{R}\{\mathbf{Z}_{\Lambda H \Sigma}^{-1}\} = \mathcal{O}
\begin{pmatrix}
1 & \omega & \omega & 1 & 1 & 1 \\
\omega & \omega & \omega & 1 & 1 & \omega^2 \\
\omega & \omega & \omega & 1 & \omega^2 & \omega^2 \\
1 & 1 & 1 & \omega^{-1} & \omega & \omega \\
1 & 1 & \omega^2 & \omega & \omega^2 & \omega^2 \\
1 & \omega^2 & \omega^2 & \omega & \omega^2 & \omega^2
\end{pmatrix}
\quad (54)
$$

$$
\mathfrak{I}\{\mathbf{Z}_{\Lambda H \Sigma}^{-1}\} = \mathcal{O}
\begin{pmatrix}
\omega^{-1} & \omega & \omega & 1 & \omega & \omega \\
\omega & \omega & \omega & 1 & \omega^2 & \omega^2 \\
\omega & \omega & \omega & 1 & \omega^2 & \omega^2 \\
1 & 1 & 1 & \omega^{-1} & \omega & \omega \\
\omega & \omega^2 & \omega^2 & \omega & \omega & \omega \\
\omega & \omega^2 & \omega^2 & \omega & \omega & \omega
\end{pmatrix}
. \quad (55)
$$

B. The Non-Skin-Effect-Dominated Eddy-Current Regime

The low-frequency asymptotic behavior of the interior
wavenumber in the non-SED ECR is ℜ { k₁L } = O (ω¹/²),
$\Re\{k_1 L\} = \mathcal{O}(\omega^{1/2})$ for $\omega \to 0$. This results in the following asymptotic behavior of the Loop-Star decomposed system and its inverse

$$ \Re\{Z_{\Lambda H \Sigma}\} = \mathcal{O} \begin{pmatrix} \omega^2 & \omega^2 & \omega^2 & \omega^{3/2} & \omega^{3/2} & 1 \\ \omega^2 & \omega^2 & \omega^2 & \omega^{3/2} & 1 & 1 \\ \omega^2 & \omega^2 & 1 & 1 & 1 & 1 \\ \omega^{3/2} & \omega^{3/2} & 1 & 1 & 1 & 1 \\ \omega^{3/2} & 1 & 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 & 1 & \omega^{-1/2} \end{pmatrix} \quad (56) $$

$$ \Im\{Z_{\Lambda H \Sigma}\} = \mathcal{O} \begin{pmatrix} \omega & \omega & \omega & \omega & \omega & \omega \\ \omega & \omega & \omega & \omega & \omega & \omega \\ \omega & \omega & \omega^{-1} & \omega & \omega & \omega \\ \omega & \omega & \omega & \omega & \omega & \omega \\ \omega & \omega & \omega & \omega & \omega & \omega^{-1} \end{pmatrix} \quad (57) $$

and

$$ \Re\{Z_{\Lambda H \Sigma}^{-1}\} = \mathcal{O} \begin{pmatrix} \omega^{-1/2} & 1 & \omega^{3/2} & 1 & 1 & 1 \\ 1 & 1 & \omega^2 & 1 & 1 & \omega^{3/2} \\ \omega^{3/2} & \omega^2 & \omega^2 & \omega^2 & \omega^2 & \omega^2 \\ 1 & 1 & \omega^2 & 1 & \omega^{3/2} & \omega^{3/2} \\ 1 & 1 & \omega^2 & \omega^{3/2} & \omega^{3/2} & \omega^{3/2} \\ 1 & \omega^{3/2} & \omega^2 & \omega^{3/2} & \omega^{3/2} & \omega^{3/2} \end{pmatrix} \quad (58) $$

$$ \Im\{Z_{\Lambda H \Sigma}^{-1}\} = \mathcal{O} \begin{pmatrix} \omega^{-1} & \omega^{1/2} & \omega & \omega^{1/2} & \omega^{1/2} & \omega^{1/2} \\ \omega^{1/2} & \omega & \omega & \omega & \omega & \omega \\ \omega & \omega & \omega & \omega & \omega^{5/2} & \omega^{5/2} \\ \omega^{1/2} & \omega & \omega & \omega & \omega & \omega \\ \omega^{1/2} & \omega & \omega^{5/2} & \omega & \omega & \omega \\ \omega^{1/2} & \omega & \omega^{5/2} & \omega & \omega & \omega \end{pmatrix} \quad (59) $$

### C. The Skin-Effect-Dominated Eddy-Current Regime

The product $(k_1 L)$ stays constant in the low-frequency limit considered in this case, that is, $\Re\{k_1 L\} = \mathcal{O}(1)$, $\Im\{k_1 L\} = \mathcal{O}(1)$ for $\omega \to 0$. The Loop-Star decomposition of matrix $\bar{Z}$ then behaves as

$$ \Re\{Z_{\Lambda H \Sigma}\} = \mathcal{O} \begin{pmatrix} \omega & \omega & \omega & 1 & 1 & 1 \\ \omega & \omega & \omega & 1 & 1 & 1 \\ \omega & \omega & 1 & 1 & 1 & 1 \\ 1 & 1 & \omega^{-1} & \omega^{-1} & \omega^{-1} & \omega^{-1} \\ 1 & 1 & \omega^{-1} & \omega^{-1} & \omega^{-1} & \omega^{-1} \\ 1 & 1 & \omega^{-1} & \omega^{-1} & \omega^{-1} & \omega^{-1} \end{pmatrix} \quad (60) $$

$$ \Im\{Z_{\Lambda H \Sigma}\} = \mathcal{O} \begin{pmatrix} \omega & \omega & \omega & 1 & 1 & 1 \\ \omega & \omega & \omega & 1 & 1 & 1 \\ \omega & \omega & \omega^{-1} & 1 & 1 & 1 \\ 1 & 1 & 1 & \omega^{-1} & \omega^{-1} & \omega^{-1} \\ 1 & 1 & 1 & \omega^{-1} & \omega^{-1} & \omega^{-1} \\ 1 & 1 & 1 & \omega^{-1} & \omega^{-1} & \omega^{-1} \end{pmatrix} \quad (61) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2: The magnetic frill excitation, represented as a blue arrow, is inductive (right) if there is a global loop of the geometry passing through the frill, capacitive (left) otherwise.</figcaption></figure>

The Loop-Star decomposition of the matrix inverse is

$$ \Re\{Z_{\Lambda H \Sigma}^{-1}\} = \mathcal{O} \begin{pmatrix} \omega^{-1} & \omega^{-1} & \omega & 1 & 1 & 1 \\ \omega^{-1} & \omega^{-1} & \omega & 1 & 1 & 1 \\ \omega & \omega & \omega^2 & \omega^2 & \omega^2 & \omega^2 \\ 1 & 1 & \omega^2 & \omega & \omega & \omega \\ 1 & 1 & \omega^2 & \omega & \omega & \omega \\ 1 & 1 & \omega^2 & \omega & \omega & \omega \end{pmatrix} \quad (62) $$

$$ \Im\{Z_{\Lambda H \Sigma}^{-1}\} = \mathcal{O} \begin{pmatrix} \omega^{-1} & \omega^{-1} & \omega & 1 & 1 & 1 \\ \omega^{-1} & \omega^{-1} & \omega & 1 & 1 & 1 \\ \omega & \omega & \omega & \omega^2 & \omega^2 & \omega^2 \\ 1 & 1 & \omega^2 & \omega & \omega & \omega \\ 1 & 1 & \omega^2 & \omega & \omega & \omega \\ 1 & 1 & \omega^2 & \omega & \omega & \omega \end{pmatrix} \quad (63) $$

### D. Excitation Definition and Resulting Solution

The behavior of the different components of the excitation vector are required to determine the behavior of the components of the solution of the PMCHWT.

In this work, we focus on the magnetic frill excitation [28], [29], as it is one of the most common ways to model an imposed potential difference in the circuitual framework that allows for a closed form expression of the impressed electric and magnetic fields (differently from other excitation types, such as the voltage gap).

The low-frequency asymptotic scalings of the quasi-Helmholtz components of the electric and magnetic fields produced by a magnetic frill have already been studied in [19] and are reported in the following tables. It is worth noting that, following the reasoning in [24], the harmonic component of the electric field shows a different behavior depending on whether there exists a global loop of the geometry passing through the magnetic frill or not, resulting in an inductive excitation or a capacitive one respectively (see Figure 2), both of which will be studied in the following.

By applying the inverse of the operator matrices to the resulting RHSs, it is possible to retrieve the low-frequency asymptotic scalings of the quasi-Helmholtz components of the surface current densities and, from them, to identify the dominant components, i.e., the components that are dominant in the low-frequency limits considered. The scattering of each of the current components, in their real and imaginary parts, results in a component of the electromagnetic fields. As a result of this analysis, the current components required to correctly reconstruct the fields in the three regimes for the two types of excitations considered are summarized in tables I and II, in the right column of sections (c), (d), (e), and (f).
When solving numerically the standard PMCHWT in float-
ing point arithmetic, the real parts (or imaginary parts) of
the loop and star components of the currents are stored in
the same floating point number. In the low-frequency limit,
this determines a loss of information about the non-dominant
current components, lying outside the dynamic range of the
floating point number determined by the larger current con-
tributions. Nevertheless, the lost components may be crucial
in the reconstruction of the fields. Indeed, by comparing the
current components required for the correct evaluation of the
scattered fields with the ones that can be retrieved without
proper rescaling of the quasi-Helmholtz components of the
system in the right column of section (b), one can notice
that in many cases the current resulting from a naive solution
of the PMCHWT equation does not allow for the correct
reconstruction of the fields. For example, from Table I we
infer that, in the non-SED ECR, the current component $\mathfrak{F}(j\boldsymbol{\Sigma})$
is required for a correct reconstruction of the exterior electric
field; however this component is lost in numerical cancellation
in the low-frequency limit if a proper rescaling of the current
coefficients is not applied, resulting in a deterioration of the
accuracy of the scattered field.

Now that the critical components of the current have been
identified, we can design our preconditioning strategy with
the objective of preserving them in addition to curing the ill-
conditioning of the matrix.

V. PRECONDITIONING OF THE PMCHWT IN SEVERAL
LOW-FREQUENCY REGIMES

This section aims at defining a preconditioned formulation
of the PMCHWT equation with the following objectives: (i)
stabilizing the condition number with respect to both the
electrical length and the conductivity of the scatterer in the
regimes analyzed above, and (ii) preventing the loss of dom-
inant components of the current in the low-frequency limit,
in order to guarantee the accuracy of the solution and fields.
This can be obtained by properly rescaling the quasi-Helmholtz
decomposed components of the system, by means of primal
and dual quasi-Helmholtz projectors. We propose here a quasi-
Helmholtz projectors based preconditioning strategy tailored to
the magnetic frill (or loop) excitations analyzed in tables I and
II, referred to as excitations of inductive (ind.) or capacitive
(cap.) types. The proposed formulation follows from a left-
right preconditioning of the original one, i.e.,

$$
LZRY = L \begin{pmatrix} e \\ h \end{pmatrix}, \qquad (64)
$$

where the left and right preconditioning matrices take the form

$$
L = \begin{cases}
\begin{pmatrix} \frac{1}{\sqrt{\eta_0}} \mathbf{G}^{-T} \mathbf{M}_{\text{up}} \mathbf{G}^{-1} & \mathbf{0} \\ \mathbf{0} & \sqrt{\eta_0} \mathbf{G}^{-T} \mathbf{M}_{\text{low}} \mathbf{G}^{-1} \end{pmatrix} & \text{ind. type} \\
\begin{pmatrix} \frac{1}{\sqrt{\eta_0}} \mathbf{M}_{\text{up}} & \mathbf{0} \\ \mathbf{0} & \sqrt{\eta_0} \mathbf{M}_{\text{low}} \end{pmatrix} & \text{cap. type}
\end{cases}
\tag{65}
$$

$$
R = \begin{cases}
\begin{pmatrix} \frac{1}{\sqrt{\eta_0}} \mathbf{M}_{\text{left}} & \mathbf{0} \\ \mathbf{0} & \sqrt{\eta_0} \mathbf{M}_{\text{right}} \end{pmatrix} & \text{ind. type} \\
\begin{pmatrix} \frac{1}{\sqrt{\eta_0}} \mathbf{G}^{-1} \mathbf{M}_{\text{left}} G^{-T} & \mathbf{0} \\ \mathbf{0} & \sqrt{\eta_0} \mathbf{G}^{-1} \mathbf{M}_{\text{right}} G^{-T} \end{pmatrix} & \text{cap. type}
\end{cases}
\tag{66}
$$

As already noticed in Section IV-D, the frequency behavior
of the harmonic component of the electric field produced
by a magnetic frill is different in case the excitation is of
inductive or capacitive type. To correctly enforce the required
cancellations, we employ two different, and symmetric, pre-
conditioning strategies for the inductive and capacitive cases
(equations (65), (66)). Indeed, the use of dual projectors in
the left preconditioning matrix and primal projectors in the
right preconditioning matrix for inductive excitations, and vice
versa for capacitive excitations, allows to correctly enforce the
cancellation of the static part of the electric component of the
RHS projected on the rotated non-solenoidal subspace, non
including the quasi-harmonic subspace, if the excitation is of
inductive type, and the cancellation of the static part of the
electric component of the RHS projected on the solenoidal
subspace, including the quasi-harmonic subspace, if the exci-
tation is capacitive.

The preconditioning blocks used for inductive type excita-
tions are

$$
\mathbf{M}_{\text{up}} = a_L \mathbf{P}^{\Sigma H} + b_L \mathbf{P}^{\Lambda}, \quad (67)
$$

$$
\mathbf{M}_{\text{low}} = c_L \mathbf{P}^{\Sigma H} + d_L \mathbf{P}^{\Lambda}, \quad (68)
$$

$$
\mathbf{M}_{\text{left}} = a_R \mathbf{P}^{\Lambda H} + b_R \mathbf{P}^{\Sigma}, \quad (69)
$$

$$
\mathbf{M}_{\text{right}} = c_R \mathbf{P}^{\Lambda H} + d_R \mathbf{P}^{\Sigma}. \quad (70)
$$

Symmetrically, the ones used in case of capacitive type exci-
tations are

$$
\mathbf{M}_{\text{up}} = a_L \mathbf{P}^{\Lambda H} + b_L \mathbf{P}^{\Sigma}, \quad (71)
$$

$$
\mathbf{M}_{\text{low}} = c_L \mathbf{P}^{\Lambda H} + d_L \mathbf{P}^{\Sigma} \quad (72)
$$

$$
\mathbf{M}_{\text{left}} = a_R \mathbf{P}^{\Sigma H} + b_R \mathbf{P}^{\Lambda}, \quad (73)
$$

$$
\mathbf{M}_{\text{right}} = c_R \mathbf{P}^{\Sigma H} + d_R \mathbf{P}^{\Lambda}. \quad (74)
$$

The scalar coefficients {$a_L$ - $d_R$} should be set as a function
of the frequency-conductivity regime considered, in order to
stabilize the formulation with respect to the frequency and the
conductivity independently.

The coefficients to be used in case of inductive and capaci-
tive excitations are reported in Table III. In addition, arguments
justifying these choices are delineated in the following.

The conditioning analysis of the preconditioned formulation
will be based on the asymptotic scalings of the Loop-Star
decomposition of the preconditioned matrix, defined as

$$
Z_{P,\Lambda H\Sigma} := \begin{pmatrix} \mathbf{A}^T & \mathbf{0} \\ \mathbf{0} & \mathbf{A}^T \end{pmatrix} LZR \begin{pmatrix} \mathbf{A} & \mathbf{0} \\ \mathbf{0} & \mathbf{A} \end{pmatrix}. \quad (75)
$$

A. The Quasi-Static Regime

In the QSR, the contribution of losses inside the material
does not influence the overall scalings of the quasi-Helmholtz
decomposed formulation. Indeed, the Loop-Star decomposi-
tion identified in (38) corresponds to the one of the PMCHWT
equation applied to lossless dielectric media, already identified
in [12].
<table><caption>TABLE I: Scalings of the real and imaginary parts of quantities of interest when $\omega \to 0$ for an inductive magnetic frill excitation</caption><thead><tr><th colspan="7">(a) Right hand side</th></tr></thead><tbody><tr><td>Excitation</td><td>(R, S)(E<sub>λ</sub><sup>i</sup>)</td><td>(R, S)(E<sub>H</sub><sup>i</sup>)</td><td>(R, S)(H<sub>λ</sub><sup>i</sup>)</td><td>(R, S)(H<sub>H</sub><sup>i</sup>)</td><td>(R, S)(H<sub>S</sub><sup>i</sup>)</td><td></td></tr><tr><td>Inductive</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)</td><td>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)</td></tr><tr><th colspan="7">(b) Surface current density</th></tr><tr><td>Regime</td><td>(R, S)(j<sub>λ</sub>)</td><td>(R, S)(j<sub>H</sub>)</td><td>(R, S)(j<sub>S</sub>)</td><td>(R, S)(m<sub>λ</sub>)</td><td>(R, S)(m<sub>H</sub>)</td><td>(R, S)(m<sub>S</sub>)</td></tr><tr><td>QSR</td><td>(ω, ω)</td><td>(ω, ω)</td><td>(ω, ω)</td><td>(1, 1)</td><td>(1, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)</td></tr><tr><td>non-SED ECR</td><td>(1, ω<sup>1/2</sup>)</td><td>(1, ω)</td><td>(ω<sup>2</sup>, ω)</td><td>(1, ω)</td><td>(1, ω)</td><td>(ω<sup>3/2</sup>, ω)</td></tr><tr><td>SED ECR</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td><td>(ω, ω)</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td></tr><tr><th colspan="7">(c) Electric interior near field</th></tr><tr><td>Regime</td><td>E(Rj<sub>λ</sub>)<br>E(Sj<sub>λ</sub>)</td><td>E(Rj<sub>H</sub>)<br>E(Sj<sub>H</sub>)</td><td>E(Rj<sub>S</sub>)<br>E(Sj<sub>S</sub>)</td><td>E(Rm<sub>λ</sub>)<br>E(Sm<sub>λ</sub>)</td><td>E(Rm<sub>H</sub>)<br>E(Sm<sub>H</sub>)</td><td>E(Rm<sub>S</sub>)<br>E(Sm<sub>S</sub>)</td></tr><tr><td>QSR</td><td>(ω<sup>4</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>4</sup>)</td><td>(ω<sup>4</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>4</sup>)</td><td>(1, 1)<br>(1, 1)</td><td>(1, ω<sup>2</sup>)<br>(ω<sup>2</sup>, 1)</td><td>(1, ω<sup>2</sup>)<br>(ω<sup>4</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω<sup>2</sup>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>2</sup>, ω)<br>(ω<sup>3/2</sup>, ω<sup>5/2</sup>)</td><td>(ω<sup>2</sup>, ω)<br>(ω<sup>2</sup>, ω<sup>3</sup>)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(1, ω)<br>(ω<sup>2</sup>, ω)</td><td>(1, ω)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>3/2</sup>, ω<sup>5/2</sup>)<br>(ω<sup>2</sup>, ω)</td></tr><tr><td>SED ECR</td><td>(1, 1)<br>(1, 1)</td><td>(1, 1)<br>(1, 1)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(1, 1)<br>(1, 1)</td><td>(1, 1)<br>(1, 1)</td><td>(1, 1)<br>(1, 1)</td></tr><tr><th colspan="7">(d) Magnetic interior near field</th></tr><tr><td>Regime</td><td>H(Rj<sub>λ</sub>)<br>H(Sj<sub>λ</sub>)</td><td>H(Rj<sub>H</sub>)<br>H(Sj<sub>H</sub>)</td><td>H(Rj<sub>S</sub>)<br>H(Sj<sub>S</sub>)</td><td>H(Rm<sub>λ</sub>)<br>H(Sm<sub>λ</sub>)</td><td>H(Rm<sub>H</sub>)<br>H(Sm<sub>H</sub>)</td><td>H(Rm<sub>S</sub>)<br>H(Sm<sub>S</sub>)</td></tr><tr><td>QSR</td><td>(ω, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω)</td><td>(ω, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω)</td><td>(ω, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω)</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω, ω)<br>(ω<sup>3</sup>, ω<sup>3</sup>)</td><td>(ω<sup>3</sup>, ω)<br>(ω, ω<sup>3</sup>)</td></tr><tr><td>non-SED ECR</td><td>(1, ω)<br>(ω<sup>3/2</sup>, ω<sup>1/2</sup>)</td><td>(1, ω)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(1, ω)<br>(ω<sup>2</sup>, ω)</td><td>(1, ω)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>3/2</sup>, ω<sup>1/2</sup>)<br>(1, ω)</td></tr><tr><td>SED ECR</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)<br>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)<br>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)<br>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)<br>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td><td>(ω<sup>-1</sup>, ω<sup>-1</sup>)<br>(ω<sup>-1</sup>, ω<sup>-1</sup>)</td></tr><tr><th colspan="7">(e) Electric exterior near field</th></tr><tr><td>Regime</td><td>E(Rj<sub>λ</sub>)<br>E(Sj<sub>λ</sub>)</td><td>E(Rj<sub>H</sub>)<br>E(Sj<sub>H</sub>)</td><td>E(Rj<sub>S</sub>)<br>E(Sj<sub>S</sub>)</td><td>E(Rm<sub>λ</sub>)<br>E(Sm<sub>λ</sub>)</td><td>E(Rm<sub>H</sub>)<br>E(Sm<sub>H</sub>)</td><td>E(Rm<sub>S</sub>)<br>E(Sm<sub>S</sub>)</td></tr><tr><td>QSR</td><td>(ω<sup>5</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>5</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>3</sup>, 1)<br>(1, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>3</sup>, 1)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>9</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>4</sup>, ω)<br>(ω<sup>3/2</sup>, ω<sup>9/2</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(1, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>3/2</sup>, ω<sup>9/2</sup>)<br>(ω<sup>4</sup>, ω)</td></tr><tr><td>SED ECR</td><td>(ω<sup>3</sup>, 1)<br>(1, ω<sup>3</sup>)</td><td>(ω<sup>3</sup>, 1)<br>(1, ω<sup>3</sup>)</td><td>(ω<sup>3</sup>, 1)<br>(1, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>3</sup>, 1)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>3</sup>, 1)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>3</sup>, 1)</td></tr><tr><th colspan="7">(f) Magnetic exterior near field</th></tr><tr><td>Regime</td><td>H(Rj<sub>λ</sub>)<br>H(Sj<sub>λ</sub>)</td><td>H(Rj<sub>H</sub>)<br>H(Sj<sub>H</sub>)</td><td>H(Rj<sub>S</sub>)<br>H(Sj<sub>S</sub>)</td><td>H(Rm<sub>λ</sub>)<br>H(Sm<sub>λ</sub>)</td><td>H(Rm<sub>H</sub>)<br>H(Sm<sub>H</sub>)</td><td>H(Rm<sub>S</sub>)<br>H(Sm<sub>S</sub>)</td></tr><tr><td>QSR</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)<br>(ω, ω<sup>4</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω, ω<sup>4</sup>)</td><td>(ω, ω<sup>4</sup>)<br>(ω, ω<sup>4</sup>)</td></tr><tr><td>non-SED ECR</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>7/2</sup>, ω<sup>1/2</sup>)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>7/2</sup>, ω<sup>1/2</sup>)<br>(1, ω<sup>3</sup>)</td></tr><tr><td>SED ECR</td><td>(ω<sup>-1</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>-1</sup>)</td><td>(ω<sup>-1</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>-1</sup>)</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)<br>(ω, ω<sup>4</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω, ω<sup>4</sup>)</td><td>(ω<sup>2</sup>, ω<sup>-1</sup>)<br>(ω<sup>-1</sup>, ω<sup>2</sup>)</td></tr><tr><th colspan="7">(g) Rescaled current density</th></tr><tr><td>Regime</td><td>a<sub>R</sub><sup>-1</sup>j<sub>λ</sub></td><td>a<sub>R</sub><sup>-1</sup>j<sub>H</sub></td><td>b<sub>R</sub><sup>-1</sup>j<sub>S</sub></td><td>c<sub>R</sub><sup>-1</sup>m<sub>λ</sub></td><td>c<sub>R</sub><sup>-1</sup>m<sub>H</sub></td><td>d<sub>R</sub><sup>-1</sup>m<sub>S</sub></td></tr><tr><td>QSR</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, ω<sup>2</sup>)</td><td>(1, 1)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>1/2</sup>, ω)</td><td>(ω<sup>1/2</sup>, ω<sup>3/2</sup>)</td><td>(ω, 1)</td><td>(1, ω)</td><td>(1, ω)</td><td>(ω, ω<sup>1/2</sup>)</td></tr><tr><td>SED ECR</td><td>(ω<sup>-1/2</sup>, ω<sup>-1/2</sup>)</td><td>(ω<sup>-1/2</sup>, ω<sup>-1/2</sup>)</td><td>(ω<sup>-1/2</sup>, ω<sup>-1/2</sup>)</td><td>(ω<sup>-1/2</sup>, ω<sup>-1/2</sup>)</td><td>(ω<sup>-1/2</sup>, ω<sup>-1/2</sup>)</td><td>(ω<sup>-1/2</sup>, ω<sup>-1/2</sup>)</td></tr></tbody></table>

1) Inductive type excitation: Different sets of coefficients leading to the Loop-Star decomposition allow for frequency and conductivity preconditioning of the formulation. One classical choice [12] is for example

$$a_L = c_L = b_R = d_R = \chi^{1/2} \text{ and } b_L = d_L = a_R = c_R = \chi^{-1/2},$$

$$Z_{P,\Lambda H\Sigma} = O \begin{pmatrix} 1 & 1 & \chi & \chi & \chi & 1 \\ \chi & \chi & \chi^2 & \chi^2 & 1 & \chi \\ \chi & \chi & 1 & 1 & 1 & \chi \\ \chi & \chi & 1 & 1 & 1 & \chi \\ \chi^2 & 1 & \chi & \chi & \chi & \chi^2 \\ 1 & 1 & \chi & \chi & \chi & 1 \end{pmatrix} \quad (76)$$
<table><caption>TABLE II: Scalings of the real and imaginary parts of quantities of interest when $\omega \to 0$ for a capacitive magnetic frill excitation</caption><thead><tr><th></th><th colspan="7">(a) Right hand side</th></tr></thead><tbody><tr><td>Excitation</td><td>(R, S)(E<sub>λ</sub><sup>i</sup>)</td><td>(R, S)(E<sub>H</sub><sup>i</sup>)</td><td>(R, S)(E<sub>Σ</sub><sup>i</sup>)</td><td>(R, S)(H<sub>λ</sub><sup>i</sup>)</td><td>(R, S)(H<sub>H</sub><sup>i</sup>)</td><td>(R, S)(H<sub>Σ</sub><sup>i</sup>)</td><td></td></tr><tr><td>Capacitive</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)</td><td>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)</td><td></td></tr><tr><th></th><th colspan="7">(b) Surface current density</th></tr><tr><td>Regime</td><td>(R, S)(j<sub>λ</sub>)</td><td>(R, S)(j<sub>H</sub>)</td><td>(R, S)(j<sub>Σ</sub>)</td><td>(R, S)(m<sub>λ</sub>)</td><td>(R, S)(m<sub>H</sub>)</td><td>(R, S)(m<sub>Σ</sub>)</td><td>Current dominant components</td></tr><tr><td>QSR</td><td>(ω, ω)</td><td>(ω, ω)</td><td>(ω, ω)</td><td>(1, 1)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>3/2</sup>, ω)</td><td>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>5/2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5/2</sup>)</td><td>∫(j<sub>λ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>SED ECR</td><td>(ω, ω)</td><td>(ω, ω)</td><td>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)</td></tr><tr><th></th><th colspan="7">(c) Electric interior near field</th></tr><tr><td>Regime</td><td>E(Rj<sub>λ</sub>)<br>E(Sj<sub>λ</sub>)</td><td>E(Rj<sub>H</sub>)<br>E(Sj<sub>H</sub>)</td><td>E(Rj<sub>Σ</sub>)<br>E(Sj<sub>Σ</sub>)</td><td>E(Rm<sub>λ</sub>)<br>E(Sm<sub>λ</sub>)</td><td>E(Rm<sub>H</sub>)<br>E(Sm<sub>H</sub>)</td><td>E(Rm<sub>Σ</sub>)<br>E(Sm<sub>Σ</sub>)</td><td>Current dominant components</td></tr><tr><td>QSR</td><td>(ω<sup>4</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>4</sup>)</td><td>(ω<sup>4</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>4</sup>)</td><td>(1, 1)<br>(1, 1)</td><td>(1, ω<sup>2</sup>)<br>(ω<sup>2</sup>, 1)</td><td>(ω<sup>2</sup>, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>Σ</sub>, m<sub>λ</sub>)<br>∫(j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>7/2</sup>, ω<sup>5/2</sup>)<br>(ω<sup>2</sup>, ω<sup>3</sup>)</td><td>(ω<sup>4</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω<sup>3</sup>)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>7/2</sup>, ω<sup>5/2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>7/2</sup>, ω<sup>5/2</sup>)</td><td>∫(j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>SED ECR</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>3</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)</td></tr><tr><th></th><th colspan="7">(d) Magnetic interior near field</th></tr><tr><td>Regime</td><td>H(Rj<sub>λ</sub>)<br>H(Sj<sub>λ</sub>)</td><td>H(Rj<sub>H</sub>)<br>H(Sj<sub>H</sub>)</td><td>H(Rj<sub>Σ</sub>)<br>H(Sj<sub>Σ</sub>)</td><td>H(Rm<sub>λ</sub>)<br>H(Sm<sub>λ</sub>)</td><td>H(Rm<sub>H</sub>)<br>H(Sm<sub>H</sub>)</td><td>H(Rm<sub>Σ</sub>)<br>H(Sm<sub>Σ</sub>)</td><td>Current dominant components</td></tr><tr><td>QSR</td><td>(ω, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω)</td><td>(ω, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω)</td><td>(ω, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω)</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω<sup>3</sup>, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω<sup>3</sup>)</td><td>(ω<sup>3</sup>, ω)<br>(ω<sup>3</sup>, ω)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>Σ</sub>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>3/2</sup>, ω<sup>5/2</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>3</sup>)<br>(ω<sup>7/2</sup>, ω<sup>5/2</sup>)</td><td>(ω<sup>2</sup>, ω)<br>(ω<sup>3/2</sup>, ω<sup>5/2</sup>)</td><td>∫(m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>SED ECR</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω<sup>2</sup>, ω<sup>2</sup>)<br>(ω, ω)</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω, ω)<br>(ω, ω)</td><td>(ω, ω)<br>(ω, ω)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)</td></tr><tr><th></th><th colspan="7">(e) Electric exterior near field</th></tr><tr><td>Regime</td><td>E(Rj<sub>λ</sub>)<br>E(Sj<sub>λ</sub>)</td><td>E(Rj<sub>H</sub>)<br>E(Sj<sub>H</sub>)</td><td>E(Rj<sub>Σ</sub>)<br>E(Sj<sub>Σ</sub>)</td><td>E(Rm<sub>λ</sub>)<br>E(Sm<sub>λ</sub>)</td><td>E(Rm<sub>H</sub>)<br>E(Sm<sub>H</sub>)</td><td>E(Rm<sub>Σ</sub>)<br>E(Sm<sub>Σ</sub>)</td><td>Current dominant components</td></tr><tr><td>QSR</td><td>(ω<sup>5</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>5</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>3</sup>, 1)<br>(1, ω<sup>3</sup>)</td><td>(1, ω<sup>3</sup>)<br>(ω<sup>3</sup>, 1)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>Σ</sub>, m<sub>λ</sub>)<br>∫(j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>11/2</sup>, ω<sup>5/2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>6</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(1, ω<sup>3</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>11/2</sup>, ω<sup>5/2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>11/2</sup>, ω<sup>5/2</sup>)</td><td>∫(j<sub>Σ</sub>)</td></tr><tr><td>SED ECR</td><td>(ω<sup>5</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>5</sup>, ω<sup>2</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(1, ω<sup>3</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>5</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>Σ</sub>)</td></tr><tr><th></th><th colspan="7">(f) Magnetic exterior near field</th></tr><tr><td>Regime</td><td>H(Rj<sub>λ</sub>)<br>H(Sj<sub>λ</sub>)</td><td>H(Rj<sub>H</sub>)<br>H(Sj<sub>H</sub>)</td><td>H(Rj<sub>Σ</sub>)<br>H(Sj<sub>Σ</sub>)</td><td>H(Rm<sub>λ</sub>)<br>H(Sm<sub>λ</sub>)</td><td>H(Rm<sub>H</sub>)<br>H(Sm<sub>H</sub>)</td><td>H(Rm<sub>Σ</sub>)<br>H(Sm<sub>Σ</sub>)</td><td>Current dominant components</td></tr><tr><td>QSR</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>4</sup>, ω)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>6</sup>, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω<sup>6</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω, ω<sup>4</sup>)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>Σ</sub>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>3/2</sup>, ω<sup>9/2</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>6</sup>, ω<sup>3</sup>)<br>(ω<sup>2</sup>, ω<sup>5</sup>)</td><td>(ω<sup>6</sup>, ω<sup>3</sup>)<br>(ω<sup>7/2</sup>, ω<sup>13/2</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω<sup>3/2</sup>, ω<sup>9/2</sup>)</td><td>∫(m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>)</td></tr><tr><td>SED ECR</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω, ω<sup>4</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>2</sup>, ω<sup>5</sup>)<br>(ω<sup>4</sup>, ω)</td><td>(ω<sup>6</sup>, ω<sup>3</sup>)<br>(ω<sup>3</sup>, ω<sup>6</sup>)</td><td>(ω<sup>4</sup>, ω)<br>(ω, ω<sup>4</sup>)</td><td>(ω, ω<sup>4</sup>)<br>(ω, ω<sup>4</sup>)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, m<sub>λ</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)</td></tr><tr><th></th><th colspan="7">(g) Rescaled current density</th></tr><tr><td>Regime</td><td>b<sub>R</sub><sup>-1</sup>j<sub>λ</sub></td><td>a<sub>R</sub><sup>-1</sup>j<sub>H</sub></td><td>a<sub>R</sub><sup>-1</sup>j<sub>Σ</sub></td><td>d<sub>R</sub><sup>-1</sup>m<sub>λ</sub></td><td>c<sub>R</sub><sup>-1</sup>m<sub>H</sub></td><td>c<sub>R</sub><sup>-1</sup>m<sub>Σ</sub></td><td>Recovered components</td></tr><tr><td>QSR</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td><td>(1, 1)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)</td></tr><tr><td>non-SED ECR</td><td>(ω<sup>2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>5/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>5/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>5/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>3/2</sup>, ω<sup>2</sup>)</td><td>(ω<sup>3/2</sup>, ω<sup>2</sup>)</td><td>∫(j<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>)</td></tr><tr><td>SED ECR</td><td>(ω<sup>3/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>3/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>5/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>3/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>3/2</sup>, ω<sup>3/2</sup>)</td><td>(ω<sup>3/2</sup>, ω<sup>3/2</sup>)</td><td>∫(j<sub>λ</sub>, j<sub>H</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)<br>∫(j<sub>λ</sub>, j<sub>H</sub>, j<sub>Σ</sub>, m<sub>λ</sub>, m<sub>H</sub>, m<sub>Σ</sub>)</td></tr></tbody></table>

leading to a constant condition number behavior in the low-frequency limit [12].

However, for the kind of excitation considered here, this choice would not allow for the preservation of all the relevant current components. Indeed, the components $j_λ$, $j_H$, and $m_Σ$, both in their real and imaginary parts, which are useful to reconstruct the magnetic field inside and outside the scatterer,

would be lost in numerical cancellation in the low-frequency limit.

To preserve the accuracy of all the required current components, we opt instead to the choice of coefficients reported in
<table><caption>TABLE III: Scalar multiplicative coefficients defining the proposed preconditioning strategy.</caption><thead><tr><th></th><th>QSR-ind</th><th>QSR-cap</th><th>non-SED ECR-ind</th><th>non-SED ECR-cap</th><th>SED ECR-ind</th><th>SED ECR-cap</th></tr></thead><tbody><tr><td>a<sub>L</sub></td><td>1</td><td>χ<sup>-2</sup></td><td>χ<sup>1/2</sup>γ<sup>-1</sup></td><td>χ<sup>-1/2</sup></td><td>χ<sup>-1/2</sup></td><td>χ<sup>-1/2</sup></td></tr><tr><td>b<sub>L</sub></td><td>χ<sup>-2</sup></td><td>1</td><td>χ<sup>-1/2</sup></td><td>χ<sup>3/2</sup></td><td>χ<sup>-1/2</sup></td><td>χ<sup>3/2</sup></td></tr><tr><td>c<sub>L</sub></td><td>χ<sup>-1</sup></td><td>χ<sup>-1</sup></td><td>χ<sup>1/2</sup></td><td>χ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td></tr><tr><td>d<sub>L</sub></td><td>χ<sup>-1</sup></td><td>χ<sup>-1</sup></td><td>χ<sup>-1/2</sup>γ</td><td>χ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td></tr><tr><td>a<sub>R</sub></td><td>χ</td><td>χ</td><td>χ<sup>-1/2</sup></td><td>χ<sup>-1/2</sup></td><td>χ<sup>-1/2</sup></td><td>χ<sup>-1/2</sup></td></tr><tr><td>b<sub>R</sub></td><td>χ</td><td>χ</td><td>χ<sup>1/2</sup>γ</td><td>χ<sup>-1/2</sup></td><td>χ<sup>3/2</sup></td><td>χ<sup>-1/2</sup></td></tr><tr><td>c<sub>R</sub></td><td>1</td><td>χ<sup>2</sup></td><td>χ<sup>-1/2</sup>γ</td><td>χ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td></tr><tr><td>d<sub>R</sub></td><td>χ<sup>2</sup></td><td>1</td><td>χ<sup>1/2</sup></td><td>χ<sup>-3/2</sup>γ<sup>2</sup></td><td>γ<sup>1/2</sup></td><td>γ<sup>1/2</sup></td></tr></tbody></table>

<caption>Table III, resulting in the Loop-Star decomposition</caption>

$$Z_{P,\Lambda H\Sigma} = \mathcal{O} \begin{pmatrix} 1 & 1 & 1 & 1 & 1 & 1 \\ \chi^2 & \chi^2 & \chi^2 & \chi^2 & 1 & \chi^2 \\ \chi^2 & \chi^2 & 1 & 1 & 1 & \chi^2 \\ \chi^2 & \chi^2 & 1 & 1 & 1 & \chi^2 \\ \chi^2 & 1 & 1 & 1 & 1 & \chi^2 \\ 1 & 1 & 1 & 1 & 1 & 1 \end{pmatrix}, \quad (77)$$

indicating the favorable condition properties of the preconditioned matrix and that the matrix has a well defined low-frequency limit converging to the static problem. The preserved current components correspond to the ones needed for the correct fields reconstruction, as shown in table I.

2) *Capacitive type excitation:* Similarly as above, we determine from table I the components to be preserved and design the multiplicative preconditioner accordingly. Our choice of coefficients, in table III, results in

$$Z_{P,\Lambda H\Sigma} = \mathcal{O} \begin{pmatrix} 1 & 1 & 1 & 1 & \chi^2 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \chi^2 & \chi^2 & 1 & 1 & \chi^2 & \chi^2 \\ \chi^2 & \chi^2 & 1 & 1 & \chi^2 & \chi^2 \\ \chi^2 & 1 & 1 & 1 & \chi^2 & \chi^2 \\ 1 & 1 & 1 & 1 & \chi^2 & 1 \end{pmatrix}, \quad (78)$$

which indicates that the resulting matrix has favorable conditioning properties, and allows to preserve the required current components (see Table II).

## B. The Non-Skin-Effect-Dominated Eddy-Current Regime

1) *Inductive type excitation:* One first intuitive choice for the scalar coefficients {$a_L - d_R$} is

$$a_L = c_L = b_R = d_R = \chi^{1/2}, \quad b_L = a_R = \chi^{-1/2}, \\ \text{and } d_L = c_R = \chi^{-1/2}\gamma.$$

This would allow to retrieve all the dominant current components. Moreover, this would result in the decomposition

$$Z_{P,\Lambda H\Sigma} = \mathcal{O} \begin{pmatrix} 1 & 1 & \chi & \xi & \xi & 1 \\ \chi & \chi & \chi^2 & \chi\xi & \gamma & \chi \\ \chi & \chi & 1 & \gamma & \gamma & \chi \\ \xi & \xi & \gamma & 1 & 1 & \xi \\ \xi^2 & 1 & \chi & \xi & \xi & \xi^2 \\ 1 & 1 & \chi & \xi & \xi & 1 \end{pmatrix}, \quad (79)$$

In presence of multiply connected geometries with finite genus *g* instead, the preconditioned matrices exhibit a finite-dimensional nullspace in the static limit, of size $2g$. This

can be seen for example by studying the Loop-Star decomposition for a genus-1 geometry and separating the toroidal from the poloidal global loop contributions (see [24] for further details on this technique), i.e., by replacing matrix $A$ with $A' := (\Lambda H^T H^P \Sigma)$, where $H^T$ and $H^P$ are the transformation matrices from the toroidal and poloidal global loops to the RWG subspace. The decomposition of the preconditioned matrix reads

$$\begin{pmatrix} A'^T & 0 \\ 0 & A'^T \end{pmatrix} L Z R \begin{pmatrix} A' & 0 \\ 0 & A' \end{pmatrix} = \mathcal{O} \begin{pmatrix} 1 & 1 & 1 & \chi & \xi & \xi & \xi & 1 \\ \chi & \chi & \chi & \chi^2 & \chi\xi & \xi^2\gamma & \gamma & \chi \\ \chi & \chi & \chi & \chi^2 & \chi\xi & \xi^2\gamma & \xi^2\gamma & \chi \\ \chi & \chi & \chi & 1 & \gamma & \gamma & \gamma & \chi \\ \xi & \xi & \xi & \gamma & 1 & 1 & 1 & \xi \\ \xi^2 & \xi^2 & 1 & \chi & \xi & \xi & \xi & \xi^2 \\ \xi^2 & \xi^2 & \xi^2 & \chi & \xi & \xi & \xi & \xi^2 \\ 1 & 1 & 1 & \chi & \xi & \xi & \xi & 1 \end{pmatrix}, \quad (80)$$

where the terms in bold take into account the fact that the static part of the MFIO matrix vanishes when the basis loop does not pass through the test loop [24]. The application of (80) to the $V$-dimensional (with $V$ denoting the number of vertices of the mesh) space of magnetic currents composed of local RWG loops and the toroidal global loop, that is, the vector

$$(\mathbf{0} \quad \mathbf{0} \quad \mathbf{0} \quad \mathbf{0} \quad \mathbf{m}_{\Lambda} \quad \mathbf{m}_{H^T} \quad \mathbf{0} \quad \mathbf{0})^{\mathrm{T}} \quad (81)$$

results in

$$(\mathcal{O}(\xi) \quad \mathcal{O}(\chi\xi) \quad \mathcal{O}(\chi\xi) \quad \mathcal{O}(\gamma) \quad \mathbf{v} \quad \mathcal{O}(\xi) \quad \mathcal{O}(\xi) \quad \mathcal{O}(\xi))^{\mathrm{T}}, \quad (82)$$

where $\mathbf{v}$ is given by the product between a rank-deficient $(V-1) \times V$ block and the $V$-dimensional current. As $\chi \ll 1$, $\gamma \ll 1$, and $\xi \ll 1$ in the non-SED ECR, the nullspace of such rank-deficient block corresponds to the nullspace of the entire matrix in the low-frequency limit. Similar results hold for the $V$-dimensional magnetic current space composed of local RWG loops and the poloidal global loop, resulting, overall, in a two-dimensional nullspace of the preconditioned formulation, which, numerically, takes the form of two isolated singular values decaying toward lower frequencies. Note that a similar analysis applied to the preconditioned formulation in the QSR confirms the absence of a nullspace related to the global loops of the geometry.
To remove the artificial nullspace introduced by this sub-optimal preconditioning, we further manipulate the set of coefficients {$a_L - d_R$}. The choice reported in table III on the one hand allows the retrieval of all the useful current components—as shown in table I—and on the other hand results in the behavior

$$Z_{P,\Lambda H\Sigma} = O \begin{pmatrix} 1 & 1 & \chi\gamma & \xi & \xi & 1 \\ \xi & \xi & \chi^2 & \xi^2 & 1 & \xi \\ \xi & \xi & 1 & 1 & 1 & \xi \\ \xi & \xi & \gamma^2 & 1 & 1 & \xi \\ \xi^2 & 1 & \chi\gamma & \xi & \xi & \xi^2 \\ 1 & 1 & \chi\gamma & \xi & \xi & 1 \end{pmatrix}, \quad (83)$$

indicating favorable conditioning properties. Moreover, by applying an analysis similar to the one above, the absence of a nullspace in the low-frequency limit related to the non-zero genus of the geometry can be verified.

2) **Capacitive type excitation:** An effective choice of coefficients to retrieve all the useful current components is the one reported in table III, resulting in the Loop-Star decomposition

$$Z_{P,\Lambda H\Sigma} = O \begin{pmatrix} 1 & 1 & 1 & 1 & \xi^2 & 1 \\ 1 & 1 & 1 & 1 & 1 & 1 \\ \chi^2 & \chi^2 & 1 & \gamma^2 & \chi^2 & \chi^2 \\ \xi^2 & \xi^2 & 1 & 1 & \xi^2 & \xi^2 \\ \xi^2 & 1 & 1 & 1 & \xi^2 & \xi^2 \\ 1 & 1 & 1 & 1 & \xi^2 & 1 \end{pmatrix}, \quad (84)$$

indicating that the formulation has a well-defined low-frequency limit and favorable conditioning properties.

### C. The Skin-Effect-Dominated Eddy-Current Regime

1) **Inductive type excitation:** As indicated in Table I, all the current components, in real and imaginary parts, are required to correctly retrieve the fields in the low-frequency limit. This requirement is satisfied by the choice of coefficients in Table III, resulting in the following behavior of the Loop-Star decomposed system

$$Z_{P,\Lambda H\Sigma} = O \begin{pmatrix} 1 & 1 & \chi^2 & \xi^{-1/2} & \xi^{-1/2} & \xi^{-1/2} \\ 1 & 1 & \chi^2 & \xi^{-1/2} & \xi^{-1/2} & \xi^{-1/2} \\ 1 & 1 & 1 & \xi^{-1/2} & \xi^{-1/2} & \xi^{-1/2} \\ \xi^{-1/2} & \xi^{-1/2} & \chi^{3/2}\gamma^{1/2} & 1 & 1 & 1 \\ \xi^{-1/2} & \xi^{-1/2} & \chi^{3/2}\gamma^{1/2} & 1 & 1 & 1 \\ \xi^{-1/2} & \xi^{-1/2} & \chi^{3/2}\gamma^{1/2} & 1 & 1 & 1 \end{pmatrix}, \quad (85)$$

attesting to the favorable conditioning properties of the formulation.

2) **Capacitive type excitation:** In this case, the preconditioning strategy implemented for the inductive type excitation can be symmetrized, by exchanging the role of left and right coefficients (see Table III), resulting in the same preconditioned formulation as in (85), which also allows for the retrieval of the required current components.

Note that the asymptotic analyses of the Loop-Star decomposed preconditioned systems shown above in (77), (78), (83), (84), and (85) provide an indication of the favorable conditioning properties of the preconditioned system, but cannot

provide a formal answer to the question of whether or not the low-frequency limit of the preconditioned matrices suffer from a finite dimensional nullspace. Nevertheless, numerical results presented in Section VI (see Figures 3c and 3d) suggest the well-posedness and stability of the proposed formulations.

## VI. NUMERICAL RESULTS

In this section, numerical results will illustrate the favorable spectral properties, the stability, and the accuracy of the proposed formulation. First, the proposed formulation is applied to a multiply connected toroidal structure and its condition number is compared to that of the PMCHWT equation (matrix $Z$) and of its rescaled version (matrix $\bar{Z}$) in a wide range of conductivities and frequencies (always in the low-frequency regime), both for inductive and capacitive-type excitation versions (Figure 3). On the one hand the low-frequency breakdown of the standard PMCHWT matrices $Z$ or $\bar{Z}$ is clear (as well as other conductivity related instabilities mentioned above), on the other hand the stability of the proposed formulation, that exhibits a condition number that remains constant in frequency and conductivity in the three regimes considered. To produce these results, we applied the preconditioning strategies proposed for the QSR, non-SED ECR, and SED ECR in the corresponding frequency and conductivity regions, identified by equations (34), (39), (44) by substituting the symbols $\gg$, $\ll$ with $>$, $<$. The demarcation between these regions has been marked with black dotted lines.

Next, we verify that the quasi-Helmholtz projectors based preconditioning does not degrade the spectral properties of the formulation with respect to mesh refinement by comparing the condition number of the standard PMCHWT matrix and of the proposed formulation with the one of a PMCHWT formulation preconditioned with a Loop-Star approach. The rescaling coefficients employed in the Loop-Star preconditioning are the same as the one used in our formulation, in Table III. The results, presented in Figure 4, show that the condition number of matrix $Z$ and the new formulation based on quasi-Helmholtz projector increase as $h^{-2}$ with refinement, which confirms the fact that the use of quasi-Helmholtz preconditioning is not detrimental with respect to the PMCHWT dense discretization breakdown. The condition number of the Loop-Star preconditioned formulation, instead, grows as $h^{-4}$.

Beyond the conditioning properties of the equations under study, we next need to assess their accuracy. We consider two canonical circuital structures and employ our full-wave solver to evaluate their impedance, of both inductive and capacitive type. These values, derived from the voltages and currents themselves computed from the simulated total fields, are then compared to the values predicted from circuit theory to verify the accuracy of the schemes. The following results were obtained with a direct solver.

The first structure considered is a torus with major and minor radii $R_M = 1$ m and $R_m = 0.2$ m excited by an inductive magnetic frill passing through the global toroidal loop of the geometry, which imposes a potential difference of 1 V. The conductivity of the toroidal body $\Omega$ varies from
<figure><img src="image_4.png" /><figcaption>Fig. 3: Condition number of the formulations applied to torus with major and minor radii of 1 m and 0.3 m discretized at $h \simeq 0.44$ m. The preconditioning applied is tailored to inductive (c) and capacitive (d) excitations.</figcaption></figure>

$1 \times 10^{-3}$ S/m to $1 \times 10^3$ S/m. We estimate the resistance and inductance of the circuit from the electromagnetic field scattered by the equivalent currents retrieved from the solution of our formulation and compare these values with the approximation valid in low frequency for conductors with uniform cross section $R_{ct} := 2R_M/(R_m^2\sigma)$ [30], and $L_{ct} \simeq 1.807$ µH [31, Eq. 4.27]. The relative difference between the results obtained from our full-wave formulation and the ones known from circuit theory, $R_{ct}$ and $L_{ct}$, is reported in Figure 5 for a wide range of frequencies. In particular, the relative error of the resistance is approximately $2.5 \times 10^{-3}$, the one of the inductance $4.2 \times 10^{-4}$, both of them independently of frequency and conductivity of the line. Figure 6 represents the magnetic field induced by the current circulating in the loop. This induced field is parallel to the torus axis and symmetric with respect to the center.

The second structure instead is a parallel plates capacitor, made up of a couple of circular parallel plates of radius

$R_p = 4$ m placed at distance $d = 0.2$ m. The plates are connected by a wire of circular cross section of radius $0.2$ m. A potential difference of $1$ V is imposed across the plates using a capacitive magnetic frill around this wire. The capacitance of the structure is computed from the simulated fields, and compared to the value of capacitance expected from circuit theory $C_{ct} := \epsilon_0\pi R_p^2/d$ in a wide range of frequencies and conductivities of the circuit. The results are summarized in Figure 7, showing a relative difference between the two values of capacitance of approximately $7.7 \times 10^{-2}$, irrespectively of frequency and conductivity. Figure 8 shows a representation of the electric field generated between the plates. This electric field is normally directed and uniform at a value near to $5$ V/m.

In the final numerical example, we demonstrate that our formulation can effectively track the field penetration depth variation inside a conductor. We considered in this case an highly conductive wire ($\sigma = 1 \times 10^7$ S/m) of rectangular
<figure><img src="image_5.png" /><figcaption>Fig. 4: Condition number of the formulations applied to a sphere of radius 1 m of conductivity 10 S/m excited at 1 × 10⁵ Hz at different mesh refinements: comparison between the original PMCHWT formulation (Z), the Loop-Star preconditioned formulation (ΛΗΣ precond.), and this work (LZR).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 6: Torus of major and minor radii of 1 m and 0.2 m of conductivity 1 S/m excited by a magnetic frill imposing 1 V at 10 mHz: induced magnetic field.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5: Full-wave evaluation of the impedance of a toroidal structure with major and minor radii of 1 m and 0.2 m of different conductivities: relative difference with respect to circuit theory expectations $R_{ct}$ and $L_{ct}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Full-wave evaluation of the impedance of a capacitive structure with circular parallel plates of radius 4 m at distance 0.2 m of different conductivities: relative difference with respect to circuit theory expectation $C_{ct}$.</figcaption></figure>

cross section enclosed in square shape with axis of symmetry directed toward $\hat{y}$. After exciting the structure with a plane wave traveling along $\hat{y}$, the resulting current density distribution in a cross section of the line at $z = 0$ m is shown in Figure 9 for different frequencies of the impinging wave. Despite the use of the preconditioning strategy tailored for inductive magnetic frill excitations, the solution is expected to be accurate even in presence of plane wave excitations if the operating frequency is not extremely low. We notice that the visible penetration length, which decreases when increasing the frequency, is comparable with the skin depth $\delta$, which is also reported in the figures as the black lines for ease of comparison.

## VII. CONCLUSION

In this work we presented a novel full-wave formulation amenable to eddy-current simulations in circuital frameworks. Built upon the PMCHWT integral equation, it results from a quasi-Helmholtz preconditioning strategy of projective type. With respect to the standard Loop-Star decomposition, the quasi-Helmholtz projectors offer several advantages, including their compatibility with fast solvers and acceleration strategies and their lack of detrimental effect on the condition number with respect to dense discretizations. The proposed preconditioning strategy offers the twofold advantage of curing the frequency and conductivity instabilities of the original formulation, including the severe low-frequency breakdown preventing its application toward statics, and curing the typical
<figure><img src="image_9.png" /><figcaption>Fig. 8: Capacitive structure with circular parallel plates of radius 4 m at distance 0.2 m of conductivity $1 \times 10^{-3}$ S/m excited by a magnetic frill imposing 1 V at 1 Hz: electric field between the plates.</figcaption></figure>

loss of accuracy occurring at very low-frequency. The resulting scheme is well-conditioned and stable for materials ranging from dielectric to highly conductive, applicable to both simply and multiply-connected structures and prone to seamless transitions between lower to higher frequencies. Although the proposed scheme provides optimal accuracy in presence of magnetic frill type excitations, the procedure employed in defining the proposed preconditioning strategy can also be applied to other types of excitations, which will be topic of future investigations.

REFERENCES

[1] C. A. Balanis, *Advanced Engineering Electromagnetics*. Hoboken, N.J: John Wiley & Sons, 2nd ed., 2012.

[2] H. K. Dirks, "Quasi-stationary fields for microelectronic applications," *Electrical Engineering*, vol. 79, pp. 145-155, Apr. 1996.

[3] S. Sauter and C. Schwab, *Boundary Element Methods*, vol. 39 of Springer Series in Computational Mathematics. Berlin, Heidelberg: Springer Berlin Heidelberg, 2011.

[4] O. Steinbach, *Numerical Approximation Methods for Elliptic Boundary Value Problems: Finite and Boundary Elements*. New York: Springer, 2008.

[5] A. Poggio and E. Miller, "Integral equation solutions of three-dimensional scattering problems," in *Computer Techniques for Electromagnetics*, pp. 159-264, Elsevier, 1973.

[6] C. Muller, "Grundprobleme der mathematischen theorie elektromagnetischer schwingungen," *The American Mathematical Monthly*, vol. 65, p. 459, June 1958.

[7] I. Bogaert, K. Cools, F. Andriulli, and H. Bagci, "Low-frequency scaling of the standard and mixed magnetic field and Müller integral equations," *IEEE Transactions on Antennas and Propagation*, vol. 62, pp. 822-831, Feb. 2014.

[8] P. Ylä-Oijala, M. Taskinen, and S. Järvenpää, "Analysis of surface integral equations in electromagnetic scattering and radiation problems," *Engineering Analysis with Boundary Elements*, vol. 32, pp. 196-209, Mar. 2008.

[9] E. van 't Wout, S. Haqshenas, P. Gélat, T. Betcke, and N. Saffari, "Boundary integral formulations for acoustic modelling of high-contrast media," *Computers & Mathematics with Applications*, vol. 105, pp. 136-149, Jan. 2022.

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

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

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

<figure><img src="image_13.png" /><figcaption>Fig. 9: Electric current density inside a rectangular cross section of sides 1 cm and 0.5 cm of conductivity $1 \times 10^7$ S/m excited by an impinging plane wave travelling along $\hat{y}$ at frequencies (a) 8 kHz, (b) 16 kHz, (c) 32 kHz, and (d) 64 kHz.</figcaption></figure>
[10] S. Adrian, A. Dely, D. Consoli, A. Merlini, and F. Andriulli, "Electromagnetic integral equations: Insights in conditioning and preconditioning," *IEEE Open Journal of Antennas and Propagation*, vol. 2, pp. 1143-1174, 2021.

[11] K. Cools, F. Andriulli, and E. Michielssen, "A Calderón multiplicative preconditioner for the PMCHWT integral equation," *IEEE Transactions on Antennas and Propagation*, vol. 59, pp. 4579-4587, Dec. 2011.

[12] Y. Beghein, R. Mitharwal, K. Cools, and F. Andriulli, "On a low-frequency and refinement stable PMCHWT integral equation leveraging the quasi-Helmholtz projectors," *IEEE Transactions on Antennas and Propagation*, vol. 65, pp. 5365-5375, Oct. 2017.

[13] R. Hiptmair, "Boundary element methods for eddy current computation," pp. 213-248, 2007.

[14] E. Kriezis, T. Tsiboukis, S. Panas, and J. Tegopoulos, "Eddy currents: Theory and applications," *Proceedings of the IEEE*, vol. 80, no. 10, pp. 1559-1589, Oct./1992.

[15] 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, pp. 2525-2565, Feb. 2011.

[16] W. Rucker, R. Hoschek, and K. Richter, "Various BEM formulations for calculating eddy currents in terms of field variables," *IEEE Transactions on Magnetics*, vol. 31, pp. 1336-1341, May 1995.

[17] K. Niino and N. Nishimura, "Calderón preconditioning approaches for PMCHWT formulations for Maxwell's equations," *International Journal of Numerical Modelling: Electronic Networks, Devices and Fields*, vol. 25, pp. 558-572, Sept. 2012.

[18] S. Yan, J. Jin, and Z. Nie, "A comparative study of Calderón preconditioners for PMCHWT equations," *IEEE Transactions on Antennas and Propagation*, vol. 58, pp. 2375-2383, July 2010.

[19] T. Chhim, A. Merlini, L. Rahmouni, J. Ortiz G., and F. 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*, vol. 1, pp. 534-548, 2020.

[20] F. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, "On a well-conditioned electric field integral operator for multiply connected geometries," *IEEE Transactions on Antennas and Propagation*, vol. 61, pp. 2077-2087, Apr. 2013.

[21] F. P. Andriulli, "Loop-star and loop-tree decompositions: Analysis and efficient algorithms," *IEEE Transactions on Antennas and Propagation*, vol. 60, pp. 2347-2356, May 2012.

[22] W. Hackbusch, "A sparse matrix arithmetic based on $H$-matrices. part i: Introduction to $H$-matrices,"

[23] W. Hackbusch and B. N. Khoromskij, "A sparse $H$-matrix arithmetic: Part II: Application to multi-dimensional problems," *Computing*, vol. 64, pp. 21-47, Feb. 2000.

[24] I. Bogaert, K. Cools, F. P. Andriulli, and D. De Zutter, "Low frequency scaling of the mixed MFIE for scatterers with a non-simply connected surface," in *2011 International Conference on Electromagnetics in Advanced Applications*, (Torino, Italy), pp. 951-954, IEEE, Sept. 2011.

[25] Z. Qian, W. Chew, and R. Suaya, "Generalized impedance boundary condition for conductor modeling in surface integral equation," *IEEE Transactions on Microwave Theory and Techniques*, vol. 55, pp. 2354-2364, Nov. 2007.

[26] A. Quarteroni and A. Valli, *Numerical Approximation of Partial Differential Equations*, vol. 23 of Springer Series in Computational Mathematics. Berlin, Heidelberg: Springer Berlin Heidelberg, 2008.

[27] H. V. Henderson and S. R. Searle, "On deriving the inverse of a sum of matrices," *SIAM Review*, vol. 23, pp. 53-60, Jan. 1981.

[28] L. Tsai, "A numerical solution for the near and far fields of an annular ring of magnetic current," *IEEE Transactions on Antennas and Propagation*, vol. 20, pp. 569-576, Sept. 1972.

[29] B. Popović, M. B. Dragović, and A. R. Djordjević, *Analysis and Synthesis of Wire-Antennas*. No. 2 in Electronic & Electrical Engineering Research Studies, Chichester, Sussex; New York: Research Studies Press, 1982.

[30] R. P. Feynman, *The Feynman Lectures on Physics. Volume 1: Mainly Mechanics, Radiation, and Heat*. New York: Basic Books, the new millennium edition, paperback first published ed., 2011.

[31] C. R. Paul, *Inductance: Loop and Partial*. Hoboken, NJ: Wiley, 2010.

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

**Viviana Giunzioni** (Member, IEEE) received the Laurea degree in Electronic and Communications Engineering in 2019, the Master of Science degree in Electronic Engineering in 2021, and the Ph.D. degree in Electrical, Electronics, and Communications Engineering in 2025, all from the Politecnico di Torino, Turin, Italy. Since 2021, she has been a Research Associate at the Politecnico di Torino. Her research interests include computational electromagnetics, focusing on preconditioning techniques and fast solvers for integral equations. Dr. Giunzioni received an Honorable Mention at the IEEE Antennas and Propagation Symposium (IEEE AP-S) in 2022 and another Honorable Mention at the URSI General Assembly and Scientific Symposium (URSI GASS) in 2023. In 2024, she was a finalist for the Young Scientist Award at the International Conference on Electromagnetics in Advanced Applications (ICEAA-IEEE APWC). In 2025, she received the Young Scientist Award at the URSI International Symposium on Electromagnetic Theory (EMTS).

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

**Alberto Scazzola** (Member, IEEE) received the Master Degree in Physics from the University of Trieste in 2018, and the PhD degree in Theory and Numerical Simulation of Condensed Matter from SISSA in 2022. His doctoral work involved research on high-temperature superconductivity. From 2022 to 2024, he held a Research Fellowship from the CERL team of the Politecnico di Torino. He is currently working as a high-school teacher. His research interests focus on computational electromagnetics and computational physics, involving the development of libraries and the implementation of numerical methods.

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

**Adrien Merlini** (S'16-M'19-SM'23) received the M.Sc. Eng. degree from the École Nationale Supérieure des Télécommunications de Bretagne (Télécom Bretagne), France, in 2015 and received the Ph.D. degree from the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique), France, in 2019.

From 2018 to 2019, he was a visiting Ph.D. student at 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 to brain imaging.

Dr. Merlini received Young Scientist Awards at the URSI GASS 2020 and the EMTS 2023 meetings. He authored and co-authored award-winning papers, including the third place at the EMTS 2023 best paper competition, the Best Paper Award at the 2022 ICEAA-IEEE APWC conference, five papers that received honorable mentions (URSI/IEEE-APS 2021, 2022, and 2023), and 3 best paper finalists (URSI GASS 2020, URSI/IEEE-APS 2021 and 2022). He is a member of IEEE-HKN, the IEEE Antennas and Propagation Society, URSI France, and of the Lab-STICC laboratory. He is currently serving as associate editor and editor for the Open Source Column for the Antenna and Propagation Magazine.