On a Low-Frequency and Contrast Stabilized Full-Wave Volume Integral Equation Solver for Lossy Media
Clément Henry, Student Member, IEEE, Adrien Merlini, Member, IEEE, Lyes Rahmouni, Member, IEEE,
and Francesco P. Andriulli,  Senior Member, IEEE

In this paper we present a new regularized electric flux volume integral equation (D-VIE) for modeling high-contrast conductive dielectric objects in a broad frequency range. This new formulation is particularly suitable for modeling biological tissues at low frequencies, as it is required by brain epileptogenic area imaging, but also at higher ones, as it is required by several applications including, but not limited to, transcranial magnetic and deep brain stimulation (TMS and DBS, respectively). When modeling inhomogeneous objects with high complex permittivities at low frequencies, the traditional D-VIE is ill-conditioned and suffers from numerical instabilities that result in slower convergence and in less accurate solutions. In this work we address these shortcomings by leveraging a new set of volume quasi-Helmholtz projectors. Their scaling by the material permittivity matrix allows for the re-balancing of the equation when applied to inhomogeneous scatterers and thereby makes the proposed method accurate and stable even for high complex permittivity objects until arbitrarily low frequencies. Numerical results, canonical and realistic, corroborate the theory and confirm the stability and the accuracy of this new method both in the quasi-static regime and at higher frequencies.

Volume integral equations, preconditioning, high-contrast (HC) conductive media, bio-electromagnetism.

§ INTRODUCTION

The electromagnetic modeling of human tissues has numerous applications that include brain source localization <cit.>, dosimetry <cit.>, deep brain stimulation (DBS) <cit.>, transcranial magnetic stimulation (TMS) <cit.>, electric impedance tomography <cit.>, and hyperthermic cancer therapy <cit.>. All these procedures require an accurate modeling of the interactions between electromagnetic fields and the human body. Depending on the application, numerical solvers can be either full-wave solvers directly derived from Maxwell's equations or static solvers based on Poisson's equation, which are only valid in the quasi-static regime <cit.>. While static solvers are sufficient to model resistive effects occurring in biological tissues, full-wave solvers are required when capacitive, inductive, and propagation effects should also be taken into consideration <cit.>. For instance, full-wave models are needed when the source has a larger spectral content (e.g. the magnetic pulse emitted by TMS coils <cit.> or the electrical current injected by an electrode in neurostimulation <cit.>). For these applications, the modeling should be done at various frequencies and, hence, the solver used should be able to perform accurately at arbitrary frequency.

Two main families of numerical solvers are widely employed for frequency domain bio-electromagnetic modeling: integral equation (IE) and differential equation solvers. In particular, integral equation solvers can be either surfacic (SIE) or volumic (VIE) depending on the nature of the tissue to model <cit.>. Although these solvers give rise to dense matrices, the higher computational cost incurred can be significantly reduced using acceleration techniques such as the fast multipole method (FMM) <cit.> or the adaptive cross approximation (ACA) <cit.>. Moreover, since IE formulations automatically enforce radiation conditions, no discretization is required outside the object. Unfortunately, however, bio-electromagnetic modeling is a challenging task for solvers based on integral equations because of the high complex permittivity contrast between the different tissues and their background <cit.>. Indeed, both a high-contrast object and/or a low operating frequency introduce a severe ill-conditioning in the discretized IE operator, which yields a slower convergence and a loss of accuracy in the solution <cit.>. These two phenomena are often referred to as the low-frequency (LF) and the high-contrast (HC) breakdowns respectively <cit.>.
In addition to the conditioning issue, the ill-scaling between the different components of the solution of the discretized system causes the latter to have fewer digits of accuracy due to finite machine precision <cit.>.

The high-contrast breakdown in piecewise homogeneous scatterers has been cured in the case of the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) formulation <cit.> and in a novel single-source integral equation <cit.> by leveraging the Calderón identities while the low-frequency breakdown of the PMCHWT was tackled by its preconditioning with quasi-Helmholtz projectors <cit.>. While these stabilized surface formulations have numerous advantages, they are limited to piecewise homogeneous models of biological tissues. Volume integral equations, instead, can model objects with a high degree of inhomogeneity. Unfortunately, as their surface counterparts, they suffer from the HC breakdown <cit.> and fail to converge rapidly in applications with high-permittivity contrast scatterers. Another limitation of traditional VIE is that, even though they are immune from the low-frequency breakdown in purely dielectric objects <cit.>, a frequency ill-scaling between the different parts of the VIE can occur when the object under study has a complex permittivity which depends on the frequency <cit.>. Therefore, in these cases, the LF breakdown can be considered an intrinsic part of the HC breakdown in the VIE. These limitations prevent the standard volume formulations to perform well in realistic biomedical applications where the modeling of high-contrast conductive tissues from the quasi-static regime to the microwave regime is required.

Regularization techniques have been introduced for solving the LF and the HC problems in the electric current VIE (J-VIE), the electric field VIE (E-VIE), and the electric flux VIE (D-VIE). The LF breakdown has been cured in the J-VIE for anisotropic and inhomogeneous scatterers using a Loop-Star-Facet decomposition for re-scaling properly the unknown <cit.>. However, solving this problem using a Helmholtz discretization deteriorates the dense discretization behavior of the VIE, i.e., it causes the conditioning of the system matrix to deteriorate when the average edge length $h$ of the discretized geometry decreases <cit.>. Instead, the quasi-Helmholtz projectors <cit.> allow for the removal of the ill-scaling in the formulation while keeping its dense mesh behavior unchanged. These projectors have been adapted to the J-VIE and used for curing the HC limitations of this equation for isotropic inhomogeneous scatterers <cit.>. Another approach to solve the HC problem is presented in <cit.>, where the E-VIE is regularized using symbol calculus and the Calderón identities. Its application to the J-VIE is discussed in <cit.>. While they are free from the HC breakdown, the two above-mentioned methods do not consider the numerical stability of the J-VIE or the E-VIE at low frequencies when modeling lossy dielectric objects, which is an important feature for a solver operating in low-frequency biomedical applications.
Regarding the D-VIE, an effective solution to both the LF and HC breakdowns has been proposed in <cit.> where an additional surface integral equation is used to adjust the background permittivity and lower the dielectric contrast. This technique is quite effective when the main problem is the contrast between the background and the object but it does not remove directly the internal contrasts between the different media within the simulated object, and an ill-conditioning thus remains for inhomogeneous objects.

The contribution of this paper is twofold, on the one hand we propose a new set of quasi-Helmholtz projectors that, differently from those proposed in the past, are the first of their kind to be oblique and to be particularly suited to manipulate the solenoidal and non-solenoidal parts the D-VIE. On the other hand, we leverage these new projectors to obtain a new regularized D-VIE which is immune from both ill-conditioning and the loss of accuracy occurring at low frequencies and in high-contrast objects.

More specifically, the new projectors, built from a weighted graph Laplacian matrix, when combined with the appropriate re-scaling, allow for the re-balancing the D-VIE in both the LF and HC regimes. The regularized D-VIE obtained is free from both the HC and the LF breakdowns and exhibits a solution that is valid until arbitrarily low frequencies, unlike standard full-wave solvers. This versatility makes it an appropriate formulation for biomedical applications where solvers that can operate in a broad frequency range and in high-contrast objects are required. The reader should however note that, although in our numerical validations we focus on biomedical applications and brain modeling in particular, the solver we propose is a completely general purpose one and several other application scenarios such as the modeling of lossy interconnects in printed circuit boards <cit.> could benefit from its use. Very preliminary results of this work have been presented in the conference contribution <cit.>.

This paper is organized as follows: in <Ref> we set the background and notation, including the definition of the D-VIE along with its discrete Helmholtz decomposition. The low-frequency and high-contrast behaviors of the D-VIE are analyzed in <Ref>. These analyses are followed by the presentation of a new set of quasi-Helmholtz projectors in <Ref> which is then employed in <Ref> to regularize the D-VIE at low frequencies and for high-contrast lossy scatterers. Bounds for the condition number of this new regularized formulation are derived in <Ref>. <Ref> is dedicated to the computational considerations and the introduction of another effective scheme for regularizing the D-VIE and <Ref> presents numerical examples demonstrating the stability and accuracy of these new formulations in a broad frequency range and for high-contrast media. Finally, for the sake of readability, we have omitted some of the mathematical technicalities in the main text. However, the interested reader will find them in the Appendices.

§ NOTATION AND BACKGROUND

Let $\Omega \subset \mathbb{R}^3$ be a simply connected object composed of a lossy dielectric and illuminated by a time-harmonic incident electric field $\vt{E}_i$ in a background of permittivity $\epsilon_0$ and permeability $\mu_0$. The scatterer is characterized by its complex isotropic relative permittivity $\epsilon_r(\vt{r}) = \epsilon_r'(\vt{r}) - \mi {\sigma(\vt{r})}/{\omega \epsilon_0}$ where $\vr \in \Omega$, $\epsilon_r'(\vr)$ is the relative permittivity, $\sigma(\vr)$ is the conductivity, and $\omega$ is the angular frequency of $\vt{E}_i$. The permeability of the scatterer is further assumed to be the permeability of vacuum $\mu_0$. Leveraging the volume equivalence principle, the scatterer can be substituted by a volume current density distribution

Convention used to define the SWG and the star functions: the SWG basis function is defined on the two tetrahedra $T_n^+$ and $T_n^-$ that are formed with their common face $\Gamma_n$ and the vertices $\vr_n^+$ and $\vr_n^-$, respectively. The half SWG basis functions are defined on $T_n^+$ only.

$ \vt{J}(\vt{r}) = \mi \omega \kappa(\vt{r}) \vt{D}(\vt{r})$
where $\kappa(\vt{r}) = ({\epsilon}(\vt{r}) - \epsilon_0)/{{\epsilon}(\vt{r})}$ is the dielectric contrast and $\vt{D}$ is the electric flux density. The D-VIE which relates $\vt{E}_i$ and $\vt{D}$ is expressed as <cit.>

\begin{equation}\label{eq:VIE}
\frac{\vt{D}(\vr)}{\epsilon(\vr)} - \frac{k_0^2}{\epsilon_0} \left(\operator{T}_A \vt{D}\right)(\vr) - \frac{1}{\epsilon_0} \left(\operator{T}_\Phi \vt{D}\right)(\vr) = \vt{E}_i(\vr)\,,\quad \vr \in \Omega,
\end{equation}

where the vector potential $\operator{T}_A$ and the scalar potential $\operator{T}_\Phi$ are defined as

\begin{gather}
\left(\operator{T}_A \vt{D}\right)(\vr) = \int_\Omega G_0(\vr,\vr') \kappa(\vr') \vt{D}(\vr') \diff v' \label{eq:operators_DVIE_Ta}\,,\\
\left(\operator{T}_\Phi \vt{D}\right)(\vr) = \nabla \int_\Omega G_0(\vr,\vr') \nabla' \cdot \left( \kappa(\vr') \vt{D}(\vr') \right) \diff v'\,, \label{eq:operators_DVIE_Tphi}
\end{gather}

in which $G_0(\vr,\vr') = \exp{(-\mi k_0 |\vr - \vr'|)}/|\vr - \vr'|$ is the free space 3D Green's function and $k_0 = \omega / c_0$ is the wavenumber of $\vt{E}_i$ in free space.

In general, the electric field radiated by an electric current distribution $\vt{J}_i$, which is commonly used to model the brain's electric activity, is

\begin{equation} \label{eq:dipole_field}
\vt{E}_i(\vr) = \left(\operator{T}_\Phi' \vt{J}_i\right)(\vr) + \left(\operator{T}_A' \vt{J}_i\right)(\vr)\,,
\end{equation}

where $\operator{T}_A' \vt{J}_i = -\mi k_0 \eta_0 \int_\Omega G_0(\vr,\vr') \vt{J}_i(\vr') \diff v'$ and $\operator{T}_\Phi' \vt{J}_i = {\eta_0}/{\mi k_0} \nabla \int_\Omega G_0(\vr,\vr') \nabla' \cdot \vt{J}_i(\vr') \diff v'$.

Equation (<ref>) is numerically solved by applying a Galerkin approach on a tetrahedral discretization of the geometry. Because the unknown $\vt{D}$ must exhibit a continuous normal component through material discontinuities, it is discretized with SWG basis functions $\{\vt{f_n}\}$ <cit.>. The SWG function associated to the face $\Gamma_n$ is defined as

\begin{equation} \label{eq:refSWGbg}
\vt{f_n}(\vr) =
\begin{cases}
\phantom{-}\frac{1}{3 V_n^+} (\vr - \vt{r_{n}^+}), \quad \vr \in T_n^+ \\
-\frac{1}{3 V_n^-} (\vr -\vt{r_{n}^-}), \quad \vr \in T_n^-
\end{cases}\,,
\end{equation}

where we have used the notations of Fig. <ref> and where $V_n^\pm$ is the volume of the tetrahedron $T_n^\pm$.

AM: add the definition of the SWGs?
AM: @ clement please check the modifs above

Using these functions, the unknown electric flux can be expanded as
\vt{D}(\vt{r}) \approx \sum_{n=1}^{N_F} [\mat{\alpha}]_n \vt{f_n}(\vt{r})

in which $N_F$ is the number of faces in the discretized geometry. The set of source functions $\{\vt{f_n}\}$ also contains the half SWG basis functions (supported by a single tetrahedron) that are defined at the boundary of the object to properly model the surface charges. The resulting equation, tested with the functions $\{ \vt{f_m}\}$ yields the linear system

\begin{gather}
\mat{Z} \mat{\alpha} = \left( \mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_\Phi} \right) \mat{\alpha} = \mat{v}\,, \text{ where}
\label{eq:discretized_eq} \\
\left[\mat{G_{\epsilon}}\right]_{mn} = \langle \vt{f_m},\vt{f_n} / \epsilon \rangle_\Omega \label{eq:def_gram}\,,\\
\left[ \mat{Z_A} \right]_{mn} = - {k_0^2}{\epsilon_0}^{-1} \langle \vt{f_m},\operator{T}_A \vt{f_n} \rangle_\Omega \label{def_vec_pot}\,,\\
\left[ \mat{Z_\Phi} \right]_{mn} = - {\epsilon_0}^{-1} \langle \vt{f_m},\operator{T}_\Phi \vt{f_n} \rangle_\Omega\,,
\end{gather}

and $\left[ \mat{v} \right]_{m} = \langle \vt{f_m},\vt{E}_i \rangle_\Omega$ with $\langle \vt{a} , \vt{b} \rangle_\Omega = \int_\Omega \vt{a} \cdot \vt{b} \diff v$.

In the following, for the sake of completeness, we briefly review some results on the discrete Helmholtz decomposition of the electric flux in the D-VIE. After being expanded with divergence conforming functions, the electric flux can be further decomposed into a linear combination of solenoidal (i.e. divergence free) and non-solenoidal (i.e. non-divergence free) components <cit.>. This decomposition is crucial because the solenoidal and non-solenoidal components of the D-VIE behave differently with respect to the frequency or the permittivity, hence they will have to be separated in the analysis and the regularization of the D-VIE. Such a separation can be performed using a Loop-Star decomposition on the discretized D-VIE, similarly to what is done in surface formulations <cit.>.

AM : explain what loops and stars are. Also, the notion of scaling of loop star is not clear and not defined.
AM: @clement check modifs + what do you mean decomposing the solution? You need to decompose the whole system

The unknown $\mat{\alpha}$ in (<ref>) can thus be decomposed into a sum of its Loop (solenoidal) and Star (non-solenoidal) components such that

\mat{\alpha} = \mat{\Lambda} \mat{l} + \mat{\Sigma} \mat{s}
$ <cit.>
where $\mat{\Lambda} \in \R^{N_F \times N_L}$ is the loop-to-SWG transformation matrix, $\mat{\Sigma} \in \R^{N_F \times N_S}$ is the star-to-SWG transformation matrix, $\mat{l}$ are the expansion coefficients of the unknown in the solenoidal basis, and $\mat{s}$ are the expansion coefficients in the non-solenoidal basis. The dimensions $N_L$ and $N_S$ are the numbers of independent loops and stars in the discretized geometry, respectively. In a tetrahedral discretization of a simply connected object, the number of stars and loops are <cit.>

\begin{equation}
N_S = N_T + N_{eF} ;~
N_L = N_{iE}- N_{iV}\,,
\end{equation}

where $N_T$, $N_{eF}$, $N_{iE}$, and $N_{iV}$ are the number of tetrahedra, external faces, internal edges, and internal vertices, respectively.
The mapping $\mat{\Sigma_v}\in \R^{N_F \times N_T}$ relates the stars and the SWG functions as

\begin{equation}
\left[\mat{\Sigma_v}\right]_{mn} =
\begin{cases}
\phantom{-}1 & \text{ if tetrahedron } n \text{ is tetrahedron } T_m^+ \\
-1 & \text{ if tetrahedron } n \text{ is tetrahedron } T_m^-\\
\phantom{-}0 & \text{ otherwise,}
\end{cases}
\label{eq:def_stars}
\end{equation}

where $T_m^+$ and $T_m^-$ represent the two tetrahedra on which the basis function $\vt{f_m}$ is defined (see Fig. <ref>). For half basis functions, only the entries corresponding to tetrahedron $T_m^+$ are filled in $\mat{\Sigma_v}$. Since these specific basis functions also model surface charges at the boundary of the object, another transformation matrix $\mat{\Sigma_s}\in \R^{N_F \times N_{eF}}$ needs to be filled for these basis functions

\begin{equation}
\left[\mat{\Sigma_s}\right]_{mn} =
\begin{cases} -1 & \text{ if face } n \text{ is face } \Gamma_m \\
\phantom{-}0 & \text{ otherwise}
\end{cases}\,.
\label{eq:def_stars_hbf}
\end{equation}

Note here that, by convention, the half basis function $\vt{f_m}$ is always supported by tetrahedron $T_m^+$.
The transformation matrix from star-to-SWG is then defined as $\mat{\Sigma} = \left[ \mat{\Sigma_v}~\mat{\Sigma_s} \right]$. After eliminating one column from $\mat{\Sigma}$ (zero total charge in $\Omega$) we obtain the full-column-rank matrix $\mat{\tilde{\Sigma}}$. The Loop functions are defined on the edges of the mesh as linear combinations of SWG basis functions <cit.>. Although schemes for the creation of an independent set of loops in a tetrahedral mesh exist <cit.>, we do not build them explicitly here, they are only introduced for supporting the discussion. For the purpose of this paper, we will employ the properties $\mat{\Lambda}^\mathrm{T} \mat{\tilde{\Sigma}} = \mat{0}$ and $\mat{\tilde{\Sigma}}^\mathrm{T} \mat{\Lambda} = \mat{0}$, that $\mat{\Lambda}$ and $\mat{\tilde{\Sigma}}$ satisfy by construction (refer to <cit.> and references therein).

In this work, we use a novel approach for performing a quasi-Helmholtz decomposition of the D-VIE without being limited by this conditioning issue. It requires the creation of volume quasi-Helmholtz projectors associated to $\mat{\Lambda}$ and $\mat{\Sigma}$. The projectors $\mat{P^{\Sigma}_{\identity}}$ and $\mat{P^{\Lambda}_{\identity}}$ can respectively select the solenoidal and the non-solenoidal parts of the discretized operator. They are defined as

\begin{gather}
\mat{P^{\Sigma}_{\identity}} = \mat{\Sigma} (\mat{\Sigma}^\mathrm{T} \mat{\Sigma})^+ \mat{\Sigma}^\mathrm{T} \label{eq:proj_std_loop}\,,\\
\mat{P^{\Lambda}_{\identity}} = \identity - \mat{P^{\Sigma}_{\identity}} = \mat{\Lambda} (\mat{\Lambda}^\mathrm{T} \mat{\Lambda})^{-1} \mat{\Lambda}^\mathrm{T} \label{eq:proj_std_star},
\end{gather}

where $^+$ is the Moore-Penrose pseudo inverse. Note that, thanks the orthogonality between $\mat{\Sigma}$ and $\mat{\Lambda}$, both $\mat{P^{\Sigma}_{\identity}}$ and $\mat{P^{\Lambda}_{\identity}}$ can be constructed using only the transformation matrix $\mat{\Sigma}$.

§ ANALYSIS AND REGULARIZATION OF THE D-VIE (SYMMETRIC D-VIE: NOT TO REVIEW)

In this work, we use a novel approach for performing a quasi-Helmholtz decomposition of the D-VIE without being limited by the conditioning issue of the discrete Helmholtz decomposition. To this aim, we propose a more general form of the volume quasi-Helmholtz projectors. These new projectors, $\mat{P^{\Sigma}_{A^{-1}}}$ and $\mat{P^{\Lambda}_{A^{-1}}}$, can respectively select the non-solenoidal and the solenoidal parts of the discretized operator. They are defined as

\begin{gather}
\mat{P^{\Sigma}_{A^{-1}}} = \mat{A} \mat{\Sigma} \left( \mat{\Sigma}^\mathrm{T} \mat{A} \mat{\Sigma} \right)^{+} \mat{\Sigma}^\mathrm{T} \label{eq:scaled-projectors-star-A}\,,\\
\mat{P^{\Lambda}_{A^{-1}}} = \identity - \mat{P^{\Sigma}_{A^{-1}}} = \mat{\Lambda} \left( \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\Lambda} \right)^{-1} \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \label{eq:scaled-projectors-loop-A} ,
\end{gather}

where $^+$ is the Moore-Penrose pseudo inverse. Differently from the definitions of quasi-Helmholtz projectors in previous works <cit.>
[inline]add citation here. Please note that I have changed in (16) and (17) the name of the projectors because you were using your notation in a inconsistent and confusing way for G. See if you get a better idea because mine is not super...
we introduce here an invertible matrix scaling $\mat{A}$. Given this new definition, identity (<ref>) does not follow directly from the properties of standard projectors, but its proof is provided in Appendix <ref>. These new projectors are not necessarily symmetric and their transposes are denoted by $\mat{P^{\Lambda \mathrm{\normalfont{T}}}_{A^{-1}}}$ and $\mat{P^{\Sigma \mathrm{\normalfont{T}}}_{A^{-1}}}$ in the following. When $\mat{A} = \identity$, we obtain volume quasi-Helmholtz projectors, denoted by $\mat{P^\Lambda_{\identity}}$ and $\mat{P^\Sigma_{\identity}}$ here, that have the same form as their surface counterparts <cit.>. Instead, by choosing $\mat{A}= \mat{G_{\epsilon}}^{-1}$, we obtain the scaled projectors $\mat{P^{\Sigma}_G}$ and $\mat{P^{\Lambda}_G}$, for which it can be easily seen that the following identities hold

\begin{gather}
\mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{G_{\epsilon}} \mat{P^{\Lambda}_G} = \mat{G_{\epsilon}} \mat{P^{\Lambda}_G} \label{eq:id1_scaled_proj}\,,\\
\mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{G_{\epsilon}} \mat{P^{\Sigma}_G} = \mat{0} \label{eq:id2_scaled_proj}\,,\\
\mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \mat{G_{\epsilon}} \mat{P^{\Lambda}_G } = \mat{0} \label{eq:id3_scaled_proj}\,,\\
\mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \mat{G_{\epsilon}} \mat{P^{\Sigma}_G} =  \mat{G_{\epsilon}} \mat{P^{\Sigma}_G}
\label{eq:id4_scaled_proj}.
\end{gather}

We will further employ $\mat{P^{\Lambda}_{\identity}}$ and $\mat{P^{\Sigma}_{\identity}}$ ($\mat{A} = \identity$) for the frequency stabilization of the D-VIE and the analysis of its high-contrast behavior. However, to solve the high-contrast problem of the D-VIE in lossy dielectric objects with an inhomogeneous permittivity profile, we will leverage $\mat{P^{\Sigma}_G}$ and $\mat{P^{\Lambda}_G}$ ($\mat{A} = \mat{G_{\epsilon}}^{-1}$), which allow a re-scaling that depends on the position in the object.

§.§ Low-frequency behavior of the D-VIE

We analyze here the low-frequency behavior of the D-VIE when modeling lossy dielectric objects surrounded by free space. At low frequencies, VIEs are subject to the low-frequency breakdown due to the frequency dependence of the complex permittivity of the object <cit.>. In the discretized D-VIE, this ill-scaling causes the ill-conditioning of the system matrix and the loss of significant digits in the solution coefficients. In this section, we first expose the low-frequency behavior of the D-VIE through a Loop-Star decomposition. This behavior will then be regularized using the volume quasi-Helmholtz projectors $\mat{P^{\Lambda}_{\identity}}$ and $\mat{P^{\Sigma}_{\identity}}$ in combination with appropriate scalings determined during the Loop-Star analysis to regularize the D-VIE at low frequencies without deteriorating its dense discretization behavior.

First, the Loop-Star decomposition of the D-VIE is obtained using the mapping matrix $\mat{B_{\Lambda \Sigma}}$ as

\begin{equation} \label{eq:helmholtz_decomp}
\mat{B^\mathrm{\textnormal{T}}_{\Lambda \Sigma}}(\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi}}) \mat{B_{\Lambda \Sigma}} \mat{\alpha_{\Lambda \Sigma}} = \mat{B^\mathrm{\textnormal{T}}_{\Lambda \Sigma}} \mat{v},
\end{equation}

where $\mat{B_{\Lambda \Sigma}} \mat{\alpha_{\Lambda \Sigma}} = \mat{\alpha}$. Given that the application of $\mat{\Lambda}$ to $\mat{Z_{\Phi}}$ yields an exact cancellation of the volume parts of the scalar potential and of some of its surface parts, we further decompose $\mat{Z_{\Phi}}$ to make it possible to analyze the behavior of the remaining terms

\begin{equation}
\mat{Z_{\Phi}} = \mat{Z_{\Phi, 11}} + \mat{Z_{\Phi ,\epsilon \epsilon}} + \mat{Z_{\Phi,\epsilon 1}} + \mat{Z_{\Phi,1 \epsilon}},
\end{equation}

in which $\mat{Z_{\Phi, 11}}$ only accounts for the volumic contributions in the object and the surfacic contributions on its boundary ($\partial \Omega$)

\begin{equation} \label{eq:Z11}
\begin{split}
&[\mat{Z_{\Phi,11}}]_{nm} = \\
& \epsilon_0^{-1}\left[ \int_\Omega \kappa(\vr) \nabla \cdot \vt{f_n}(\vt{r}) \int_\Omega G_0(\vr,\vr') \kappa(\vr') \nabla \cdot \vt{f_m}(\vt{r'}) \diff v' \diff v \right. \\
&\left. - \int_\Omega \kappa(\vr) \nabla \cdot \vt{f_n}(\vt{r}) \int_{\partial \Omega} G_0(\vr,\vr') \kappa^+_m \uv{n}_m \cdot \vt{f_m}(\vt{r'}) \diff s' \diff v \right.\\
&\left. - \int_{\partial \Omega} \kappa^+_n \uv{n}_n \cdot \vt{f_n}(\vt{r}) \int_{\Omega} G_0(\vr,\vr') \kappa(\vr')\nabla \cdot \vt{f_m}(\vt{r'}) \diff v' \diff s \right.\\
&\left.+ \int_{\partial \Omega} \kappa^+_n \uv{n}_n \cdot \vt{f_n}(\vt{r}) \int_{\partial \Omega} G_0(\vr,\vr') \kappa^+_m \uv{n}_m \cdot \vt{f_m}(\vt{r'}) \diff s' \diff s \right],
\end{split}
\end{equation}

with $\kappa^+_n$ being the dielectric contrast in tetrahedron $T^+_n$ and $\uv{n}_n$ the outward unit normal of the triangle of $T^+_n$ which pertains to $\partial \Omega$. $\mat{Z_{\Phi, \epsilon, \epsilon}}$ includes the purely surfacic contributions internal to the object (at the interface $\Gamma_n$ between two tetrahedra $T^-_n$ and $T^+_n$ of different material contrasts $\kappa^-_n$ and $\kappa^+_n$)

\begin{equation} \label{eq:Zepseps}
\begin{split}
&[\mat{Z_{\Phi,\epsilon \epsilon}}]_{nm} = \epsilon_0^{-1} . \\
& \int_{\Gamma_n} \deltaKappa_n \uv{n}_n \cdot \vt{f_n}(\vt{r}) \int_{\Gamma_m} G_0(\vr,\vr') \deltaKappa_m \uv{n}_m \cdot \vt{f_m}(\vt{r'}) \diff s' \diff s,
\end{split}
\end{equation}

in which $\deltaKappa_n = \kappa^-_n - \kappa^+_n$ and $\uv{n}_m$ is a unit vector normal to $\Gamma_m$ oriented from $T^+_m$ to $T^-_m$. Finally, $\mat{Z_{\Phi, \epsilon 1}}$ and $\mat{Z_{\Phi, 1\epsilon}}$ account for the associated cross terms

\begin{equation} \label{eq:Z1eps}
\begin{split}
&[\mat{Z_{\Phi,1\epsilon}}]_{nm} = \epsilon_0^{-1} . \\
&\left[\int_\Omega \kappa(\vr) \nabla \cdot \vt{f_n}(\vt{r}) \int_{\Gamma_m} G_0(\vr,\vr') \deltaKappa_m \uv{n}_m \cdot \vt{f_m}(\vt{r'}) \diff s' \diff v \right.\\
&\left.- \int_{\partial \Omega} \kappa^+_n \uv{n}_n \cdot \vt{f_n}(\vt{r}) \int_{\Gamma_m} G_0(\vr,\vr') \deltaKappa_m \uv{n}_m \cdot \vt{f_m}(\vt{r'}) \diff s' \diff s \right],
\end{split}
\end{equation}

\begin{equation} \label{eq:Zeps1}
\begin{split}
&[\mat{Z_{\Phi,\epsilon 1}}]_{nm} = \epsilon_0^{-1} . \\
&\left[ \int_{\Gamma_n} \deltaKappa_n \uv{n}_n \cdot \vt{f_n}(\vt{r}) \int_{\Omega} G_0(\vr,\vr') \kappa(\vr') \nabla \cdot \vt{f_m}(\vt{r'}) \diff v' \diff s \right.\\
&\left.- \int_{\Gamma_n} \deltaKappa_n \uv{n}_n \cdot \vt{f_n}(\vt{r}) \int_{\partial \Omega} G_0(\vr,\vr') \kappa^+_m \uv{n}_m \cdot \vt{f_m}(\vt{r'}) \diff s' \diff s\right],
\end{split}
\end{equation}

respectively. The following cancellations can be shown to arise when applying $\mat{\Lambda}^\mathrm{T}$ and $\mat{\Lambda}$ to the scalar potential

\begin{equation}
\mat{\Lambda}^\mathrm{T} \mat{Z_{\Phi, 11}}= \mat{0}, ~ \mat{\Lambda}^\mathrm{T} \mat{Z_{\Phi, 1\epsilon}}= \mat{0}, ~ \mat{Z_{\Phi, 11}} \mat{\Lambda} = \mat{0}, ~ \mat{Z_{\Phi, \epsilon 1}} \mat{\Lambda} = \mat{0}.
\end{equation}

The fact that $\mat{\Lambda}$ cancels the surface terms defined on $\partial \Omega$ in $\mat{Z_{\Phi, 11}}$, $\mat{Z_{\Phi, 1\epsilon}}$, and $\mat{Z_{\Phi, \epsilon 1}}$ results from the absence of loops functions on $\partial \Omega$. This property is proven in Appendix <ref>.

We then represent $\mat{Z_{\Lambda \Sigma}}$ in (<ref>) as a 2-by-2 block matrix

\begin{equation} \label{eq:block_mat_loop-star}
\mat{Z_{\Lambda \Sigma}} = \mat{B^\mathrm{\textnormal{T}}_{\Lambda \Sigma}}\mat{Z} \mat{B_{\Lambda \Sigma}} = \begin{bmatrix}
\mat{\Lambda}^\mathrm{T} \mat{Z} \mat{\Lambda} & \mat{\Lambda}^\mathrm{T} \mat{Z} \mat{\Sigma}\,,\\
\mat{\Sigma}^\mathrm{T} \mat{Z} \mat{\Lambda} & \mat{\Sigma}^\mathrm{T} \mat{Z} \mat{\Sigma}
\end{bmatrix}, \\
\end{equation}

in which

\begin{gather}
\mat{\Lambda}^\mathrm{T} \mat{Z} \mat{\Lambda} = \mat{\Lambda}^\mathrm{T} (\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi ,\epsilon \epsilon}})\mat{\Lambda} \label{eq:decomp_LL}\,,\\
\mat{\Lambda}^\mathrm{T} \mat{Z} \mat{\Sigma} = \mat{\Lambda}^\mathrm{T}(\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi ,\epsilon \epsilon}} + \mat{Z_{\Phi,\epsilon 1}} ) \mat{\Sigma} \label{eq:decomp_LS}\,,\\
\mat{\Sigma}^\mathrm{T} \mat{Z} \mat{\Lambda} = \mat{\Sigma}^\mathrm{T}(\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi ,\epsilon \epsilon}} + \mat{Z_{\Phi,1 \epsilon}})\mat{\Lambda} \label{eq:decomp_SL}\,,\\
\mat{\Sigma}^\mathrm{T} \mat{Z} \mat{\Sigma} = \mat{\Sigma}^\mathrm{T}(\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi, 11}} + \mat{Z_{\Phi ,\epsilon \epsilon}} + \mat{Z_{\Phi,\epsilon 1}} + \mat{Z_{\Phi,1 \epsilon}})\mat{\Sigma}.
\label{eq:decomp_SS}
\end{gather}

We have the following frequency scalings for the real and imaginary parts of each term in the decomposition of $\mat{Z_{\Lambda \Sigma}}$ ((<ref>), (<ref>), (<ref>), and (<ref>))

\begin{gather}
\Re(\mat{G_{\epsilon}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{G_{\epsilon}}) \underrel{\omega \to 0}{=} \omega \label{eq:scaling-gram}\,,\\
\Re(\mat{Z_A}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{Z_A}) \underrel{\omega \to 0}{=} \omega^3\\
\Re(\mat{Z_{\Phi,11}}) \underrel{\omega \to 0}{=} 1, ~
\Im(\mat{Z_{\Phi,11}}) \underrel{\omega \to 0}{=} \omega\\
\Re(\mat{Z_{\Phi,\epsilon \epsilon}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{Z_{\Phi,\epsilon \epsilon}}) \underrel{\omega \to 0}{=} \omega^3\\
\Re(\mat{Z_{\Phi,\epsilon 1}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{Z_{\Phi,\epsilon 1}}) \underrel{\omega \to 0}{=} \omega\\
\Re(\mat{Z_{\Phi,1 \epsilon}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{Z_{\Phi,1 \epsilon}}) \underrel{\omega \to 0}{=} \omega \label{eq:scaling-1eps}.
\end{gather}

Using (<ref>) to (<ref>), we deduce the scalings for $\mat{Z_{\Lambda \Sigma}}$ at low frequencies

\begin{gather}
\Re\left(\mat{Z_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) &\mathcal{O}(\omega^2) \\
\mathcal{O}(\omega^2) & \mathcal{O}(1)
\end{bmatrix} \label{eq:scaling_mat_real}\,,\\
\Im\left( \mat{Z_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) &\mathcal{O}(\omega) \\
\mathcal{O}(\omega) & \mathcal{O}(\omega)
\end{bmatrix} \label{eq:scaling_mat_imag},
\end{gather}

which, according to the Gershgorin circle theorem, confirms that $\mat{Z_{\Lambda \Sigma}}$ is ill-conditioned. In the following we provide the scalings of the right hand side for plane wave and dipole excitations, which are also required in the numerical stability analysis. For a plane wave excitation, the right hand side of (<ref>) has the following scalings when $\omega\rightarrow 0$

\begin{equation} \label{eq:scaling_rhs}
\Re \left(\mat{B_{\Lambda \Sigma}^\mathrm{\textnormal{T}}} \mat{v_{\textnormal{{PW}}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\mat{B^\mathrm{\textnormal{T}}_{\Lambda \Sigma}} \mat{v_{\textnormal{{PW}}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix},
\end{equation}

in which $\left[ \mat{v_{\textnormal{PW}}} \right]_{n} = \int_\Omega \kappa(\vr) \vt{f_n}(\vr) \cdot \vt{E}_0 \exp{(-\mi \vt{k} \cdot \vr )} \diff v$ with $\vt{E}_0$ being the polarization of the plane wave and $\vt{k}$ its wave vector. In contrast to the electric field integral equation <cit.>, these scalings were obtained by considering the inhomogeneous conductivity profile of the object and the testing $\{\kappa \vt{f_n}\}$ of the D-VIE we opted for. To derive the frequency dependence of the Loop-Star decomposition of a dipole excitation, we start from the expression of the field it radiates (<ref>) in which $\vt{J}_i(\vr, \vr_0) = \mi \omega \delta(\vr - \vr_0)\vt{p}$ is the current dipole with $\delta$, $\vt{p}$, and $\vt{r_0}$ being the Dirac delta function, the dipole moment, and the dipole position, respectively. We then test (<ref>) with $\{\kappa \vt{f_n}\}$ to obtain the discretized right hand side

\begin{equation}
\begin{split}
\left[ \mat{v_{\textnormal{dip}}} \right]_{n} &= \left[ \mat{v_{\mathrm{A}}} \right]_{n} + \left[ \mat{v_{\mathrm{\Phi}}} \right]_{n}\\
&= \frac{k_0^2}{\epsilon_0} \vt{p} \cdot \int_\Omega \kappa(\vr) \vt{f_n}(\vr) G_0(\vr, \vr_0) \diff v\\
&-\frac{1}{ \epsilon_0}\int_{\Omega} \nabla \cdot \left(\kappa(\vr) \vt{f_n}(\vr)\right) \vt{p} \cdot \nabla G_0(\vr, \vr_0) \diff v.
\end{split}
\end{equation}

From the Taylor expansion of $G_0$, its gradient, and the inverse of the complex permittivity, we obtain that $\Re(G_0(\vr,\vr_0)) = \mathcal{O}(1)$, $\Im(G_0(\vr,\vr_0)) = \mathcal{O}(\omega)$, $\Re(\nabla G_0(\vr,\vr_0)) = \mathcal{O}(1)$, $\Im(\nabla G_0(\vr,\vr_0)) = \mathcal{O}(\omega^3)$, $\Re(\epsilon(\vr)^{-1}) = \mathcal{O}(\omega^2)$, and $\Im(\epsilon(\vr)^{-1}) = \mathcal{O}(\omega)$ at low frequencies. This yields the following frequency dependencies for the real and imaginary parts of the Loop-Star decomposition of $\mat{v_{\mathrm{A}}}$ and $\mat{v_{\mathrm{\Phi}}}$

\begin{gather}
\Re(\mat{\Lambda}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{\Lambda}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \omega^3\\
\Re(\mat{\Sigma}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{\Sigma}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \omega^3\\
\Re(\mat{\Lambda}^\mathrm{T} \mat{v_{\mathrm{\Phi}}}) \underrel{\omega \to 0}{=} \omega^2, ~
\Im(\mat{\Lambda}^\mathrm{T} \mat{v_{\mathrm{\Phi}}}) \underrel{\omega \to 0}{=} \omega\\
\Re(\mat{\Sigma}^\mathrm{T} \mat{v_{\mathrm{\Phi}}} ) \underrel{\omega \to 0}{=} 1, ~
\Im(\mat{\Sigma}^\mathrm{T} \mat{v_{\mathrm{\Phi}}} ) \underrel{\omega \to 0}{=} \omega,
\end{gather}

from which we obtain the scalings of $\mat{v_{\mathrm{dip}}}$ at low frequencies

\begin{equation} \label{eq:scaling_rhs_diple}
\Re \left(\mat{B_{\Lambda \Sigma}^\mathrm{\textnormal{T}}} \mat{v_{\textnormal{dip}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\mat{B^\mathrm{\textnormal{T}}_{\Lambda \Sigma}} \mat{v_{\textnormal{dip}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix}.
\end{equation}

[inline]CH: @Francesco: We obtain numerically and analytically the same scalings for the PW and the dipole, but I couln't find a intuitive reason for this. Shall will add something in this paragraph?

Note that since the scalings of the right hand side of the dipole are the same as for the plane wave, the following analysis will apply for both right hand sides. To make this clear we will use a unique right hand side denomination $\mat{v}$.

Subsequently, to retrieve the scalings of the solution, we need to derive the scalings of the inverse of $\mat{Z_{\Lambda \Sigma}}$ by inverting (<ref>) using Schur complement formulas <cit.>

\begin{gather}
\Re\left(\mat{Z_{\Lambda \Sigma}}^{-1} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) &\mathcal{O}(1) \\
\mathcal{O}(1) & \mathcal{O}(1)
\end{bmatrix} \label{eq:scaling_mat_inv_real}\,,\\
\Im\left( \mat{Z_{\Lambda \Sigma}}^{-1} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1/\omega) &\mathcal{O}(\omega) \\
\mathcal{O}(\omega) & \mathcal{O}(\omega)
\end{bmatrix} \label{eq:scaling_mat_inv_imag}.
\end{gather}

Finally by multiplying the scaling matrix of $\mat{Z_{\Lambda \Sigma}}^{-1}$ ((<ref>) and (<ref>)) and the scaling vector of $\mat{B^\mathrm{\normalfont{T}}_{\Lambda \Sigma}} \mat{v}$ (<ref>), we obtain the same scalings of $\mat{\alpha_{\Lambda \Sigma}}$ for dipole and plane wave excitations

\begin{equation} \label{eq:scaling_solution}
\Re \left(\mat{\alpha_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\mat{\alpha_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix}.
\end{equation}

From the scalings obtained above, we observe that the real part of the Loop component of the right hand side ((<ref>) and (<ref>)) scales as $\mathcal{O}(\omega^2)$ while the real part of the Star component scales as $\mathcal{O}(1)$, which means that the real part of the Loop component of the right hand side is lost at low frequencies, when these two components are not separated and stored together in finite precision arithmetic. By propagating the scalings of the right hand side to the solution leveraging (<ref>) and (<ref>), it results that the real part of the Loop component of the right hand side contributes to the dominant imaginary part of the Loop component of the solution $\mat{\alpha_{\Lambda \Sigma}}$ (<ref>), the right hand side must be re-scaled appropriately to retrieve correctly the electric field inside the object.

To solve the two above-mentioned issues, we regularize (<ref>) by re-scaling appropriately the loop/star-to-SWG transformation matrices, resulting in the following Loop-Star D-VIE

\begin{equation} \label{eq:helmholtz_decomp_precond}
\begin{bmatrix}
C_1 \mat{\Lambda}^\mathrm{T} \\
C_2 \mat{\Sigma}^\mathrm{T}
\end{bmatrix} (\mat{G_{\epsilon}}+ \mat{Z_A} + \mat{Z_{\Phi}}) [C_3 \mat{\Lambda}~C_4\mat{\Sigma}] \mat{\alpha_{\Lambda \Sigma}} = \begin{bmatrix}
C_1 \mat{\Lambda}^\mathrm{T} \\
C_2 \mat{\Sigma}^\mathrm{T}
\end{bmatrix} \mat{v},
\end{equation}

in which the coefficients $C_1$, $C_2$, $C_3$ and $C_4$ have to be determined to make the Loop and Star parts of the preconditioned matrix ($\mat{\tilde{Z}_{\Lambda \Sigma}}$), its right hand side ($\mat{\tilde{B}^\mathrm{\textnormal{T}}_{\Lambda \Sigma}} \mat{v}$), and its solution ($\mat{\alpha_{\Lambda \Sigma}}$) free from loss of accuracy. Since the loop and Star components of the solution of (<ref>) have the same frequency dependence for both the real and the imaginary parts (<ref>), we can set $C_3$ and $C_4$ to $1$ and attempt to leverage the two remaining coefficients to regularize the rest of the equation. Therefore, we now need to determine $C_1$ and $C_2$ for $\mat{Z_{\Lambda \Sigma}}$ to be well-conditioned and to avoid numerical loss. By choosing $C_1 = \mi / \omega $ and $C_2 = 1$, the matrix becomes well-conditioned and the right hand side is no longer subject to numerical loss of accuracy for both plane wave and dipole excitations. The resulting Loop-Star matrix $\mat{\tilde{Z}_{\Lambda \Sigma}} = \mat{\tilde{B}^\mathrm{\normalfont{T}}_{\Lambda \Sigma}} (\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi}}) \mat{B_{\Lambda \Sigma}}$ has the following frequency scalings

\begin{gather} \label{eq:scaling_mat_precond}
\Re\left(\mat{\tilde{Z}_{\Lambda \Sigma}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) &\mathcal{O}(1) \\
\mathcal{O}(\omega^2) & \mathcal{O}(1)
\end{bmatrix}\,,\\
\Im\left(\mat{\tilde{Z}_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) &\mathcal{O}(\omega) \\
\mathcal{O}(\omega) & \mathcal{O}(\omega)
\end{bmatrix},
\end{gather}

where $\mat{\tilde{B}_{\Lambda \Sigma}} = [C_1 \mat{\Lambda}~C_2 \mat{\Sigma}]$. The preconditioned right hand side (plane wave or dipole) scales as

\begin{equation} \label{eq:scaling_rhs_precond_dipole}
\Re \left(\mat{\tilde{B}^{\textnormal{T}}_{\Lambda \Sigma}} \mat{v}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\mat{\tilde{B}^{\textnormal{T}}_{\Lambda \Sigma}} \mat{v}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix},
\end{equation}

when $\omega$ goes to $0$. Therefore, the coefficients $C_1$, $C_2$, $C_3$, and $C_4$ together with the Loop-Star decomposition have successfully cured the ill-scaling in the matrix and the right hand side.

Subsequently, we can logically relate the above-mentioned coefficients and the quasi-Helmholtz projectors defined in (<ref>) and (<ref>) with $\mat{A} = \identity$ to build a preconditioned equation that does not suffer from the dense discretization breakdown. We obtain the following left and right preconditioners

\begin{gather}
\mat{L} = C_1 \mat{P^{\Lambda}_{\identity}} + C_2 \mat{P^\Sigma_{\identity}} \label{eq:preconditioner_left_definition}\,,\\
\mat{R} = C_3 \mat{P^\Lambda_{\identity}} + C_4 \mat{P^\Sigma_{\identity}}.
\label{eq:preconditioner_right_definition}
\end{gather}

Since $C_3=C_4=1$, we have $\mat{R} = \identity$, and the resulting preconditioned equation to be solved reads

\begin{equation} \label{eq:regularized_sys_proj_lf}
\mat{L}^\mathrm{T} \mat{Z} \mat{\alpha} = \mat{L}^\mathrm{T} \mat{v},
\end{equation}

in which

\begin{equation}
\begin{split}
\mat{L} &= C_1 \gamma_\Lambda\mat{P^{\Lambda}_{\identity}} + C_2 \gamma_\Sigma \mat{P^\Sigma_{\identity}} \\
&= \frac{\mi \gamma_\Lambda}{ \omega} \mat{P^{\Lambda}_{\identity}} + \gamma_\Sigma \mat{P^\Sigma_{\identity}}. \label{eq:preconditioner_left}
\end{split}
\end{equation}

Note that the frequency-independent coefficients $\gamma_\Lambda$ and $\gamma_\Sigma$ were added in (<ref>) to make the norms of $\left(\mi \gamma_\Lambda/\omega\right) \mat{P^{\Lambda }_{\identity}} \mat{Z} \mat{P^{\Lambda}_{\identity}}$ and of $ \gamma_\Sigma \mat{P^{\Sigma}_{\identity}} \mat{Z} \mat{P^{\Sigma}_{\identity}}$ of order $\mathcal{O}(1)$ regardless of the geometry. They are defined as

\begin{gather}
\gamma_\Lambda = ({\epsilon_0 \| \mat{P^{\Lambda}_{\identity}} \mat{G} \mat{P^{\Lambda}_{\identity}} \|})^{-1} \label{eq:coeff_lambda_1}\,,\\
\gamma_\Sigma = ({\| \mat{P^{\Sigma}_{\identity}} \mat{Z_\Phi} \mat{P^{\Sigma}_{\identity}} \|})^{-1},
\label{eq:coeff_sigma_1}
\end{gather}

with $\norm {\cdot} $ denoting the $l^2\text{-norm}$ and $\left[\mat{G}\right]_{nm} = \langle \vt{f_n},\vt{f_m} \rangle_\Omega$ being the SWG-SWG Gram matrix.

§.§ High-contrast behavior of the D-VIE

In the following we analyze the high-contrast behavior of the D-VIE. Note that, since the frequency ill-scaling originates from the frequency dependence of the complex permittivity, the preconditioner derived in this section will also solve the LF breakdown, which is an intrinsic part of the HC breakdown. This section is complementary with the previous one in the sense that we treat the ill-scaling of the D-VIE due to the inhomogeneity inside the object, which requires a re-balancing of the equation with a scaling that is a function of the position in the object. We address at first the HC problem of the D-VIE in the quasi-static regime ($k_0 D\ll 1$ with $D$ being the diameter of the object), in which the vector potential matrix $\mat{Z_A}$ can be neglected. The discretized D-VIE operator without the vector potential matrix is expressed as

\begin{equation}
\mat{Z} = \mat{G_{\epsilon}} + \mat{Z_\Phi}.
\label{eq:discretized_eq_noVec}
\end{equation}

To identify the scalings of (<ref>) with HC objects in which $\epsilon_r(\vr) \rightarrow \infty$ and $\kappa(\vr) \rightarrow 1$, we use the volume quasi-Helmholtz projectors $\mat{P^{\Lambda}_{\identity}}$ and $\mat{P^{\Sigma}_{\identity}}$ to decompose $\mat{Z}$ into four matrices representing the solenoidal and non-solenoidal parts of $\mat{Z}$ and the associated cross terms

\begin{equation}
\begin{split}
\mat{Z} &= \left( \mat{P^{\Sigma}_{\identity}} + \mat{P^{\Lambda}_{\identity}} \right)^\mathrm{T} \mat{Z} \left( \mat{P^{\Sigma}_{\identity}} + \mat{P^{\Lambda}_{\identity}} \right) \\
&= \mat{P^{\Sigma}_{\identity}} \mat{Z} \mat{P^{\Sigma}_{\identity}} + \mat{P^{\Sigma}_{\identity}} \mat{Z} \mat{P^{\Lambda}_{\identity}}\\
\mat{P^{\Lambda}_{\identity}} \mat{Z} \mat{P^{\Sigma}_{\identity}} + \mat{P^{\Lambda}_{\identity}} \mat{Z} \mat{P^{\Lambda}_{\identity}},
\label{eq:discretized_eq_decomp}
\end{split}
\end{equation}

where we used the fact that $\mat{P^{\Sigma}_{\identity}}$ and $\mat{P^{\Lambda}_{\identity}}$ are symmetric projectors.

We note that this decomposition follows the same approach as the decomposition performed in (<ref>) of the previous section, except that we now employ projectors instead of the loop/star-to-SWG tranformation matrix. By taking the limit of each term in (<ref>) when the relative permittivity goes to infinity and using the fact that $\mat{G_{\epsilon}}= \mathcal{O}(\epsilon_r(\vr)^{-1})$, $\mat{Z_\Phi} = \mathcal{O}(1)$, $\mat{P^{\Lambda}_{\identity}} \mat{Z_\Phi} = \mathcal{O}(\epsilon_r(\vr)^{-1})$, and $ \mat{Z_\Phi} \mat{P^{\Lambda}_{\identity}} = \mathcal{O}(\epsilon_r(\vr)^{-1})$ (surface contributions in $\mat{Z^\Phi}$) when $\epsilon_r(\vr) \rightarrow \infty$, we get the following scalings

\begin{gather}
\mat{P^{\Lambda}_{\identity}} \mat{Z} \mat{P^{\Lambda}_{\identity}} = \mathcal{O}(\epsilon_r(\vr)^{-1})\quad \epsilon_r(\vr) \rightarrow \infty, \label{eq:p_block_LL}\,,\\
\mat{P^{\Lambda}_{\identity}} \mat{Z} \mat{P^{\Sigma}_{\identity}} =\mathcal{O}(\epsilon_r(\vr)^{-1})\quad \epsilon_r(\vr) \rightarrow \infty, \label{eq:p_block_LS}\,,\\
\mat{P^{\Sigma}_{\identity}} \mat{Z} \mat{P^{\Lambda}_{\identity}} =\mathcal{O}(\epsilon_r(\vr)^{-1})\quad \epsilon_r(\vr) \rightarrow \infty, \label{eq:p_block_SL}\,,\\
\mat{P^{\Sigma}_{\identity}} \mat{Z} \mat{P^{\Sigma}_{\identity}} =\mathcal{O}(1)\quad \epsilon_r(\vr) \rightarrow \infty .
\label{eq:p_block_SS}\end{gather}

From the scalings obtained above and using the Gershgorin circle theorem, we can conclude that the standard D-VIE is also ill-scaled for high permittivities. As for the LF breakdown, the solenoidal part of the system matrix should be re-scaled independently of the rest of the equation. However, due to the inhomogeneity in the object, we need to use another form of the volume quasi-Helmholtz projectors defined in (<ref>) and (<ref>) to properly re-scale the equation. Here we set $\mat{A}= \mat{G_{\epsilon}}^{-1}$ in (<ref>) and (<ref>) and obtain the following quasi-Helmholtz projectors

\begin{gather}
\mat{P^{\Sigma}_G} = \mat{G_{\epsilon}}^{-1}\mat{\Sigma} \left( \mat{\Sigma}^\mathrm{T} \mat{G_{\epsilon}}^{-1} \mat{\Sigma} \right)^{+} \mat{\Sigma}^\mathrm{T} \label{eq:scaled-projectors-star-G}\,,\\
\mat{P^{\Lambda}_G} = \identity - \mat{P^{\Sigma}_G} = \mat{\Lambda} \left( \mat{\Lambda}^\mathrm{T} \mat{G_{\epsilon}} \mat{\Lambda} \right)^{-1} \mat{\Lambda}^\mathrm{T} \mat{G_{\epsilon}} \label{eq:scaled-projectors-loop-G}.
\end{gather}

From identities (<ref>) to (<ref>), we observe that these new projectors enable the proper re-scaling of the dominant term of the solenoidal part of the equation, that is the Gram matrix ($ \| \mat{P^{\Lambda}_{\identity}} \mat{G_{\epsilon}}\| \gg \|\mat{P^{\Lambda}_{\identity}} \mat{Z_{\Phi}}\|$). By leveraging them, we can build preconditioners to re-scale the D-VIE with the coefficients $\mat{C}_1$, $\mat{C}_2$, $\mat{C}_3$, and $\mat{C}_4$ similarly as in (<ref>) and (<ref>), except that these coefficients are now matrices. As in (<ref>), it is sufficient to employ only $\mat{C}_1$ and $\mat{C}_2$ to re-scale the D-VIE and set $\mat{C}_3$ and $\mat{C}_4$ to $\identity$. The resulting preconditioner is defined as follows

\begin{equation} \label{eq:precond_hc}
\mat{L_G} = \mat{C}_1 \mat{P^{\Lambda}_G} + \mat{C}_2 \mat{P^{\Sigma}_G},
\end{equation}

in which $\mat{C}_1 = {1}/{\| \mat{P^{\Lambda}_G} \|} \mat{G_{\epsilon}}^{-1} $ was determined using identities (<ref>) to (<ref>) to obtain an identity in the solenoidal part of the regularized D-VIE and $\mat{C}_2 = \gamma_\Sigma \identity = {1}/{\| \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} \mat{P^{\Sigma}_G} \|} \identity$ was chosen in the same way as in (<ref>). The resulting preconditioned matrix $\mat{M} = \mat{L_G}^\mathrm{T}\mat{Z}$ can be decomposed as in (<ref>)

\begin{equation}
\begin{split}
\mat{M} &= \mat{L_G}^\mathrm{T}\mat{Z}\\
&= \left( \mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} + \mat{C}_2 \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \right)\mat{Z}
\left( \mat{P^{\Lambda}_G} + \mat{P^{\Sigma }_G} \right)\\
&= \mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi})\mat{P^{\Lambda}_G} + \mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Sigma}_G} \\
&+{\gamma_\Sigma} \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Lambda}_G} + {\gamma_\Sigma} \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Sigma}_G} .
\label{eq:discretized_eq_decomp_precond}
\end{split}
\end{equation}

To derive the scalings of the four terms obtained in (<ref>), we first expand each of the terms as follows

\begin{gather}
\mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi})\mat{P^{\Lambda}_G} = \left( \mat{P^{\Lambda}_G} + \mat{G_{\epsilon}}^{-1} \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} \mat{P^{\Lambda}_G} \right)/\| \mat{P^{\Lambda}_G} \| , \label{eq:expansion_4blocks_precond_LL}\,,\\
\mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Sigma}_G} = \left(\mat{0} + \mat{G_{\epsilon}}^{-1} \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} \mat{P^{\Sigma}_G}\right)/\| \mat{P^{\Lambda}_G} \| , \label{eq:expansion_4blocks_precond_LS}\,,\\
{\gamma_\Sigma} \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Lambda}_G} ={\gamma_\Sigma} \left(\mat{0} + \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} \mat{P^{\Lambda}_G} \right), \label{eq:expansion_4blocks_precond_SL}\,,\\
{\gamma_\Sigma} \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Sigma}_G} = {\gamma_\Sigma} \left(\mat{G_{\epsilon}} \mat{P^{\Sigma}_G}
+ \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} \mat{P^{\Sigma}_G}\right) .
\label{eq:expansion_4blocks_precond_SS}
\end{gather}

Leveraging (<ref>) to (<ref>) and recalling that $\mat{P^{\Lambda}_G} = \mathcal{O}(1) $, $ \mat{G_{\epsilon}} \mat{P^{\Sigma}_G} = \mathcal{O}(\epsilon_r(\vr)^{-1})$, $\mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} = \mathcal{O}(\epsilon_r(\vr)^{-1})$, $ \mat{Z_\Phi} \mat{P^{\Lambda}_G} = \mathcal{O}(\epsilon_r(\vr)^{-1})$, $\mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} \mat{Z_\Phi} = \mathcal{O}(1)$, and $ \mat{Z_\Phi} \mat{P^{\Sigma}_G} = \mathcal{O}(1)$ when $\epsilon_r \rightarrow \infty$, we obtain that the four terms in (<ref>), representing the solenoidal and non-solenoidal parts of $\mat{M}$ and the associated cross terms, scale as

\begin{gather}
\mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi})\mat{P^{\Lambda}_G} = \mathcal{O}(1)\quad \epsilon_r(\vr) \rightarrow \infty, \label{eq:p_block_LL_precond}\,,\\
\mat{C}_1 \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Sigma}_G} =\mathcal{O}(1)\quad \epsilon_r(\vr) \rightarrow \infty, \label{eq:p_block_LS_precond}\,,\\
{\gamma_\Sigma} \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Lambda}_G} =\mathcal{O}(\epsilon_r(\vr)^{-1})\quad \epsilon_r(\vr) \rightarrow \infty, \label{eq:p_block_SL_precond}\,,\\
{\gamma_\Sigma} \mat{P^{\Sigma \mathrm{\normalfont{T}}}_G} (\mat{G_{\epsilon}} + \mat{Z_\Phi}) \mat{P^{\Sigma}_G} =\mathcal{O}(1)\quad \epsilon_r(\vr) \rightarrow \infty .
\label{eq:p_block_SS_precond}
\end{gather}

From (<ref>) to (<ref>), it results that the preconditioned matrix $\mat{M}$ obtained is well-conditioned. The final formulation reads

\begin{equation} \label{eq:precond_equation_G}
\mat{L_G}^\mathrm{T} \mat{Z} \mat{\alpha} = \mat{L_G}^\mathrm{T} \mat{v}.
\end{equation}

In the quasi-static regime, the condition number of the preconditioned system is independent of the frequency and the permittivity. However, if we now move away from this regime, we also need to consider the vector potential in our analysis. The solenoidal part of the preconditioned vector potential $\mat{L_G}^\mathrm{T} \mat{Z_{A}}$ scales as $\epsilon_r k_0^2$ and perturbs the conditioning of the preconditioned formulation if it becomes dominant compared to the solenoidal part of the preconditioned Gram matrix $\mat{L_G}^\mathrm{T} \mat{G_\epsilon}$, which scales as $\mathcal{O}(1)$. Therefore, to obtain a well-conditioned system matrix, we need the following inequality to hold true

\begin{equation}
\norm {\mat{G_{\epsilon}}^{-1} \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{G_{\epsilon}}\mat{P^{\Lambda}_G}} > \epsilon_r k_0^2 \norm {\mat{G}^{-1} \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{Z_{A,1}} \mat{P^{\Lambda}_G}},
\end{equation}

where $\mat{Z_{A}} = k_0^2 \mat{Z_{A,1}}$. This means that $\epsilon_r$ should be bounded as

\begin{equation}
\epsilon_r(\vr) < \frac{c}{k_0^2},
\end{equation}

where $c =\norm {\mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} } / \norm {\mat{G}^{-1} \mat{P^{\Lambda \mathrm{\normalfont{T}}}_G} \mat{Z_{A,1}} \mat{P^{\Lambda}_G}}$ is a constant.

It should be noted that in dielectric objects with particular geometries, Mie resonances can occur before the bound provided above. For conductive materials in which the dielectric relative permittivity $\epsilon_r'$ is negligible (i.e. $(\sigma / \omega \epsilon_0 \epsilon_r') \gg 1$), the bound for $\sigma$ can be expressed as

\begin{equation}
\sigma(\vr) < \frac{c}{k_0 \eta_0} .
\end{equation}

These ranges for the permittivity and the conductivity cover a large number of biological tissues at low frequencies ($f < \SI{1}{\mega \hertz}$). Above this frequency, we can set the coefficient matrices $\mat{C}_1$ and $\mat{C}_2$ to $\identity$ in (<ref>), which is equivalent to use the standard D-VIE. As the permittivity of biological tissues decreases when the frequency increases <cit.>, the standard D-VIE has a reasonably low condition number at high frequencies.

§ LOW-FREQUENCY AND HIGH-CONTRAST ANALYSES OF THE D-VIE

In this section we will first present an analysis of the low-frequency problems of the D-VIE. This will then be followed by an analysis of the high-contrast problems of the D-VIE. Leveraging these analyses and the nature of the criticalities observed, a cure for both problems will then be presented in <Ref>.

§.§ Low-frequency Analysis of the D-VIE

We analyze here the low-frequency behavior of the D-VIE when modeling lossy dielectric objects surrounded by free space. At low frequencies VIEs are subject to the low-frequency breakdown due to the frequency dependence of the complex permittivity of the object <cit.>. In this subsection we will confirm this with an analysis coherent with our framework that will then be used in the next section to construct the new formulation we propose in this work.
Here we will also show that for the discretized D-VIE, the equation ill-scaling not only causes the ill-conditioning of the associated linear system matrix, but also the loss of significant digits in the solution coefficients.
Our strategy will be to expose the low-frequency behavior of the D-VIE through a suitably normalized Loop-Star decomposition.

We propose here a new generalization of the Loop-Star decomposition, that differently from the standard decomposition, is not coefficient orthogonal, but rather an oblique decomposition. In its normalized form, it is defined as follows

\begin{equation} \label{eq:normDecomposition}
\mat{B^{A^{-1}}_{\Lambda \Sigma}} = \left[ \mat{\Lambda} (\transpose{\mat{\Lambda}} \mat{\Lambda})^{-\frac{1}{2}} ,	\mat{A} \mat{\tilde{\Sigma}} (\transpose{\mat{\tilde{\Sigma}}} \mat{A}^{2}\mat{\tilde{\Sigma}})^{-\frac{1}{2}}\right]\,,
\end{equation}

You do not explain why this is a decomposition: the proof for invertibility of this matrix works for real matrices only
AM: the usage of U like this is a bit sloppy, and you talk about 1/2 and you have only -1/2

where $\mat{A}$ is an invertible real symmetric matrix and for any invertible matrix $\mat{U}$, $\mat{U}^{-\frac{1}{2}}$ is the inverse of one of the (potentially many) square root matrices of $\mat{U}$. The fact that $\mat{B^{A^{-1}}_{\Lambda \Sigma}}$ is invertible is proven in Appendix <ref>. The reader should note that, due to the presence of the matrix square roots, the above decomposition is quite inefficient to implement. In this work, however, the decomposition will only be used as a theoretical tool in the conditioning analysis and, subsequently, as the starting point for defining our new oblique projectors. In other words, the computation of (<ref>) will never be required for the methods proposed in this paper and thus the often inefficient-to-compute matrix square roots will have no computational impact.
Since the matrix (<ref>) will be used to study the low-frequency limit, it must have a well-defined static limit. This is obtained by selecting $\mat{A}^{-1} = \mat{\tilde{G}_{\epsilon}}$ with $ \mat{\tilde{G}_\epsilon}= \mat{G_\epsilon}/\left(\mi \omega \right)$. It should also be noted that $\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}$ block-diagonalizes the D-VIE Gram matrix $\mat{G_\epsilon}$. In fact, since $(\transpose{\mat{\Lambda}} \mat{\Lambda})^{-\frac{1}{2}} \transpose{\mat{\Lambda}}
j\omega\mat{\tilde{G}_{\epsilon}}\mat{\tilde{G}_\epsilon}^{-1} \mat{\tilde{\Sigma}} (\transpose{\mat{\tilde{\Sigma}}} \mat{\tilde{G}_\epsilon}^{-2}\mat{\tilde{\Sigma}})^{-\frac{1}{2}}=\matO$ (because $\transpose{\mat{\Lambda}}\mat{\tilde{\Sigma}}=\matO$), we have

\begin{equation} \label{eq:identity_scaled_ls}
\transpose{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{{G}_{\epsilon}} \mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}
\transpose{\matN{\Lambda}}\mat{{G}_{\epsilon}}\matN{\Lambda}&\matO \\
\matO &\transpose{\matN{\Sigma}} \mat{{G}_{\epsilon}}\matN{\Sigma}
\end{bmatrix}\,,
\end{equation}

where $\matN{\Lambda} = \mat{\Lambda} (\transpose{\mat{\Lambda}} \mat{\Lambda})^{-\frac{1}{2}}$ and
$\matN{\Sigma} = \mat{\tilde{G}_\epsilon}^{-1} \mat{\tilde{\Sigma}} (\transpose{\mat{\tilde{\Sigma}}} \mat{\tilde{G}_\epsilon}^{-2}\mat{\tilde{\Sigma}})^{-\frac{1}{2}}$. Using $\mat{B^{\tilde{G}_{{\epsilon}}}_{\Lambda \Sigma}}$, the Loop-Star decomposition of the D-VIE can be expressed as

\begin{equation} \label{eq:decompositio}
\mat{Z_{\Lambda \Sigma}} \mat{\alpha_{\Lambda \Sigma}} =\transpose{\mat{B^{\tilde{G}_{{\epsilon}}}_{\Lambda \Sigma}}}(\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi}})\mat{B^{\tilde{G}_{{\epsilon}}}_{\Lambda \Sigma}} \mat{\alpha_{\Lambda \Sigma}}\,,
\end{equation}

where $\mat{B^{\tilde{G}_{{\epsilon}}}_{\Lambda \Sigma}} \mat{\alpha_{\Lambda \Sigma}} = \mat{\alpha}$. We further decompose $\mat{Z_{\Phi}}$ as
\mat{Z_{\Phi}} = \mat{Z_{\Phi, 11}} + \mat{Z_{\Phi,1 \epsilon}}

in which $\mat{Z_{\Phi, 11}}$ only accounts for the volumic contributions in the object and the surfacic contributions on its boundary ($\partial \Omega$)

\begin{equation} \label{eq:Z11-ns}
\begin{split}
&[\mat{Z_{\Phi,11}}]_{mn} =  \epsilon_0^{-1} . \\
&\left[ \int_\Omega \nabla \cdot \vt{f_m}(\vt{r}) \int_\Omega G_0(\vr,\vr') \kappa(\vr') \nabla \cdot \vt{f_n}(\vt{r'}) \diff v' \diff v \right. \\
&\left. - \int_\Omega \nabla \cdot \vt{f_m}(\vt{r}) \int_{\partial \Omega} G_0(\vr,\vr') \kappa^+_n \uv{n}_n \cdot \vt{f_n}(\vt{r'}) \diff s' \diff v \right.\\
&\left. - \int_{\partial \Omega} \uv{n}_m \cdot \vt{f_m}(\vt{r}) \int_{\Omega} G_0(\vr,\vr') \kappa(\vr')\nabla \cdot \vt{f_n}(\vt{r'}) \diff v' \diff s \right.\\
&\left.+ \int_{\partial \Omega} \uv{n}_m \cdot \vt{f_m}(\vt{r}) \int_{\partial \Omega} G_0(\vr,\vr') \kappa^+_n \uv{n}_n \cdot \vt{f_n}(\vt{r'}) \diff s' \diff s \right]\,,
\end{split}
\end{equation}

with $\kappa^+_n$ being the dielectric contrast in tetrahedron $T^+_n$ and $\uv{n}_n$ the outward unit normal of the triangle of $T^+_n$ which pertains to $\partial \Omega$; instead $\mat{Z_{\Phi, 1 \epsilon}}$ includes the purely surfacic contributions internal to the object (at the interface $\Gamma_n$ between two tetrahedra $T^-_n$ and $T^+_n$ of different material contrasts $\kappa^-_n$ and $\kappa^+_n$)

\begin{equation} \label{eq:Z1eps-ns}
\begin{split}
&[\mat{Z_{\Phi,1\epsilon}}]_{mn} = \epsilon_0^{-1} . \\
&\left[\int_\Omega \nabla \cdot \vt{f_m}(\vt{r}) \int_{\Gamma_n} G_0(\vr,\vr') \deltaKappa_n \uv{n}_n \cdot \vt{f_n}(\vt{r'}) \diff s' \diff v \right.\\
&\left.- \int_{\partial \Omega} \uv{n}_m \cdot \vt{f_m}(\vt{r}) \int_{\Gamma_n} G_0(\vr,\vr') \deltaKappa_n \uv{n}_n \cdot \vt{f_n}(\vt{r'}) \diff s' \diff s \right]\,,
\end{split}
\end{equation}

in which $\deltaKappa_n = \kappa^-_n - \kappa^+_n$ and $\uv{n}_n$ is a unit vector normal to $\Gamma_n$ oriented from $T^+_n$ to $T^-_n$.

When applying $\matN{\Lambda}^\mathrm{T}$ and $\matN{\Lambda}$ to the scalar potential matrices, we obtain the relations
$ \matN{\Lambda}^\mathrm{T} \mat{Z_{\Phi, 11}}= \matO$, $\matN{\Lambda}^\mathrm{T} \mat{Z_{\Phi, 1\epsilon}}= \matO$, and $\mat{Z_{\Phi, 11}} \matN{\Lambda} = \matO$.

The fact that $\mat{\Lambda}$ (and thus $\matN{\Lambda}$) cancels the surface terms defined on $\partial \Omega$ in $\mat{Z_{\Phi, 11}}$ and $\mat{Z_{\Phi, 1\epsilon}}$ results from the absence of Loop functions on $\partial \Omega$. This property is proven in Appendix <ref>.

We then represent $\mat{Z_{\Lambda \Sigma}}$ in (<ref>) as a 2-by-2 block matrix

\begin{gather}
\mat{Z_{\Lambda \Sigma}}= \transpose{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{Z} \mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}
= \begin{bmatrix}
\matN{\Lambda}^\mathrm{T} \mat{Z} \matN{\Lambda} & \matN{\Lambda}^\mathrm{T} \mat{Z} \matN{\Sigma}\\
\matN{\Sigma}^\mathrm{T}\mat{Z} \matN{\Lambda} & \matN{\Sigma}^\mathrm{T} \mat{Z} \matN{\Sigma}
\end{bmatrix}\,, \text{ in which}\label{eq:block_mat_loop-star-ns} \\
\matN{\Lambda}^\mathrm{T} \mat{Z} \matN{\Lambda} = \matN{\Lambda}^\mathrm{T} (\mat{G_{\epsilon}} + \mat{Z_A})\matN{\Lambda} \label{eq:decomp_LL_ns}\,,\\
\matN{\Lambda}^\mathrm{T} \mat{Z}\matN{\Sigma} = \matN{\Lambda}^\mathrm{T} \mat{Z_A} \matN{\Sigma} \label{eq:decomp_LS_ns}\,,\\
\matN{\Sigma}^\mathrm{T} \mat{Z} \matN{\Lambda} = \matN{\Sigma}^\mathrm{T} (\mat{Z_A} + \mat{Z_{\Phi,1 \epsilon}}) \matN{\Lambda} \label{eq:decomp_SL_ns}\,,\\
\matN{\Sigma}^\mathrm{T} \mat{Z} \matN{\Sigma} = \matN{\Sigma}^\mathrm{T} (\mat{G_{\epsilon}} + \mat{Z_A} + \mat{Z_{\Phi}}) \matN{\Sigma}\,.
\label{eq:decomp_SS_ns}
\end{gather}

To identity the frequency behavior of the terms in (<ref>), (<ref>), (<ref>), and (<ref>), we first need to introduce the low-frequency behavior of the real and imaginary parts of the material parameters $\epsilon$, $\kappa$, and $\deltaKappa$
AM: introduce the notation $\Re$ and $\Im$?

\begin{gather}
\Re(\epsilon) \underrel{\omega \to 0}{=} \mathcal{O}(1), ~
\Im(\epsilon) \underrel{\omega \to 0}{=} \mathcal{O}(1/\omega) \label{eq:eps_scaling}\,,\\
\Re(1/\epsilon) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(1/\epsilon) \underrel{\omega \to 0}{=} \mathcal{O}(\omega) \label{eq:eps_inv_scaling}\,,\\
\Re(\kappa) \underrel{\omega \to 0}{=} \mathcal{O}(1), ~
\Im(\kappa) \underrel{\omega \to 0}{=} \mathcal{O}(\omega) \label{eq:kappa_scaling}\,,\\
\Re(\deltaKappa) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(\deltaKappa)\underrel{\omega \to 0}{=} \mathcal{O}(\omega)\,, \label{eq:delta_kappa_scaling}
\end{gather}

which derives from the definitions of the complex permittivity and the dielectric contrast. Then, using the definitions of $\mat{G_\epsilon}$, $\mat{Z_A}$, $\mat{Z_{\Phi,11}}$, and $\mat{Z_{\Phi,1\epsilon}}$ in (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>), we can deduce the frequency scalings of these matrices from (<ref>), (<ref>), (<ref>), and (<ref>). We finally obtain the following low-frequency behavior for the real and imaginary parts of $\mat{G_\epsilon}$, $\mat{Z_A}$, $\mat{Z_{\Phi,11}}$, and $\mat{Z_{\Phi,1\epsilon}}$
these are not the decomposition terms since you do not have $\Lambda$ and $\Sigma$

\begin{gather}
\Re(\mat{G_{\epsilon}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(\mat{G_{\epsilon}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega) \label{eq:scaling-gram-ns}\,,\\
\Re(\mat{Z_A}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(\mat{Z_A}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^3) \label{eq:scaling-A-ns}\,,\\\
\Re(\mat{Z_{\Phi,11}}) \underrel{\omega \to 0}{=} \mathcal{O}(1), ~
\Im(\mat{Z_{\Phi,11}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega) \label{eq:scaling-11-ns}\,,\\
\Re(\mat{Z_{\Phi,1 \epsilon}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(\mat{Z_{\Phi,1 \epsilon}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega)\,. \label{eq:scaling-1eps-ns}
\end{gather}

Using (<ref>) to (<ref>), we deduce the scalings for $\mat{Z_{\Lambda \Sigma}}$ at low frequencies

\begin{gather}
\Re\left(\mat{Z_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) &\mathcal{O}(\omega^2) \\
\mathcal{O}(\omega^2) & \mathcal{O}(1)
\end{bmatrix} \label{eq:scaling_mat_real_ns}\,,\\
\Im\left( \mat{Z_{\Lambda \Sigma}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) &\mathcal{O}(\omega^3) \\
\mathcal{O}(\omega) & \mathcal{O}(\omega)
\end{bmatrix} \label{eq:scaling_mat_imag_ns}\,,
\end{gather}

which, according to the Gershgorin circle theorem, confirm that $\mat{Z_{\Lambda \Sigma}}$ is ill-conditioned for $\omega\to 0$. Since
$\mat{B^{\tilde{G}_{{\epsilon}}}_{\Lambda \Sigma}}$ is a well-conditioned matrix (refer to Appendix <ref>), the ill-conditioning of $\mat{Z_{\Lambda \Sigma}}$ implies the ill-conditioning of $\mat{Z}$ for $\omega\to 0$.

Besides the conditioning of the system matrix, it is also important to determine whether or not the solution coefficients and the right hand side vectors are preserved in the static limit. We first provide the scalings of the right hand side for plane wave and dipole excitations, which are frequently employed in bioelectromagnetic applications. For a plane wave excitation, the Loop and Star components of the right hand side $\mat{v_{\textnormal{{PW}}}}$ have the following scalings when $\omega\rightarrow 0$

\begin{equation} \label{eq:scaling_rhs_pw_ns}
\Re \left(\transpose{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v_{\textnormal{{PW}}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\transpose{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v_{\textnormal{{PW}}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix},
\end{equation}

AM: In 3d, the integration of a loop is still 0, correct?
in which $\left[ \mat{v_{\textnormal{PW}}} \right]_{m} = \int_\Omega \vt{f_m}(\vr) \cdot \vt{E}_0 \exp{(-\mi \vt{k} \cdot \vr )} \diff v$ with $\vt{E}_0$ being the polarization of the plane wave and $\vt{k}$ its wave vector. Note that the right hand side in the D-VIE for a plane wave excitation scales similarly to the plane wave right hand side of the surface electric field integral equation <cit.>. To derive the frequency dependence of the Loop-Star decomposition of a dipole excitation, we start from the expression of the field it radiates (<ref>) in which $\vt{J}_i(\vr, \vr_0) = \mi \omega \delta(\vr - \vr_0)\vt{p}$ is the current dipole with $\delta$, $\vt{p}$, and $\vt{r_0}$ being the Dirac delta function, the dipole moment, and the dipole position, respectively. We then test (<ref>) with $\{\vt{f_m}\}$ to obtain the discretized right hand side

\begin{equation}
\begin{split}
\left[ \mat{v_{\textnormal{dip}}} \right]_{m} &= \left[ \mat{v_{\mathrm{A}}} \right]_{m} + \left[ \mat{v_{\mathrm{\Phi}}} \right]_{m}\\
&= \frac{k_0^2}{\epsilon_0} \vt{p} \cdot \int_\Omega \vt{f_m}(\vr) G_0(\vr, \vr_0) \diff v\\
&-\frac{1}{ \epsilon_0}\int_{\Omega} \nabla \cdot \vt{f_m}(\vr) \vt{p} \cdot \nabla G_0(\vr, \vr_0) \diff v\,.
\end{split}
\end{equation}

From the Taylor expansion of $G_0$ and its gradient, we obtain that $\Re(G_0(\vr,\vr_0)) = \mathcal{O}(1)$, $\Im(G_0(\vr,\vr_0)) = \mathcal{O}(\omega)$, $\Re(\nabla G_0(\vr,\vr_0)) = \mathcal{O}(1)$, and $\Im(\nabla G_0(\vr,\vr_0)) = \mathcal{O}(\omega^3)$ at low frequencies. This yields the following frequency dependencies for the real and imaginary parts of the Loop-Star decomposition of $\mat{v_{\mathrm{A}}}$ and $\mat{v_{\mathrm{\Phi}}}$

\begin{gather}
\Re(\matN{\Lambda}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(\matN{\Lambda}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^5)\,,\\
\Re(\matN{\Sigma}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^2), ~
\Im(\matN{\Sigma}^\mathrm{T} \mat{v_{\mathrm{A}}}) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^3)\,,\\
\Re(\matN{\Lambda}^\mathrm{T} \mat{v_{\mathrm{\Phi}}}) = 0, ~
\Im(\matN{\Lambda}^\mathrm{T} \mat{v_{\mathrm{\Phi}}}) = 0\,,\\
\Re(\matN{\Sigma}^\mathrm{T} \mat{v_{\mathrm{\Phi}}} ) \underrel{\omega \to 0}{=} \mathcal{O}(1), ~
\Im(\matN{\Sigma}^\mathrm{T} \mat{v_{\mathrm{\Phi}}} ) \underrel{\omega \to 0}{=} \mathcal{O}(\omega^3)\,,
\end{gather}

AM: $\matN{\Lambda}^\mathrm{T} \mat{v_{\mathrm{\Phi}}}$ is not an exact 0? Why the limit?

from which we obtain the scalings of $\mat{v_{\mathrm{dip}}}$ at low frequencies

\begin{equation} \label{eq:scaling_rhs_dipole_ns}
\Re \left(\transpose{\mat{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}}} \mat{v_{\textnormal{dip}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\transpose{\mat{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}}} \mat{v_{\textnormal{dip}}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^5) \\
\mathcal{O}(\omega^3)
\end{bmatrix}.
\end{equation}

This concludes the scaling analysis of the right hand side for plane wave and dipole excitations. The scalings obtained are summarized in <Ref>, in which $\mat{v_\Lambda}$ denotes the Loop part of the right hand side vector and $\mat{v_\Sigma}$ its Star part. In the Loop-Star D-VIE, all the terms in (<ref>) are preserved since they are stored separately. However, in the standard D-VIE, the terms recovered correspond only to the dominant terms (real and imaginary) in (<ref>), all the other terms are lost due to finite precision arithmetic.

AM: this does not sound correct, if the scalings are not the same, the dominant one remains
AM: You should reming the reader that this will only happen if you do not use loop-star

To identify the impact of this loss on the solution coefficients, we subsequently retrieve the scalings of $\mat{\alpha_{\Lambda \Sigma}}$ for plane wave and dipole excitations. This requires the knowledge of the frequency scalings of the inverse of $\mat{Z_{\Lambda \Sigma}}$, which are derived by inverting (<ref>) using Schur complement formulas <cit.>

\begin{gather}
\Re\left(\mat{Z_{\Lambda \Sigma}}^{-1} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) &\mathcal{O}(\omega^2) \\
\mathcal{O}(1) & \mathcal{O}(1)
\end{bmatrix} \label{eq:scaling_mat_inv_real_ns}\,,\\
\Im\left( \mat{Z_{\Lambda \Sigma}}^{-1} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1/\omega) &\mathcal{O}(\omega) \\
\mathcal{O}(\omega) & \mathcal{O}(\omega)
\end{bmatrix} \label{eq:scaling_mat_inv_imag_ns}.
\end{gather}

Finally by multiplying the scaling matrix of $\mat{Z_{\Lambda \Sigma}}^{-1}$ ((<ref>) and (<ref>)) and the scaling vector of $\transpose{\mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v}$ ((<ref>) for the plane wave or (<ref>) for the dipole), we obtain the following scalings of $\mat{\alpha_{\Lambda \Sigma}}$

\begin{gather}
\Re \left(\mat{\alpha_{\Lambda \Sigma}^\text{PW}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\mat{\alpha_{\Lambda \Sigma}^\text{PW}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix}\,,  \label{eq:scaling_solution_pw}\\
\Re \left(\mat{\alpha_{\Lambda \Sigma}^\text{dip}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^2) \\
\mathcal{O}(1)
\end{bmatrix};
\Im \left(\mat{\alpha_{\Lambda \Sigma}^\text{dip}} \right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix}\label{eq:scaling_solution_dipole}\,,
\end{gather}

for plane wave and dipole excitations, respectively.

Frequency scalings of the real and imaginary parts of the right hand side vectors $\mat{v}$ and $\mat{\tilde{v}}$ for plane wave and dipole excitations.

[Right hand side $\mat{v}$]

2*($\Re$ , $\Im$) $(\mat{v_\Lambda} )$

2*$(\Re , \Im)( \mat{v_\Sigma})$

Terms     Terms required

recovered     for a correct solution

2*Plane Wave

2*$(\omega^2 , \omega)$

2*$(1 ,\omega)$

$\Im(\mat{v_\Lambda}),\Im(\mat{v_\Sigma})$    $\Im(\mat{v_\Lambda}),\Im(\mat{v_\Sigma})$

$(\omega^2 , \omega^5)$

$(1 , \omega^3)$

[Solution $\mat{\alpha_{\Lambda \Sigma}}$]

Source     ($\Re$ , $\Im$) $(\mat{\alpha_{\Lambda }})$    $(\Re , \Im)(\mat{\alpha_{\Sigma}})$     Terms recovered

2*Plane Wave     2*$(1 ,\omega)$     2*$(1 ,\omega)$      $\Re(\mat{\alpha_{\Lambda }}),\Re(\mat{\alpha_{\Sigma}})$
$\Im(\mat{\alpha_{\Lambda }}),\Im(\mat{\alpha_{ \Sigma}})$

Dipole     2*$(\omega^2 ,\omega)$      2*$(1 ,\omega^3)$     $\Re(\mat{\alpha_{ \Sigma}})$

[Scaled right hand side $\mat{\tilde{v}}$ ]

2*($\Re$ , $\Im$) $(\mat{\tilde{v}_\Lambda} )$

2*$(\Re , \Im)( \mat{\tilde{v}_\Sigma})$

Terms     Terms required

recovered     for a correct solution

2*Plane Wave

2*$(1 ,\omega)$

2*$(1 ,\omega)$

$\Re(\mat{\tilde{v}_\Lambda}),\Re(\mat{\tilde{v}_\Sigma})$     $\Re(\mat{\tilde{v}_\Lambda}),\Re(\mat{\tilde{v}_\Sigma})$

$\Im(\mat{\tilde{v}_\Lambda}),\Im(\mat{\tilde{v}_\Sigma})$     $\Im(\mat{\tilde{v}_\Lambda}),\Im(\mat{\tilde{v}_\Sigma})$

$(\omega^4 ,\omega)$

$(1 ,\omega^3)$

At low frequencies, the real (or imaginary) part of the Loop and Star components of the solution coefficients in (<ref>) have the same frequency scalings, hence all the components of the solution for a plane wave excitation are preserved in this regime. However, the real part of the Loop component of the solution for a dipole excitation in (<ref>) becomes much smaller than the real part of its Star component at low frequencies and is lost due to finite precision arithmetic. In the scope of this paper, the D-VIE is applied in scenarios where the electric field inside the object is required. Since there is a simple scalar relation between the unknown of the D-VIE (i.e. the electric flux) and the electric field, the re-amplification of a lost (Loop or Star) part is not possible, and thus it is not required to preserve the components of the solution lost due to finite precision at low frequencies.

Nevertheless, a possible loss of accuracy in the solution occurs when a term of the right hand side, which is lost at low frequencies, contributes to one of the dominant terms of the solution.
AM: we might want to make this discussion more vague, since it might not be the only effect into play: Re-ask In <ref>, the terms contributing to the solution for both types of excitation are provided. For the dipole and the plane wave, the real part of the Loop component of the right hand side is lost due to finite precision while it is supposed to contribute to the dominant terms of the solution. For this reason, there is a loss of accuracy in the solution at low frequencies with plane wave and dipole excitations.

§.§ High-Contrast Analysis of the D-VIE

In the following we analyze the high-contrast behavior of the D-VIE for complex permittivity objects with dominant imaginary part ($\sigma / (\omega \epsilon_0 \epsilon_r') \gg 1$) at low frequencies. In this regime, the dielectric contrast $\kappa$ between the object and the background is approximately 1 and the Green's function $G_0$ can be bounded from above by a frequency-independent term.

In our theoretical treatment, the scatterer $\Omega$ is supposed to be a piecewise homogeneous object composed of a region $\Omega_M$ which has the maximum conductivity $\sigma_{\text{max}}$ of $\Omega$ (associated to the complex permittivity $\epsilon_{\text{max}}$) while the conductivity is bounded in its complementary domain $\Omega \backslash \Omega_M$. In the following, the boundary between $\Omega_M$ and the rest of the object is denoted by $\partial \Omega_M$ and the maximum conductivity ratio is defined as $r_\sigma = \sigma_\text{max}/\sigma_\text{min}$ with $\sigma_\text{min}$ being the minimum conductivity in $\Omega$. The following high-contrast analysis will be performed for $r_\sigma$ going toward infinity with $\sigma_\text{min}$ fixed. Note that this scenario represents the high internal contrast occurring in lossy dielectric objects in several application scenarios.

AM: it sounds a bit fishy to claim that the ratio going to infinity is relevant for biomedical application.
Also, in this section the frequency is fixed, correct? If so it mght need clarifying: SHould be fixed now (removed source loc application)

Similarly to the low-frequency breakdown, the HC problem originates from the solenoidal part of $\mat{Z}$ that is ill-scaled due to the permittivity scaling of the Gram matrix $\mat{G_\epsilon}$. While the approach carried out in the low-frequency analysis only required the determination of the frequency scaling of $\|\mat{G_\epsilon}\|$, some additional considerations on the material dependence of the minimum singular value of the Gram matrix are needed in the HC analysis. To this aim, we leverage the block structure of $\mat{G_\epsilon}$ for the above-mentioned scatterer. The entries of $\mat{G_\epsilon}$ can be written as

\left[ \mat{G_{\epsilon}} \right]_{mn} = \int_{T_m^+} \vt{f_m} \cdot \epsilon^{-1} \vt{f_n} \diff v +
\int_{T_m^-} \vt{f_m} \cdot \epsilon^{-1}  \vt{f_n} \diff v\,,

with $T^\pm_m$ being the tetrahedra on which $\vt{f_m}$ is defined. In the case of a half basis function, only $T^+_m$ is required.

From this definition, we can build the following block matrix

\begin{gather}
\mat{G_\epsilon}
= \begin{bmatrix}
\mat{G_M} & \transpose{\mat{B_{M}}} \\
\mat{B_{M}} & \mat{G_R}
\end{bmatrix}\,, \text{ in which} \label{eq:block-structure} \\
\mat{G_{M}}  = \mat{G_\epsilon}[a_1,\dots, a_{N_{FM}};a_1,\dots, a_{N_{FM}}]\,, \\
\mat{B_{M}}  = \mat{G_\epsilon}[a_{N_{FM} + 1},\dots, a_{N_F};a_1,\dots, a_{N_{FM}}] \,,\\
\mat{G_{R}}  = \mat{G_\epsilon}[a_{N_{FM} + 1},\dots, a_{N_F};a_{N_{FM} + 1},\dots, a_{N_F}]\,,
\end{gather}

with $\{a_1,\dots, a_{N_{FM}}\}$ being the indices of the rows (and columns) of $ \mat{G_\epsilon}$ corresponding to the $N_{FM}$ SWG basis functions that have both of their supporting tetrahedra in $\Omega_M$ and $\{a_{N_{FM} + 1},\dots, a_{N_F}\}$ being the ${N_F-N_{FM}}$ indices of the remaining SWG basis functions.

AM: the matrix indices are not consistent
AM: much better, ask francesco if the matlab notation should be introduced though (CH: source: https://proofwiki.org/wiki/Definition:Submatrix/Notation)

Using the block structure given in (<ref>) and the fact $\Re(1/\epsilon_{\text{max}} ) = \mathcal{O}(r_{\sigma}^{-2})$, $\Im(1/\epsilon_{\text{max}} ) = \mathcal{O}(r_{\sigma}^{-1})$, $\mat{G_R} =\mathcal{O}(1/\epsilon)$ (bounded complex permittivity), $\mat{B_M} =\mathcal{O}(1/\epsilon_{\text{max}} )$, and $\mat{G_M} =\mathcal{O}(1/\epsilon_{\text{max}} )$ when $r_\sigma \rightarrow \infty$, we obtain the following scalings for the blocks of $\mat{G_\epsilon}$ when the contrast goes to infinity

\begin{gather}
\Re \left(\mat{G_\epsilon}\right) \underrel{r_\sigma \to \infty}{=} \begin{bmatrix}
\mathcal{O}(1 / r_\sigma^2) & \mathcal{O}(1 / r_\sigma^2) \\
\mathcal{O}(1 / r_\sigma^2) & \mathcal{O}(1)
\end{bmatrix} \label{eq:scaling_Gram_Re}\,,\\
\Im \left(\mat{G_\epsilon}\right) \underrel{r_\sigma \to \infty}{=} \begin{bmatrix}
\mathcal{O}(1 / r_\sigma) & \mathcal{O}(1 / r_\sigma) \\
\mathcal{O}(1 / r_\sigma) & \mathcal{O}(1)
\end{bmatrix} ,\label{eq:scaling_Gram_Im}
\end{gather}

which, according to the Gershgorin circle theorem, show that $\mat{G_\epsilon}$ suffers from ill-conditioning in objects with high internal contrast. This HC problem in the Gram matrix can be a source of ill-conditioning in the discretized D-VIE, and hence should be considered when regularizing the D-VIE for high-contrast.

AM: this statement seems a bit abrupt.

§ OBLIQUE QUASI-HELMHOLTZ PROJECTORS

Just like the orthogonal quasi-Helmholtz projectors <cit.> arise from the inversion of the standard Loop-Star decomposition, the oblique Loop-Star decomposition introduced in equation (<ref>) can give rise to an entire new family of projectors that, in general, will not be orthogonal, but oblique.
So starting from the decomposition described in (<ref>) (omitting the normalization which is not required for the following and using $\mat{\Sigma}$ instead of $\mat{\tilde{\Sigma}}$), we write the solution coefficient vector as

\begin{equation}\label{eq:obliquedec}
\mat{\alpha} = \mat{\Lambda} \mat{\tilde{l}} + \mat{A} \mat{\Sigma} \mat{\tilde{s}}\,,
\end{equation}

pb coherence when building projectors: $\mat{\Sigma}$ has its columns linearly independent

in which as before $\mat{A}$ is an invertible symmetric matrix.
Note that (<ref>) is a valid Helmholtz decomposition, as demonstrated in Appendix <ref>.
We can now obtain a new set of projectors by solving for $\mat{\Lambda} \mat{\tilde{l}}$ and for $\mat{A} \mat{\Sigma} \mat{\tilde{s}}$ separately. To this end, we first multiply (<ref>) by $\mat{\Lambda}^\mathrm{T} \mat{A}^{-1}$ and $\mat{\Sigma}^\mathrm{T} $, and we obtain the following equations

\begin{equation}
\mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\alpha} = \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\Lambda} \mat{\tilde{l}};~
\mat{\Sigma}^\mathrm{T} \mat{\alpha} = \mat{\Sigma}^\mathrm{T} \mat{A} \mat{\Sigma} \mat{\tilde{s}}\label{eq:app1_decomp_lp}\,,
\end{equation}

in which we used the properties that $\mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{A} \mat{\Sigma} = \mat{0}$ and $\mat{\Sigma}^\mathrm{T} \mat{\Lambda} = \mat{0}$. From (<ref>), we can now express the coefficients of the oblique Loop and Star basis functions from the coefficients of the original basis functions

\begin{equation}
\mat{\tilde{l}} = (\mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\Lambda})^{-1} \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\alpha};~
\mat{\tilde{s}} = (\mat{\Sigma}^\mathrm{T} \mat{A} \mat{\Sigma})^+ \mat{\Sigma}^\mathrm{T} \mat{\alpha} \label{eq:coeff-loopstar}\,,
\end{equation}

where $^+$ is the Moore-Penrose pseudo inverse, which is required in (<ref>) due to the one-dimensional null space of $\mat{\Sigma}^\mathrm{T} \mat{A} \mat{\Sigma}$.

AM: pinv notation not introduced
The next step is to obtain the solenoidal and non-solenoidal parts of $\mat{\alpha}$ in terms of $\mat{\alpha}$ itself by applying $\mat{\Lambda} $ and $\mat{A} \mat{\Sigma} $ to the equations in (<ref>)

\begin{gather}
\mat{\Lambda} \mat{\tilde{l}} = \mat{\Lambda} (\mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\Lambda})^{-1} \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\alpha} \label{eq:dec1}\,,\\
\mat{A} \mat{\Sigma} \mat{\tilde{s}} = \mat{A} \mat{\Sigma} (\mat{\Sigma}^\mathrm{T} \mat{A} \mat{\Sigma})^+ \mat{\Sigma}^\mathrm{T} \mat{\alpha}\,.\label{eq:dec2}
\end{gather}

What we have obtained is a family of two complementary oblique projectors

\begin{gather}
\mat{P^{\Lambda}_{A^{-1}}} = \mat{\Lambda} \left( \mat{\Lambda}^\mathrm{{\mathrm{T}}} \mat{A}^{-1} \mat{\Lambda} \right)^{-1} \mat{\Lambda}^\mathrm{{\mathrm{T}}} \mat{A}^{-1}\,,\\
\mat{P^{\Sigma}_{A^{-1}}} = \mat{A} \mat{\Sigma} \left( \mat{\Sigma}^\mathrm{{\mathrm{T}}} \mat{A} \mat{\Sigma} \right)^{+} \mat{\Sigma}^\mathrm{{\mathrm{T}}}\,,
\end{gather}

which are the oblique quasi-Helmholtz projectors we propose in this work. From (<ref>), (<ref>), and (<ref>) it follows that

\mat{\alpha} = (\mat{P^{\Lambda}_{A^{-1}}} + \mat{P^{\Sigma}_{A^{-1}}}) \mat{\alpha}

and this proves the complementarity property

\mat{P^{\Lambda}_{A^{-1}}} + \mat{P^{\Sigma}_{A^{-1}}} = \identity

These new projectors are, in general, not symmetric and their transposes will be denoted by $\mat{P^{\Lambda \mathrm{\normalfont{T}}}_{A^{-1}}}$ and $\mat{P^{\Sigma \mathrm{\normalfont{T}}}_{A^{-1}}}$. The following properties can be easily proven

\begin{gather}
\mat{P^{\Lambda \mathrm{\normalfont{T}}}_{A^{-1}}} \mat{A^{-1}} \mat{P^{\Lambda}_{A^{-1}}} = \mat{A^{-1}} \mat{P^{\Lambda}_{A^{-1}}} ;~
\mat{P^{\Lambda \mathrm{\normalfont{T}}}_{A^{-1}}} \mat{A^{-1}} \mat{P^{\Sigma}_{A^{-1}}} = \matO \label{eq:id2_scaled_proj_ns_A}\,,\\
\mat{P^{\Sigma \mathrm{\normalfont{T}}}_{A^{-1}}} \mat{A^{-1}} \mat{P^{\Lambda}_{A^{-1}} } = \matO ;~
\mat{P^{\Sigma \mathrm{\normalfont{T}}}_{A^{-1}}} \mat{A^{-1}} \mat{P^{\Sigma}_{A^{-1}}} =  \mat{A^{-1}} \mat{P^{\Sigma}_{A^{-1}}}\,.
\label{eq:id4_scaled_proj_ns_A}
\end{gather}

When $\mat{A} = \identity$, we obtain the standard Loop-Star decomposition and the associated projectors are the standard quasi-Helmholtz projectors <cit.>. The above formulas define in general oblique projectors for all formulations where a valid Loop-Star decomposition can be defined, be it 2D, 3D surface, or 3D volume. For specializing to the volume equation of interest in this work, however, we set $\mat{A} = \mat{G_\epsilon}^{-1}$ and we obtain the scaled projectors $\mat{P^{\Sigma}_{G_\epsilon}}$ and $\mat{P^{\Lambda}_{G_\epsilon}}$ defined as

\begin{gather}
\mat{P^{\Sigma}_{G_\epsilon}} = \mat{G_{\epsilon}}^{-1} \mat{\Sigma} \left( \mat{\Sigma}^\mathrm{T} \mat{G_{\epsilon}}^{-1} \mat{\Sigma} \right)^{+} \mat{\Sigma}^\mathrm{T} \label{eq:scaled-projectors-star-G-ns}\,,\\
\mat{P^{\Lambda}_{G_\epsilon}} = \identity - \mat{P^{\Sigma}_{G_\epsilon}} = \mat{\Lambda} \left( \mat{\Lambda}^\mathrm{T} \mat{G_{\epsilon}} \mat{\Lambda} \right)^{-1} \mat{\Lambda}^\mathrm{T} \mat{G_{\epsilon}}\,, \label{eq:scaled-projectors-loop-G-ns}
\end{gather}

for which the associated cancellation properties become

\begin{gather}
\mat{P^{\Lambda \mathrm{\normalfont{T}}}_{G_\epsilon}} \mat{G_{\epsilon}} \mat{P^{\Lambda}_{G_\epsilon}} = \mat{G_{\epsilon}} \mat{P^{\Lambda}_{G_\epsilon}} ;~
\mat{P^{\Lambda \mathrm{\normalfont{T}}}_{G_\epsilon}} \mat{G_{\epsilon}} \mat{P^{\Sigma}_{G_\epsilon}} = \matO \,, \label{eq:id1_scaled_proj_ns}\\
\mat{P^{\Sigma \mathrm{\normalfont{T}}}_{G_\epsilon}} \mat{G_{\epsilon}} \mat{P^{\Lambda}_{G_\epsilon} } = \matO ;~
\mat{P^{\Sigma \mathrm{\normalfont{T}}}_{G_\epsilon}} \mat{G_{\epsilon}} \mat{P^{\Sigma}_{G_\epsilon}} =  \mat{G_{\epsilon}} \mat{P^{\Sigma}_{G_\epsilon}}\,. \label{eq:id4_scaled_proj_ns}
\end{gather}

Moreover, when applying $\mat{P^{\Lambda}_{G_\epsilon}}$ to the scalar potential matrices, we obtain the relations

\mat{P^{\Lambda}_{G_\epsilon}} \mat{G_\epsilon}^{-1} \mat{Z_{\Phi, 11}}= \matO$, $\mat{P^{\Lambda}_{G_\epsilon}} \mat{G_\epsilon}^{-1}  \mat{Z_{\Phi, 1\epsilon}}= \matO$, and $\mat{Z_{\Phi, 11}} \mat{P^{\Lambda}_{G_\epsilon}}  = \matO$.

AM: recall the cancellation properties on the scalar potentials ?

§ THE NEW REGULARIZED VOLUME INTEGRAL EQUATION

To find the correct parameters of a regularizer for the D-VIE based on the new oblique quasi-Helmholtz projectors just introduced, we will employ an approach often and successfully used with projectors at low-frequency: first a regularization is obtained by determining the parameters of a, impractical-to-implement but theoretically useful, normalized Loop-Star decomposition; then the same coefficients are used for linearly combining solenoidal and non-solenoidal projectors. In the present context, however, we not only have to regularize the low-frequency breakdown, but also the high-contrast breakdown. Thus, after finding the frequency-regularizing coefficients with a normalized Loop-Star analysis we will propose a formulation that also solves the high-contrast breakdown by using an ansatz whose effectiveness will then be proved in the following section.

Consider the following normalized and re-scaled Loop-Star decomposition of the D-VIE

\begin{equation}
\label{eq:helmholtz_decomp_precond_ns}
\mat{\tilde{Z}_{\Lambda \Sigma}}\mat{{\alpha}_{\Lambda \Sigma}}=
\begin{bmatrix}
C_1 \matN{\Lambda}^\mathrm{T} \\
C_2 \matN{\Sigma}^\mathrm{T}
\end{bmatrix} \mat{Z} [C_3 \matN{\Lambda},~C_4 \matN{\Sigma}] \mat{{\alpha}_{\Lambda \Sigma}}
= \begin{bmatrix}
C_1 \matN{\Lambda}^\mathrm{T} \\
C_2 \matN{\Sigma}^\mathrm{T}
\end{bmatrix} \mat{v}=\mat{\tilde{v}}\,,
\end{equation}

in which the coefficients $C_1$, $C_2$, $C_3$ and $C_4$ have to be determined to make the Loop and Star parts of the preconditioned matrix ($\mat{\tilde{Z}_{\Lambda \Sigma}}$) and its right hand side $\mat{\tilde{v}}$ free from ill-conditioning and loss of accuracy. Since the numerical loss comes from the right hand side and not from the solution directly, we can set $C_3$ and $C_4$ to $1$ and attempt to leverage the two remaining coefficients to regularize the rest of the equation. Therefore, we now need to determine $C_1$ and $C_2$ for $\mat{\tilde{Z}_{\Lambda \Sigma}}$ to be well-conditioned and to avoid numerical loss. It will be now shown that, by choosing

\begin{equation}\label{eq:props}
C_1\propto 1 / (\mi \omega \epsilon_0);~
C_2\propto 1\,,
\end{equation}

the matrix becomes well-conditioned and the right hand side is no longer subject to numerical loss of accuracy in the dominant terms of the solution for both plane wave and dipole excitations.

AM: this is not true, the dipole RHS still is

In fact, the resulting Loop-Star matrix $\mat{\tilde{Z}_{\Lambda \Sigma}}$ in this case has the following frequency scalings

\begin{align} \label{eq:scaling_mat_precond_ns_real}
\Re\left(\mat{\tilde{Z}_{\Lambda \Sigma}}\right) &\underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) &\mathcal{O}(\omega^2) \\
\mathcal{O}(\omega^2) & \mathcal{O}(1)
\end{bmatrix} \,, \\ \label{eq:scaling_mat_precond_ns_imag}
\Im\left(\mat{\tilde{Z}_{\Lambda \Sigma}} \right) &\underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) &\mathcal{O}(\omega) \\
\mathcal{O}(\omega) & \mathcal{O}(\omega)
\end{bmatrix}\,,
\end{align}

from which it is evident that, while the diagonal blocks of (<ref>) scale as $\mathcal{O}(1)$, all the other blocks get to zero at low frequencies, thus it results that $\mat{\tilde{Z}_{\Lambda \Sigma}}$ is well-conditioned.

The preconditioned right hand side on the other hand scales as

\begin{gather}
\Re \left(\transpose{\mat{\tilde{B}^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v_\textnormal{PW}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(1) \\
\mathcal{O}(1)
\end{bmatrix},
\Im \left(\transpose{\mat{\tilde{B}^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v_\textnormal{PW}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega)
\end{bmatrix}, \label{eq:scaling_rhs_precond_pw_ns} \\
\Re \left(\transpose{\mat{\tilde{B}^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v_\textnormal{dip}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega^4) \\
\mathcal{O}(1)
\end{bmatrix},
\Im \left(\transpose{\mat{\tilde{B}^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v_\textnormal{dip}}\right) \underrel{\omega \to 0}{=} \begin{bmatrix}
\mathcal{O}(\omega) \\
\mathcal{O}(\omega^3)
\end{bmatrix},  \label{eq:scaling_rhs_precond_dipole_ns}
\end{gather}

for plane wave and dipole excitations, respectively. From <ref>, it results that all the terms of the right hand side vectors for the plane wave and the dipole contributing to the dominants parts of the solutions are preserved at low frequencies. Therefore, the Loop-Star coefficients $C_1$, $C_2$, $C_3$, and $C_4$ would cure both matrix ill-conditioning and right hand side cancellations.

From this analysis we learn the following guidelines to produce an ansatz for the D-VIE regularizer: (i) the regularization should be a left preconditioner only, (ii) a frequency scaling proportional to $\omega^{-1}$ should be applied to the solenoidal part $\mat{P^{\Lambda}_{G_\epsilon}}$, (iii) a clear source of high-contrast breakdown originates from the ($\epsilon$-dependent) Gram matrix. By noticing that the Gram matrix $\mat{G_{\epsilon}}$ is asymptotically proportional to $\omega$, guidelines (i)-(iii) are satisfied if the solenoidal part is multiplied by $\mat{G_{\epsilon}}^{-1}$. Otherwise said, the regularizer we propose is defined as

By leveraging them, we can build preconditioners to re-scale the D-VIE with the coefficients $\mat{C}_1$, $\mat{C}_2$, $\mat{C}_3$, and $\mat{C}_4$ similarly as in (<ref>) and (<ref>), except that these coefficients are now matrices. As in (<ref>), it is sufficient to employ only $\mat{C}_1$ and $\mat{C}_2$ to re-scale the D-VIE and set $\mat{C}_3$ and $\mat{C}_4$ to $\identity$.

\begin{equation} \label{eq:precond_hc-ns}
\mat{L_{G_\epsilon}} = \gamma_\Lambda \mat{P^{\Lambda }_{G_\epsilon}} \mat{G_\epsilon}^{-1} + \gamma_\Sigma \mat{P^{\Sigma }_{G_\epsilon}} \mat{G}^{-1}\,,
\end{equation}

in which the proportionality factors $\gamma_\Lambda = {1}/{\| \mat{P^{\Lambda}_{G_\epsilon}} \|}$ and $\gamma_\Sigma = {1}/{\| \mat{P^{\Sigma}_{G_\epsilon}} \mat{G}^{-1} \mat{Z_\Phi} \mat{P^{\Sigma}_{G_\epsilon}} \|} $ are chosen to ensure that both solenoidal and non-solenoidal parts contribute with unitary weight to the final operator. Note that since the solenoidal part is re-scaled by $\mat{G_\epsilon}^{-1}$, the non-solenoidal part is re-scaled by the $\epsilon$-independent Gram matrix $\mat{G}^{-1}$ for the sake of consistency. The final formulation reads

\begin{equation} \label{eq:precond_equation_G_ns}
\mat{L_{G_\epsilon}} \mat{Z} \mat{\alpha} = \mat{L_{G_\epsilon}} \mat{v}\,.
\end{equation}

Although (<ref>) cures the low-frequency breakdown of the D-VIE by construction, it is only an ansatz regarding the cure of the HC breakdown. However, in the next section, we will prove also the effectiveness of the formulation for the HC breakdown.

§ THEORETICAL FRAMEWORK AND CONDITIONING BOUNDS

As already delineated in the previous section, given that $\gamma_\Lambda \mat{P^{\Lambda }_{G_\epsilon}} \mat{G_\epsilon}^{-1}=O(1/\omega)$ and that $ \gamma_\Sigma \mat{P^{\Sigma }_{G_\epsilon}} \mat{G}^{-1}=O(1)$, the low-frequency stability (in the limit $\omega\to 0$) follows from (<ref>)-(<ref>). Let's now focus on the high-contrast regime studying the behavior of the new preconditioned operator $\mat{M} = \mat{L_{G_\epsilon}}\mat{Z}$ in (<ref>) for high conductivity ratio. Using equations (<ref>) and (<ref>), the resulting preconditioned matrix $\mat{M}$ can be decomposed as

\begin{equation}
\begin{split}
\mat{M} &= \mat{L_{G_\epsilon}}\mat{Z}
= \gamma_\Lambda \mat{P^{\Lambda}_{G_\epsilon}} + \gamma_\Sigma \mat{P^{\Sigma}_{G_\epsilon}} \mat{G}^{-1} \mat{G_\epsilon} \\
&+ \gamma_\Lambda \mat{P^{\Lambda}_{G_\epsilon}} \mat{G_\epsilon}^{-1} \mat{Z_A}
+ \gamma_\Sigma \mat{P^{\Sigma}_{G_\epsilon}} \mat{G}^{-1}
\mat{Z_A}\\
&+ \gamma_\Sigma \mat{P^{\Sigma}_{G_\epsilon}} \mat{G}^{-1}
\mat{Z_{\Phi,1 \epsilon}}
+ \gamma_\Sigma \mat{P^{\Sigma }_{G_\epsilon}} \mat{G}^{-1}
\mat{Z_{\Phi,1 1}} \mat{P^{\Sigma }_{G_\epsilon}}\,.
\label{eq:discretized_eq_decomp_precond_ns}
\end{split}
\end{equation}

To investigate the conditioning of $\mat{M}$, we leverage the following auxiliary normalized Loop-Star decomposition matrices $\mat{Q_L}= [\mat{Q_{\Lambda L}}, \mat{Q_{\Sigma L}} ]$ and $\mat{Q_R}= [\mat{Q_{\Lambda R}}, \mat{Q_{\Sigma R}} ]$ such that

\begin{gather}
[\mat{Q_{\Lambda L}}, \mat{Q_{\Sigma L}} ] = [ \mat{G_\epsilon} \mat{\Lambda} (\transpose{\mat{\Lambda}} \mat{G_\epsilon}^2 \mat{\Lambda})^{-\frac{1}{2}} ,~ \mat{\tilde{\Sigma}} (\transpose{\mat{\tilde{\Sigma}}} \mat{\tilde{\Sigma}})^{-\frac{1}{2}}],\\
[\mat{Q_{\Lambda R}}, \mat{Q_{\Sigma R}} ] =[ \mat{\Lambda} (\transpose{\mat{\Lambda}} \mat{\Lambda})^{-\frac{1}{2}}
,~ \mat{G_\epsilon}^{-1} \mat{\tilde{\Sigma}} (\transpose{\mat{\tilde{\Sigma}}} \mat{G_\epsilon}^{-2}\mat{\tilde{\Sigma}})^{-\frac{1}{2}}],
\end{gather}

AM: You are taking the squre root of a singular matrix. Fine since 0 singular value is not here with Sigma truncated

which allow decomposing $\mat{M}$ into a block matrix $\mat{M_B}$ whose diagonal blocks correspond to the solenoidal and non-solenoidal parts of $\mat{M}$ and the off-diagonal blocks to its associated cross terms. Note that the following identities hold with these normalized decomposition matrices

\begin{gather}
\transpose{\mat{Q_{\Lambda L}}} \mat{P^{\Lambda }_{G_\epsilon}} = \transpose{\mat{Q_{\Lambda L}}};~
\transpose{\mat{Q_{\Sigma L}}} \mat{P^{\Sigma }_{G_\epsilon}} = \transpose{\mat{Q_{\Sigma L}}}\,,\\
\mat{P^{\Lambda }_{G_\epsilon}} \mat{Q_{\Lambda R}} = \mat{Q_{\Lambda R}};~
\mat{P^{\Sigma }_{G_\epsilon}} \mat{Q_{\Sigma R}} = \mat{Q_{\Sigma R}}\,.
\end{gather}

AM: it looks like you are using $\mat{\Lambda}^\mathrm{T} \mat{G_{\epsilon}} \mat{\Lambda} \left( \mat{\Lambda}^\mathrm{T} \mat{G_{\epsilon}} \mat{\Lambda} \right)^{+} = \mat I$ (among others), are we sure this is true?: Fine since 0 singular value is not here with Sigma truncated

The decomposed matrix then reads

\begin{gather}
\mat{M_B}
= \transpose{\mat{Q_L}} \mat{M} \mat{Q_R} = \begin{pmatrix}
\mat{M_{\Lambda \Lambda}} & \mat{M_{\Lambda \Sigma}}\\
\mat{M_{\Sigma \Lambda}} & \mat{M_{\Sigma \Sigma}}
\end{pmatrix},\,\text{ in which}\\
\mat{M_{\Lambda \Lambda}} = \gamma_\Lambda \transpose{\mat{Q_{\Lambda L}}} ( \identity + \mat{G_\epsilon}^{-1} \mat{Z_A}) \mat{Q_{\Lambda R}} \label{eq:diagup}\,,\\
\mat{M_{\Lambda \Sigma}} = \gamma_\Lambda \transpose{\mat{Q_{\Lambda L}}} \mat{G_\epsilon}^{-1} \mat{Z_A} \mat{Q_{\Sigma R}}\,,\\
\mat{M_{\Sigma \Lambda}} = \gamma_\Sigma \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} (\mat{G_\epsilon} + \mat{Z_A} + \mat{Z_{\Phi,1 \epsilon}} ) \mat{Q_{\Lambda R}}\,,\\
\mat{M_{\Sigma \Sigma}} = \gamma_\Sigma \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} (\mat{G_\epsilon} + \mat{Z_A} + \mat{Z_{\Phi}}) \mat{Q_{\Sigma R}}\,.\label{eq:diagdown}
\end{gather}

AM: Check the last 2 results?

Using the fact that the product of two square matrices $\mat{U}$ and $\mat{V}$ can be bounded as $\cond{\mat{U} \mat{V}} \leq \cond{\mat{U}}\cond{\mat{V}}$, we obtain the following upper bound for $\cond{\mat{M}}$

\begin{equation}
\begin{alignedat}{2}
\cond{\mat{M}} \leq \cond{\mat{Q_L}} \cond{\mat{M_B}} \cond{\mat{Q_R}}\,.\label{eq:bound_cound_M}
\end{alignedat}
\end{equation}

In Appendix <ref> we show that, in an object with $\sigma / (\omega \epsilon_0 \epsilon_r') \gg 1$, the condition number of $\mat{Q_L}$ ($\mat{A} = \mat{G_\epsilon}^{-1}$) can be bounded as

\begin{equation}
\begin{alignedat}{2}
\cond{\mat{Q_L}} \leq \left(\frac{1 + \sqrt{1 - {\|\mat{P^{\Sigma}_{\mat{G_\epsilon}}}\|^{-2}}}}{1 - \sqrt{1 - {\|\mat{P^{\Sigma}_{\mat{G_\epsilon}}}\|^{-2}}}}\right)^{\frac{1}{2}}\,.
\end{alignedat}
\end{equation}

Then, using the fact $\|\mat{P^{\Sigma}_{\mat{G_\epsilon}}}\| = \mathcal{O}(1)$ when $r_{\sigma} \rightarrow \infty$ for the geometry defined above (see Appendix <ref>), we finally obtain that the condition number of $\mat{Q_L}$ is bounded.

Similarly, the condition number of $\mat{Q_R}$ can be shown to be bounded.

Subsequently, we investigate the conditioning of $\mat{M_B}$. To this end, $\mat{M_B}$ is first split into two matrices

\begin{equation}\label{eq:bound1-block}
\mat{M_B}
= \mat{M_D} + \mat{M_O} = \begin{pmatrix}
\mat{M_{\Lambda \Lambda}} & \mat{0}\\
\mat{0} & \mat{M_{\Sigma \Sigma}}
\end{pmatrix}
\mat{0} & \mat{M_{\Lambda \Sigma}}\\
\mat{M_{\Sigma \Lambda}} & \mat{0}
\end{pmatrix},
\end{equation}

representing its diagonal blocks and its off-diagonal blocks, respectively. To avoid that the off-diagonal blocks of $\mat{M_O} $ render $\mat{M_B}$ singular, the following inequality should hold true

AM: what does it mean to perturb?

\begin{equation}\label{eq:final-bound2}
\|\mat{M_O}\| \leq \alpha s_\text{min}(\mat{M_D})\,,
\end{equation}

in which $ s_\text{min}(\mat{M_D})$ is the minimum singular value of $\mat{M_D}$ and the scalar $\alpha$ ($0<\alpha< 1$) determines how strictly the inequality should be respected. Using the fact that $
\|\mat{M_O}\| = \max(\|\mat{M_{\Lambda \Sigma}}\|,\|\mat{M_{\Sigma \Lambda}}\|) \label{eq:diag-bound}$ and $s_\text{min}(\mat{M_D}) = \min ( s_\text{min}(\mat{M_{\Lambda \Lambda}}), s_\text{min}(\mat{M_{\Sigma \Sigma}}))$, the inequality (<ref>) becomes

\begin{equation}\label{eq:final-bound}
\max(\|\mat{M_{\Lambda \Sigma}}\|,\|\mat{M_{\Sigma \Lambda}}\|) \leq \alpha \min ( s_\text{min}(\mat{M_{\Lambda \Lambda}}), s_\text{min}(\mat{M_{\Sigma \Sigma}}))\,.
\end{equation}

To know under which condition (<ref>) is satisfied, bounds for $\|\mat{M_{\Lambda \Sigma}}\|$, $\|\mat{M_{\Sigma \Lambda}}\|$, $ s_\text{min}(\mat{M_{\Lambda \Lambda}})$, and $ s_\text{min}(\mat{M_{\Sigma \Sigma}})$ are investigated next. The following result, proven in Appendix <ref>, is employed to determine a lower bound for the minimum singular value of a sum of matrices

\begin{equation} \label{eq:lb_sigmamin}
s_\text{min}(\mat{A} + \mat{B}) \geq \|\mat{A}\| \frac{\text{cond}(\mat{A}) - 1}{\text{cond}(\mat{A})^2} \frac{ s_\text{min}(\mat{A}) - \| \mat{B} \| }{s_\text{min}(\mat{A}) + \| \mat{B} \| }\,,
\end{equation}

it applies when $\mat{A}$ is an invertible matrix and $\|\mat{B} \| < s_\text{min}(\mat{A})$.

Before employing this identity on the diagonal blocks of $\mat{M_D}$, we introduce the following matrices

\begin{gather}
\mat{A_\Lambda} = \gamma_\Lambda \transpose{\mat{Q_{\Lambda L}}} \mat{Q_{\Lambda R}}\,,\\
\mat{B_\Lambda} = \gamma_\Lambda \transpose{\mat{Q_{\Lambda L}}} \mat{G_\epsilon}^{-1} \mat{Z_A} \mat{Q_{\Lambda R}}\,,\\
\mat{A_\Sigma} = \gamma_\Sigma \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \mat{Z_\Phi} \mat{Q_{\Sigma R}}\,,\\
\mat{B_\Sigma} = \gamma_\Sigma \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} (\mat{G_\epsilon} + \mat{Z_A})\mat{Q_{\Sigma R}}\,,
\end{gather}

which represent the terms in the diagonal blocks of $\mat{M_D}$ ((<ref>) and (<ref>)). Then, using the fact that $\mat{A_\Lambda}$ and $\mat{A_\Sigma}$ are invertible ($\mat{G}$ and $\mat{Z_\Phi}$ being non-singular) and restricting ourselves to the case in which

\begin{equation}
s_\text{min}(\mat{A_\Lambda}) > \|\mat{B_\Lambda}\| \label{eq:cond_diagUp};~
s_\text{min}(\mat{A_\Sigma}) > \|\mat{B_\Sigma}\|\,,
\end{equation}

AM: enforcing how? By playing on the $\gamma$ ? If so it is not clear, especially since they popped up out of nowhere

the minimum singular values of $\mat{M_{\Lambda \Lambda}}$ and $\mat{M_{\Sigma \Sigma}}$ can be bounded as

\begin{gather}
s_\text{min}(\mat{M_{\Lambda \Lambda}})
\geq f(\mat{A_\Lambda}) \frac{s_\text{min}(\mat{A_\Lambda}) - \| \mat{B_\Lambda} \|}{s_\text{min}(\mat{A_\Lambda}) + \| \mat{B_\Lambda} \| }\label{eq:lb_blockL1}\,,\\
s_\text{min}(\mat{M_{\Sigma \Sigma}}) \geq f(\mat{A_\Sigma}) \frac{s_\text{min}(\mat{A_\Sigma}) - \| \mat{B_{\Sigma}} \|}{ s_\text{min}(\mat{A_\Sigma}) + \| \mat{B_{\Sigma}} \|}\,,\label{eq:lb_blockS1}
\end{gather}

in which $f(\mat{A}) = \|\mat{A}\| (\text{cond}(\mat{A}) - 1)/(\text{cond}(\mat{A})^2)$. Next, it is shown in Appendix <ref> that the minimum and maximum singular values of the scaled Gram matrix $\mat{G_\epsilon}$ ($s_\text{min}(\mat{G_\epsilon})$ and $\|\mat{G_\epsilon}\|$, respectively) can be bounded as $s_\text{min}(\mat{G_\epsilon}) = \| \mat{G_\epsilon}^{-1} \|^{-1} \geq { s_\text{min}(\mat{G}) \omega}/{\sigma_\text{max}}$ and $\|\mat{G_\epsilon}\| \leq {\|\mat{G}\| \omega}/{\sigma_\text{min}}$. Using these bounds and the fact that $\mat{Z_A}$ is only frequency-dependent when $\kappa \approx 1$ (high-contrast regime), we can extract the conductivity and frequency dependence from $\mat{B_\Lambda}$ and $\mat{B_\Sigma}$ AM: is this enough B lambda and B simgma contain Z A too in (<ref>) and (<ref>). The norm of the matrix $\mat{B_\Lambda}$ can be bounded as

\begin{equation} \label{eq:bound_BLambda}
\begin{split}
\|\mat{B_\Lambda} \| &\leq
\gamma_\Lambda \|\transpose{\mat{Q_{\Lambda L}}}\| \|\mat{G_\epsilon}^{-1}\| \|\mat{Z_A} \mat{Q_{\Lambda R}}\| \\
&\leq \gamma_\Lambda k_0^2 \frac{\sigma_\text{max} \|\mat{G}^{-1}\|}{\omega} \| \mat{Z^1_A} \mat{Q_{\Lambda R}}\|\,,
\end{split}
\end{equation}

AM: you are using a result on $\|\mat{G_\epsilon}^{-1}\|$ that is stated nowhere

where $\mat{Z_A} = k_0^2 \mat{Z^1_A}$. Then, $\|\mat{B_\Sigma}\|$ can be bounded from above as

\begin{equation} \label{eq:bound_BSigma}
\begin{split}
\|\mat{B_\Sigma} \| & \leq \gamma_\Sigma \left(\| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \mat{G_\epsilon} \mat{Q_{\Sigma R}}\| + \| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \mat{Z_A}\mat{Q_{\Sigma R}}\|\right) \\
&\leq \gamma_\Sigma \left( \|\transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \|\frac{\|\mat{G}\| \omega}{\sigma_\text{min}} + k_0^2\| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \mat{Z^1_A}\mat{Q_{\Sigma R}}\|\right).
\end{split}
\end{equation}

Leveraging (<ref>), (<ref>), (<ref>), (<ref>), and the fact that the inequalities $\|\mat{U}\| \leq c$ ($c \in \R_+$), $s \geq ({a - \|\mat{U}\|})/({a + \|\mat{U}\|})$ ($a \in \R_+^*$), and $\|\mat{U}\| < a$ imply $  s \geq {(a - c)}/{(a + c)}$, the following lower bounds for the minimum singular values of $\mat{M_{\Lambda \Lambda}}$ and $\mat{M_{\Sigma \Sigma}}$, which are a function of $r_\sigma$ and $\omega$, are obtained

\begin{gather}
s_\text{min}(\mat{M_{\Lambda \Lambda}})
\geq f(\mat{A_\Lambda}) \frac{s_\text{min}(\mat{A_\Lambda}) - r_\sigma \omega \tau_\Lambda}{s_\text{min}(\mat{A_\Lambda}) + r_\sigma \omega \tau_\Lambda }\label{eq:lb_blockL}\,,\\
s_\text{min}(\mat{M_{\Sigma \Sigma}}) \geq f(\mat{A_\Sigma}) \frac{s_\text{min}(\mat{A_\Sigma}) - \omega^2 \tau_{\Sigma A} - \omega \tau_{\Sigma G} }{ s_\text{min}(\mat{A_\Sigma}) + \omega^2\tau_{\Sigma A} + \omega \tau_{\Sigma G} }\,, \text{ in which}\label{eq:lb_blockS} \\
\tau_\Lambda = \gamma_\Lambda \sigma_\text{min} / c_0^2 \|\mat{G}^{-1}\| \| \mat{Z^1_A} \mat{Q_{\Lambda R}}\|\,, \\
\tau_{\Sigma A} =\gamma_\Sigma / c_0^2 \| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \mat{Z^1_{A}} \mat{Q_{\Sigma R}}\| \,,\\
\tau_{\Sigma G} =\gamma_\Sigma \sigma_\text{min}^{-1} \| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \| \|\mat{G}\|\,,
\end{gather}

are parameters that do not depend of $\omega$ and $r_\sigma$ ($\kappa \approx 1$ in this regime).

Provided that (<ref>) holds, we thus obtain lower bounds for the minimum singular values of $\mat{M_{\Lambda \Lambda}}$ and $\mat{M_{\Sigma \Sigma}}$ ((<ref>) and (<ref>)) AM: you need an argument to use the bounds in the orignial inequality (equlity inside equlatity) which translates into the non-singularity of $\mat{M_D}$ in (<ref>). In that case, identity (<ref>) can be applied to (<ref>) to find a lower bound for the minimum singular value of $\mat{M_B}$. Before applying (<ref>), we also need to enforce that the norms of $\mat{M_{\Lambda \Sigma}}$ and $\mat{M_{\Sigma \Lambda}}$ and the minimum singular values of $\mat{M_{\Lambda \Lambda}}$ and $\mat{M_{\Sigma \Sigma}}$ respect inequality (<ref>). Using a similar approach as before, the off-diagonal blocks $\mat{M_{\Lambda \Sigma}}$ and $\mat{M_{\Sigma \Lambda}}$ can be bounded as

\begin{gather}
\|\mat{M_{\Lambda \Sigma}}\| \leq r_\sigma \omega \tau_{\Lambda \Sigma A} \label{eq:ub-blockL}\,,\\
\|\mat{M_{\Sigma \Lambda}}\| \leq \omega \tau_{\Sigma G} + \omega^2 \tau_{\Sigma \Lambda A} + \omega \tau_{\Sigma \Lambda \Phi}\,, \text{ in which}\label{eq:ub-blockS} \\
\tau_{\Lambda \Sigma A} = \gamma_\Lambda \sigma_\text{min} / c_0^2 \|\mat{G}^{-1}\| \| \mat{Z^1_A} \mat{Q_{\Sigma R}}\| \,,\\
\tau_{\Sigma \Lambda A} = \gamma_\Sigma / c_0^2 \| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \mat{Z^1_{A}} \mat{Q_{\Lambda R}}\| \,,\\
\tau_{\Sigma \Lambda \Phi} =\gamma_\Sigma \sigma_\text{min}^{-1} \| \transpose{\mat{Q_{\Sigma L}}} \mat{G}^{-1} \| \|\mat{Z^1_{\Phi,1\epsilon}}\|\,,
\end{gather}

such that $\|\mat{Z_{\Phi,1\epsilon}}\| \leq \omega/\sigma_\text{min}\|\mat{Z^1_{\Phi,1\epsilon}}\|$ in which $\mat{Z^1_{\Phi,1\epsilon}}$ is a permittivity-independent matrix that follows the relation $\mat{Z_{\Phi,1\epsilon}} = \mat{Z^1_{\Phi,1\epsilon}} \mat{\Upsilon_{\deltaKappa}}$ with $\mat{\Upsilon_{\deltaKappa}}$ being a diagonal matrix filled as $[\mat{\Upsilon_{\deltaKappa}}]_{nn} = \epsilon_0^{-1} \deltaKappa_n$. This last inequality can be obtained by observing that $\|\mat{\Upsilon_{\deltaKappa}}\| \leq \omega / \sigma_{\text{min}} $.

Finally, we can enforce the conditions (<ref>) and (<ref>) using the bounds derived in (<ref>), (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>). These conditions being enforced, the resulting condition number of $\mat{M_B}$ can be bounded from above as

\begin{equation} \label{eq:cond_MB}
\begin{split}
\cond{\mat{M_B}}
&\leq \cond{\mat{M_D}\left( \identity + \mat{M_D}^{-1} \mat{M_O} \right)} \\
&\leq \cond{\mat{M_D}} \frac{1 + \|\mat{M_D}^{-1} \mat{M_O}\|}{1 - \|\mat{M_D}^{-1} \mat{M_O}\|} \\
&\leq \frac{\max(\|\mat{M_{\Lambda \Lambda}}\|,\|\mat{M_{\Sigma \Sigma}}\|)}{\min(s_\text{min}(\mat{M_{\Lambda \Lambda}}),s_\text{min}(\mat{M_{\Sigma \Sigma}}))} \frac{1 + \alpha}{1 - \alpha}\,,
\end{split}
\end{equation}

in which the condition (<ref>) was employed. An upper bound for the condition number of $\mat{M}$ can then be obtained from (<ref>) using (<ref>).

The variables in these bounds are the frequency and the maximum conductivity ratio in the object $r_\sigma$ which are dictated by the electromagnetic problem to solve, and the parameter $\alpha$ that should be chosen depending on the range of validity requirement.

AM: $r_\sigma$ is not really free, is it? It is dictated by the scenario
AM: this thing on $\alpha$ is still not clear: Fixed

The values of $\omega$, $r_\sigma$, and $\alpha$ for which the bounds above-mentioned are respected define the range in which the formulation is well behaving. Note that the parameter $\alpha$ allows tuning simultaneously the bound for the condition number of $\mat{M}$ and the range of values of $\omega$ and $r_\sigma$ for which the conditions are met, e.g., a low value of $\alpha$ results in a smaller upper bound for the condition number in (<ref>) but it also makes the condition (<ref>) harder to respect.AM: I fear this is not super clear...

The bounds previously derived establish that, in the static limit, the conditions (<ref>) and (<ref>) are met regardless of the maximum conductivity ratio $r_\sigma$ and of the parameter $\alpha$ ($0<\alpha< 1$) and the upper bound for the condition number of $\mat{M_B}$ given in (<ref>) also becomes independent of $r_\sigma$ and $\alpha$. Moreover, since the condition number of $\mat{Q_L}$ and $\mat{Q_R}$ is bounded when $r_\sigma \rightarrow \infty$, it results from (<ref>) that in the static limit the condition number of $\mat{M}$ is bounded as the maximum conductivity ratio increases. Away from that limit, $\mat{M}$ will still be properly conditioned as long as the bounds above-mentioned are enforced. To show the proper behavior of this formulation beyond the static limit, some values of $r_\sigma$ and the frequency for which the bounds are respected are given in the numerical results section. Finally, the reader should note that at high-frequency, when low-frequency regularization is not required, the standard D-VIE can be used without changing the implementation by simply setting the coefficients and matrices multiplying the projectors in (<ref>) to 1 and $\mat{\identity}$, respectively.

§ IMPLEMENTATION RELATED DETAILS

This section presents the details related to the implementation of the proposed preconditioner. Starting from an existing D-VIE, the quasi-Helmholtz projectors must be computed efficiently not to deteriorate the overall complexity of the solver. The main difficulty to build the projectors resides in the inversion of the stiffness matrix $\transpose{\mat{\Sigma}} \mat{G_\epsilon}^{-1} \mat{\Sigma}$, which has a condition number that grows when the discretization of the geometry is refined. To remedy this issue, this operation can be done leveraging the algebraic multigrid (AMG) method as preconditioner <cit.> together with an iterative solver. Besides, since the Gram matrix $\mat{G_\epsilon}$ is sparse and symmetric, its direct inverse can be obtained in an efficient manner using a multifrontal solver <cit.>.

To address the low-frequency breakdown in surface formulations such as the electric field integral equation (EFIE), the construction of the projectors requires inverting the graph Laplacian matrix $\mat{\Sigma_{RWG}}^\mathrm{T} \mat{\Sigma_{RWG}}$ in which $\mat{\Sigma_{RWG}}$ denotes the surface star-to-RWG transformation matrix <cit.>. This operation can be performed in linear complexity by leveraging the algebraic multigrid (AMG) method as preconditioner <cit.> together with an iterative solver.

However, to solve the high-constrast and the low-frequency problems in the D-VIE, we use scaled projectors, which involve the inversion of a weighted graph Laplacian. Since the scatterer is inhomogeneous, $\mat{G_{\epsilon}}$ is ill-conditioned and its inverse cannot always be computed efficiently.
[inline]CH: @Francesco: Taking the absolute value of $\mathrm{diag}(\mat{G_{\epsilon}})$ as $\mat{D}$ seems to work similarly for the condition number. However, the reason behind the current definition of $\mat{D}$ is to avoid numerical cancellations in the precond RHS when the object is conductive (use of imaginary part).

In addition to this, however, we have found also a very effective strategy that does not require the inversion of $\mat{G_{\epsilon}}$, resulting in less computational and implementation efforts. In fact, we observed numerically that using the real part or the imaginary part of the diagonal of the Gram matrix in the stiffness matrix also makes the formulation stable. The set of projectors in this case is defined as follows

\begin{gather}
\mat{P^{\Sigma}_D} = \mat{D}^{-1} \mat{\Sigma} (\mat{\Sigma}^\mathrm{T} \mat{D}^{-1} \mat{\Sigma})^{+} \mat{\Sigma}^\mathrm{T}\,,\\
\mat{P^{\Lambda}_D} = \identity - \mat{P^{\Sigma}_D} = \mat{\Lambda} (\mat{\Lambda}^\mathrm{T} \mat{D} \mat{\Lambda})^{-1} \mat{\Lambda}^\mathrm{T} \mat{D}\,,
\end{gather}

where $\mat{D} \in \Ci^{N_F \times N_F}$ is a diagonal matrix constructed as

\begin{equation}
\mat{D}=
\begin{cases} \phantom{\mi} \Re \left(\mat{D_{G}} \right) & \text{ if } \|\Re \left(\mat{D_{G}}^{-1} \right)\| > \|\Im \left( \mat{D_{G}}^{-1} \right) \|\\
\mi \Im\left(\mat{D_{G}}\right) & \text{ if } \|\Re \left(\mat{D_{G}}^{-1} \right)\| < \|\Im \left( \mat{D_{G}}^{-1} \right) \|\,,
\end{cases}
\label{eq:def_D}
\end{equation}

with $\mat{D_{G}}$ being a diagonal matrix filled with the diagonal of $\mat{G_{\epsilon}}$. Similarly as in (<ref>), the preconditioned D-VIE with $\mat{P^{\Sigma}_D}$ and $\mat{P^{\Lambda}_D}$ reads

\begin{gather}
\mat{L_D} \mat{Z} \mat{\alpha} = \mat{L_D} \mat{v}\,, \text{ in which} \label{eq:precond_equation_D}\\
\mat{L_D} = \frac{\mat{P^{\Lambda}_D} \mat{D}^{-1}}{\| \mat{P^{\Lambda}_D} \mat{D}^{-1} \mat{G_\epsilon} \mat{P^{\Lambda}_D} \|} \ + \frac{\mat{P^{\Sigma}_D} \mat{D_0}^{-1}}{\| \mat{P^{\Sigma}_D} \mat{D_0}^{-1} \mat{Z_\Phi} \mat{P^{\Sigma}_D} \|}\,,
\end{gather}

and $\mat{D_0}$ is a diagonal matrix filled with the diagonal of $\mat{G}$.

This new preconditioner gives rise to a slightly higher condition number than with the preconditioner obtained from $\mat{G_{\epsilon}}$, but has been numerically verified to be stable in a broad frequency range and for high permittivity objects, as illustrated in next section.
Moreover the stiffness matrix $\mat{\Sigma}^\mathrm{T} \mat{D}^{-1} \mat{\Sigma}$, which is also a weighted graph Laplacian matrix, can be inverted efficiently using an aggregation-based AMG method from <cit.> together with a conjugate gradient (CG) algorithm.

§ OLD STUFF

Besides the loop/star-to-SWG mapping matrix, a new form of volume quasi-Helmholtz projectors can be introduced from (<ref>). This new set of projectors, $\mat{P^{\Sigma}_{A^{-1}}}$ and $\mat{P^{\Lambda}_{A^{-1}}}$, can respectively select the non-solenoidal and the solenoidal parts of the discretized operator without being limited by the conditioning issue of the mapping matrix $\mat{B^{A^{-1}}_{\Lambda \Sigma}}$. They are defined as

\begin{gather}
\mat{P^{\Sigma}_{A^{-1}}} = \mat{A} \mat{\Sigma} \left( \mat{\Sigma}^\mathrm{T} \mat{A} \mat{\Sigma} \right)^{+} \mat{\Sigma}^\mathrm{T} \label{eq:scaled-projectors-star-A-ns}\,,\\
\mat{P^{\Lambda}_{A^{-1}}} = \identity - \mat{P^{\Sigma}_{A^{-1}}} = \mat{\Lambda} \left( \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \mat{\Lambda} \right)^{-1} \mat{\Lambda}^\mathrm{T} \mat{A}^{-1} \label{eq:scaled-projectors-loop-A-ns} ,
\end{gather}

where $^+$ is the Moore-Penrose pseudo inverse. Differently from the definitions of quasi-Helmholtz projectors in previous works <cit.>, we introduce here the invertible and symmetric scaling matrix $\mat{A}$ in the projectors. Given this new definition, identity (<ref>) does not follow directly from the properties of standard projectors, but its proof is provided in Appendix <ref>.

This behavior will then be regularized using the volume quasi-Helmholtz projectors $\mat{P^{\Lambda}_{A^{-1}}}$ and $\mat{P^{\Sigma}_{A^{-1}}}$ in combination with appropriate scalings determined during the Loop-Star analysis to regularize the D-VIE at low frequencies without deteriorating its dense discretization behavior. To benefit from the cancellations of identity (<ref>) in the Loop-Star decomposition of the D-VIE and (<ref>) to (<ref>) for its projector counterpart, the scaling matrix $\mat{A}$ could be set to $\mat{G_\epsilon}^{-1}$.

, we obtain the full rank loop/star-to-SWG basis $\mat{B^{\identity}_{\Lambda \Sigma}} = \left[\mat{\Lambda}~\mat{\Sigma} \right]\in \R^{N_{F} \times N_F}$. Applying $\mat{B^{\identity}_{\Lambda \Sigma}}$ to the discretized D-VIE results in the following 2-by-2 block matrix

\begin{equation} \label{eq:ls_decomp_bg}
\transpose{\mat{B^{\identity}_{\Lambda \Sigma}}} \mat{Z} \mat{B^{\identity}_{\Lambda \Sigma}} =\begin{bmatrix}
\mat{\Lambda}^\mathrm{T}\mat{Z} \mat{\Lambda} & \mat{\Lambda}^\mathrm{T}\mat{Z} \mat{\Sigma} \\
\mat{\Sigma}^\mathrm{T}\mat{Z} \mat{\Lambda} & \mat{\Sigma}^\mathrm{T} \mat{Z} \mat{\Sigma}
\end{bmatrix}.
\end{equation}

In this loop/star-to-SWG block matrix, the solenoidal and non-solenoidal parts are separated and can be re-scaled appropriately to cure the ill-conditioning of the original matrix. However, the matrix $\mat{B^{\identity}_{\Lambda \Sigma}}$ is ill-conditioned with a condition number that increases when the discretization of the geometry is refined.

The final equation for the normalized Loop-Star D-VIE reads

\begin{equation} \label{eq:ls-DVIE}
\transpose{\mat{\tilde{B}^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{Z} \mat{B^{\tilde{G}_\epsilon}_{\Lambda \Sigma}} \mat{\alpha_{\Lambda \Sigma}} = \transpose{\mat{\tilde{B}^{\tilde{G}_\epsilon}_{\Lambda \Sigma}}} \mat{v}.
\end{equation}

We note that, since the computation of the square roots present in the definition of the normalized loop/star-to-SWG matrix is not efficient, the non-normalized version of this mapping matrix should be employed in a numerical implementation of the Loop-Star D-VIE.

Subsequently, we can logically relate the above-mentioned coefficients and the quasi-Helmholtz projectors defined in (<ref>) and (<ref>) with $\mat{A} = \mat{\tilde{G}_\epsilon}^{-1}$ to build a preconditioned equation that does not suffer from the dense discretization breakdown. We obtain the following left and right preconditioners

\begin{gather}
\mat{L} = C_1 \mat{P^{\Lambda \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} + C_2 \mat{P^{\Sigma \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} \label{eq:preconditioner_left_definition_ns}\,,\\
\mat{R} = C_3 \mat{P^\Lambda_{\mat{\tilde{G}_\epsilon}}} + C_4 \mat{P^\Sigma_{\mat{\tilde{G}_\epsilon}}}.
\label{eq:preconditioner_right_definition_ns}
\end{gather}

Since $C_3=C_4=1$, we have $\mat{R} = \identity$, and the resulting preconditioned equation to be solved reads

\begin{equation} \label{eq:regularized_sys_proj_lf_ns}
\mat{L} \mat{Z} \mat{\alpha} = \mat{L} \mat{v},
\end{equation}

in which

\begin{equation}
\begin{split}
\mat{L} &= C_1 \gamma_\Lambda\mat{P^{\Lambda \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} + C_2 \gamma_\Sigma \mat{P^{\Sigma \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} \label{eq:preconditioner_left_ns} \\
&= \frac{ \gamma_\Lambda}{\mi \omega \epsilon_0} \mat{P^{\Lambda \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} + \gamma_\Sigma \mat{P^{\Sigma \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}}.
\end{split}
\end{equation}

Note that the frequency-independent coefficients $\gamma_\Lambda$ and $\gamma_\Sigma$ were added in (<ref>) to make the norms of $\gamma_\Lambda/(\mi \omega \epsilon_0) \mat{P^{\Lambda \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} \mat{Z} \mat{P^{\Lambda}_{\mat{\tilde{G}_\epsilon}}}$ and $ \gamma_\Sigma \mat{P^{\Sigma \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} \mat{Z} \mat{P^{\Sigma}_{\mat{\tilde{G}_\epsilon}}}$ of order $\mathcal{O}(1)$ regardless of the geometry. They are defined as

\begin{gather}
\gamma_\Lambda = ({\| \epsilon_0^{-1} \mat{G} \mat{P^{\Lambda}_{\mat{\tilde{G}_\epsilon}}} \|})^{-1} \label{eq:coeff_lambda_1_ns}\,,\\
\gamma_\Sigma = ({\| \mat{P^{\Sigma \mathrm{T}}_{\mat{\tilde{G}_\epsilon}}} \mat{Z_\Phi} \mat{P^{\Sigma}_{\mat{\tilde{G}_\epsilon}}} \|})^{-1},
\label{eq:coeff_sigma_1_ns}
\end{gather}

with $\norm {\cdot} $ denoting the $L^2\text{-norm}$ and $\left[\mat{G}\right]_{nm} = \langle \vt{f_n},\vt{f_m} \rangle_\Omega$ being the SWG-SWG Gram matrix.

§ NUMERICAL RESULTS

To further corroborate the theoretical developments, the new formulation has been tested in several scenarios.
Both regularization with projectors scaled with the full Gram matrix (equation (<ref>), referred in the following as “Regularized D-VIE $\mat{G_{\epsilon}}$”) and with projectors scaled with the diagonal matrix $\mat{D}$ (equation (<ref>), referred in the following as “Regularized D-VIE $\mat{D}$”) will be considered. These two regularized formulations will be of course compared to the standard D-VIE but also, for the sake of completeness, to a Loop-Star D-VIE which can be obtained by selecting $\mat{A} = \identity$ in (<ref>) and rescaling the solenoidal and non-solenoidal parts of $\mat{B^{\identity}_{\Lambda \Sigma}}$ by the coefficients introduced in (<ref>), respectively. Note that this Loop-Star decomposition is the direct extension of the decomposition employed for surface formulations <cit.>.

In the first test, the geometry used is composed of 3 homogeneous concentric spheres (i.e. three layer head model) with radii $\SI{87}{\milli \meter}$, $\SI{92}{\milli \meter}$, and $\SI{100}{\milli \meter}$ and respective normalized conductivities $\sigma_1=1$, $\sigma_2=1/15$, and $\sigma_3 = 1$, which represent the conductivities of brain, skull, and scalp in the quasi-static regime <cit.>. First, we verify the conditioning of the system matrices of the two regularized D-VIE, the Loop-Star D-VIE, and the standard D-VIE as a function of the frequency (Fig. <ref>a). The condition number of the standard D-VIE grows as the frequency decreases while the condition number of the two regularized D-VIE and the Loop-Star D-VIE remains constant until very low frequencies, which confirms the curative effects of our preconditioners. It should also be noted that, although both are constant, the condition number of the formulations we propose here is much lower than the condition number of the Loop-Star scheme, as expected.

xlabel style=font=,
xlabel=Frequency (Hz),
ylabel style=font=,
ylabel=Condition number,
grid style=line width=.1pt, draw=gray!10,
legend style=at=(0.02,1.2), anchor=south west, legend cell align=left, align=left, draw=white!15!black,
legend columns=-1

[line width=1.5pt, color=mycolor2, mark=o, mark options=solid, mycolor2]
table[row sep=crcr]
1.77827941003892e-07	1.9851e+18

0.00316227766016838	1.7645e+13

56.2341325190349	 9.9224e+08

1000000	5.5798e+04

Standard D-VIE

[line width=1.5pt, color=mycolor4, mark=x, mark options=solid, mycolor4]
table[row sep=crcr]
1e-40	3.619505340387529e+04

1e-27	3.619505340387529e+04

1e-14	3.619505340387529e+04

1e-11	3.619505340387529e+04

1.77827941003892e-07	3.619505340387529e+04

0.00316227766016838	3.619505340387529e+04

56.2341325190349	3.619505340387529e+04

1000000	3.619505340387529e+04

Loop-star D-VIE

[line width=1.5pt, color=mycolor1, mark=+, mark options=solid, mycolor1]
table[row sep=crcr]
1e-40	15.9453860554636

1e-27	15.8434161763263

1e-14	15.9337209991274

1e-11	15.8789357843667

1.77827941003892e-07	15.8146017424083

0.00316227766016838	15.9076613430996

56.2341325190349	15.8356313235237

1000000	15.92649919883

Regularized D-VIE $\mat{G_\epsilon}$

[line width=1.5pt, color=mycolor3, mark=diamond,dashed, mark options=solid, mycolor3]
table[row sep=crcr]
1e-40	16.7343838443863

1e-27	17.0979690592484

1e-14	18.5413733805133

1e-11	16.8563856484957

1.77827941003892e-07	17.3474147721712

0.00316227766016838	17.282564198998

56.2341325190349	16.6918818053481

1000000	17.0803533596877

Regularized D-VIE $\mat{D}$

Low frequency conditioning of the standard D-VIE, the Loop-Star D-VIE, and the two proposed D-VIE computed on a 3-layer conductive sphere.