# On the Use of Quasi-Helmholtz Projectors for Modeling Structures With Junctions at Arbitrarily Low Frequencies

Johann Bourhis, Adrien Merlini, Francesco P Andriulli

► To cite this version:

Johann Bourhis, Adrien Merlini, Francesco P Andriulli. On the Use of Quasi-Helmholtz Projectors for Modeling Structures With Junctions at Arbitrarily Low Frequencies. IEEE Transactions on Antennas and Propagation, 2024, 72, pp.8097 - 8102. 10.1109/tap.2024.3439772. hal-05073217

HAL Id: hal-05073217

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

Submitted on 19 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Communication

## On the Use of Quasi-Helmholtz Projectors for Modeling Structures With Junctions at Arbitrarily Low Frequencies

Johann Bourhis¹, Adrien Merlin¹, and Francesco P. Andriulli¹

**Abstract**—This work investigates the usage of quasi-Helmholtz projectors to cure the low frequency breakdown of the electric field integral equation (EFIE) when modeling structures containing junctions. We show that the quasi-Helmholtz projectors can still efficiently perform projection into the solenoidal and nonsolenoidal subspaces in the presence of junctions. Preconditioning strategies leveraging these projectors are capable of curing the low-frequency breakdown that severely compromises the accuracy and conditioning of the EFIE. Numerical studies validate the use of algebraic multigrid approaches for achieving fast and efficient implementation of the proposed scheme for problems with multiple junctions. The effectiveness of the scheme and its applicability to problems of industrial relevance is illustrated through a series of numerical results.

**Index Terms**—Boundary element methods (BEMs), graph Laplacian, junctions, low-frequency preconditioning, quasi-Helmholtz projectors.

## I. INTRODUCTION

Time harmonic electromagnetic scattering and radiation problems for perfectly electrically conducting (PEC) objects are often solved using the electric field integral equation (EFIE) when spurious resonances cannot arise. Unfortunately, as the modeling frequency decreases, this formulation is subject to increasingly severe ill-conditioning and yields increasingly less accurate results because of loss of significance in floating point operations [1]. A classical approach to fixing the low-frequency breakdown is to leverage quasi-Helmholtz decompositions [1] to treat independently the solenoidal and nonsolenoidal contributions of the elements of the discrete system. This can be achieved, for instance, using Loop-Star decompositions [1], [2], [3], [4]. However, this approach worsens the conditioning of the system for dense meshes and requires the identification of the global cycles in the structure, which is often a burdensome operation [1]. This burden can significantly worsen in the presence of junctions since they can give rise to a significant number of global cycles associated with apertures [3], [5].

The quasi-Helmholtz projectors [4] exploit orthogonal projections over the Loop and Star subspaces to decompose and rescale the problematic contributions of the system matrices, solution, and right-hand side. This technique does not require an explicit identification

of the global cycles and does not worsen the dense discretization conditioning of the preconditioned equations but its generalization to the presence of junctions had not been investigated yet.

The contribution of this communication is twofold.

1. The usage of quasi-Helmholtz projector strategies for a structure containing junctions is investigated, addressing the questions of completeness and orthogonality of the underlying decompositions.

2. The computational cost of the approach is analyzed both by studying the effect of the presence of junctions on the inversion of the graph Laplacian that appears in the definition of the projectors and by verifying the overall numerical performance of the scheme.

In particular, we study numerically the performance of the aggregation-based algebraic multigrid (AGMG) [6], [7], which is one possible implementation of multigrid methods, well-known for providing an efficient preconditioner for graph Laplacian systems [8], [9]. Our observations show that the effectiveness of the quasi-Helmholtz projectors approach is not jeopardized in the presence of junctions and that it can be successfully used with realistic structures of industrial relevance.

This communication is organized as follows. Section II lays out the background and notations. Section III analyzes the quasi-Helmholtz projectors in the context of junctions, addressing their orthogonality, and the spectral properties of their associated graph Laplacians. Section IV provides numerical results to validate all the strategies presented here. Finally, Section V presents the conclusion of the work. Preliminary results from this communication were presented in the conference contribution [10].

## II. BACKGROUND AND NOTATION

Consider the problem of determining the scattered electric field $E^{\text{sct}}$ radiated when an incident electric field $E^{\text{inc}}$ impinges on a PEC obstacle. We denote by $\Gamma$ the boundary of the obstacle and by $\eta$ the characteristic impedance of the exterior medium. We assume the fields to be time-harmonic with wavenumber $k$. The electric current density $J$ induced over $\Gamma$ by $E^{\text{inc}}$ that radiates $E^{\text{sct}}$ can be obtained by solving the EFIE

$$ \mathcal{T} \mathbf{J} = -\hat{\mathbf{n}} \times \frac{1}{\eta} \mathbf{E}^{\text{inc}} \quad (1) $$

where $\mathcal{T} = jk\mathcal{T}_s - (jk)^{-1}\mathcal{T}_h$, with

$$ (\mathcal{T}_s \mathbf{J})(\mathbf{r}) = \hat{\mathbf{n}} \times \int_{\Gamma} G(\mathbf{r}, \mathbf{r}') \mathbf{J}(\mathbf{r}') dS(\mathbf{r}') \quad (2) $$

$$ (\mathcal{T}_h \mathbf{J})(\mathbf{r}) = \hat{\mathbf{n}} \times \nabla_{\Gamma} \int_{\Gamma} G(\mathbf{r}, \mathbf{r}') \nabla_{\Gamma} \cdot \mathbf{J}(\mathbf{r}') dS(\mathbf{r}') \quad (3) $$

Manuscript received 4 January 2024; revised 20 June 2024; accepted 28 July 2024. Date of publication 13 August 2024; date of current version 8 October 2024. This work has received funding from the European Union's Horizon 2020 (EU H2020) research and innovation programme through the Marie Skłodowska-Curie Actions (MSCA) under grant agreement n°955476 (project COMPETE), and from the EU H2020 through the European Research Council (ERC) under grant agreement n°724846 (project 321). (Corresponding author: Francesco P. Andriulli.)

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

Adrien Merlin is with the Microwave Department, IMT Atlantique, 29285 Brest, France (e-mail: adrien.merlini@imt-atlantique.fr).

Color versions of one or more figures in this article are available at https://doi.org/10.1109/TAP.2024.3439772.

Digital Object Identifier 10.1109/TAP.2024.3439772

in which $\hat{\mathbf{n}}$ is the outgoing normal vector along $\Gamma$ and the Green's function is $G(\mathbf{r}, \mathbf{r}') = e^{-jk|\mathbf{r}-\mathbf{r}'|} / (4\pi|\mathbf{r}-\mathbf{r}'|)$.

The EFIE (1) can be solved numerically via the boundary element method (BEM) by expanding the current as a linear combination of $N$ Rao-Wilton-Glisson (RWG) div-conforming basis functions [11],
<figure><img src="image_2.png" /><figcaption>Fig. 1. Notations of the vertices and cells used to define the RWGs.</figcaption></figure>

denoted $\varphi_n$, $n = 1, \dots, N$

$$J(r) \approx \sum_{n=1}^{N} [j]_n \varphi_n(r). \quad (4)$$

These functions are defined on a triangular mesh of $\Gamma$ as

$$\varphi_n(r) = \begin{cases} \pm \frac{r - r_n^{\pm}}{2A_n^{\pm}}, & \text{if } r \in c_n^{\pm} \\ 0, & \text{otherwise} \end{cases} \quad (5)$$

where the notation for the geometrical elements describing a pair of adjacent triangle cells $c_n^{\pm}$ is given in Fig. 1 and $A_n^{\pm}$ denotes their respective areas. Note that the definition above does not include edge normalization.

The EFIE can be used for both open and closed manifolds that can contain junctions, i.e., structures for which the associated mesh contains edges shared by more than two cells (see for instance Figs. 2 and 3 and [12]). A popular approach to handle structures with junctions is to associate to a mesh edge $e$ which connects together $n_e$ cells, $n_e - 1$ independent RWG functions [12], [13]. These functions are defined by considering $n_e - 1$ independent couples of cells sharing the edge [13]. Different choices of configurations that define RWG functions spanning the same space are possible, but they lead to different connectivity properties between cells. Here and in the rest of this contribution we consider that two adjacent cells $c_m$ and $c_n$ are connected if there exists an RWG function $\varphi_i$ such that $c_i^+ = c_m$ and $c_i^- = c_n$ or $c_i^+ = c_n$ and $c_i^- = c_m$. In particular, the degrees of connectivity—the number of cells connected to a given cell via RWG functions—and the graph distance—the smallest number of RWG functions needed to link a given pair of cells—depend on this choice. We provide in Fig. 2 a visual representation of two of the possible configurations along with their consequences on the connectivity of the underlying graph. The nodes of this graph, which we will refer to as the dual graph of the discretization, correspond to the cells $c_m$ of the mesh, and its edges link pairs of nodes ($c_m, c_n$) if there exists an RWG function $\varphi_i$ such that $c_i^+ = c_m$ and $c_i^- = c_n$. With this formalism, the total number of basis functions associated with a mesh containing $E$ edges is $N = \sum_{e=1}^{E} (n_e - 1)$. After replacing (4) into (1) and by Galerkin testing the equation with rotated RWG functions $\{\hat{n} \times \varphi_m\}_{m=1}^N$, we obtain the $N \times N$ linear system

$$Tj = e \quad (6)$$

where $T = jk T_s - (jk)^{-1} T_h$ with $[T_s]_{mn} = \langle \hat{n} \times \varphi_m, T_s \varphi_n \rangle$, $[T_h]_{mn} = \langle \hat{n} \times \varphi_m, T_h \varphi_n \rangle$, and $[e]_m = -\langle \hat{n} \times \varphi_m, \hat{n} \times E^{\text{inc}/\eta} \rangle$, and where $\langle f, g \rangle := \int_{\Gamma} f(r) \cdot g(r) dS(r)$.

### III. QUASI-HELMHOLTZ PROJECTORS FOR STRUCTURES WITH JUNCTIONS

The EFIE (6) is known to suffer from severe ill-conditioning and dramatic losses of accuracy at low-frequency (see [1] and references therein). These problems can be tackled starting with a Loop-Star decomposition of the solution vector, i.e., by decomposing the solution vector $j$ in (6) as $j = \Sigma s + \Lambda l + H h$, where $\Sigma$, $\Lambda$, and $H$ are the Star-to-RWG, Loop-to-RWG, and quasi-Harmonic-functions-to-RWG transformation matrices, respectively, and $s, l$, and $h$ are

<figure><img src="image_3.png" /><figcaption>Fig. 2. Two conventional ways to define an independent set of RWGs on a junction-edge and their associated graph representation. It illustrates two extreme configuration scenarios: on the left, the one giving rise to the maximal degree of connectivity, with one cell connected to all the others; on the right, the one giving rise to the minimal degree of connectivity (but with higher graph distance), with each cell connected to at most two other cells through this edge.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Dual graph (in red) associated with a mesh (in gray) which contains open boundaries, a junction, and an aperture. This graph highlights the existence of two local Loops around the vertex belonging to the junction and one global Loop around the aperture.</figcaption></figure>

the vectors containing the associated expansion coefficients. In the presence of junctions, the matrix $\Sigma$ is defined as in the junction-less case as

$$[\Sigma]_{mn} = \begin{cases} \pm 1, & \text{if } c_m^{\pm} \text{ is cell } n \\ 0, & \text{otherwise} \end{cases} \quad (7)$$

which is also referred to as the charge matrix, in the sense that it corresponds to the discretization of the divergence operator [1] and its rank equals the number of independent charges [3]. With this definition, $\Sigma$ contains the information on the connectivity between cells. This connectivity depends on the choice of the RWG configuration along the junctions (see the discussion above and Fig. 2) and impacts the spectral properties of $\Sigma$. These properties and their influence on the performance of the proposed method are further investigated in Section III-B.

Regarding the Loops, we propose the following definition based on the dual graph faces (see Fig. 3), for compactness

$$[A]_{mn} = \begin{cases} \pm 1, & \text{if } v_m^{\pm} \text{ belongs to the dual-face } n \\ 0, & \text{otherwise.} \end{cases} \quad (8)$$

The construction of the quasi-Harmonic-functions-to-RWG transformation matrix $H$ requires the identification of the global cycles associated with the genus of the geometry (see [1] and references therein).

The use of the Loop-Star basis is known to be problematic under certain circumstances because, on the one hand, it worsens the conditioning of the original EFIE [1] and, on the other hand, it requires the detection of topological Loops in $H$ that can jeopardize
the overall linear complexity of the scheme. Moreover, even the presence of open structure cycles (that in the formalism of this communication are associated with the matrix $\Lambda$) can be problematic and it can require an expensive cycle detection procedure. The reader should also note that the recovery of even standard Loops (around the vertices of the junction edges) can be problematic when several junctions are present.

The above drawbacks can be tackled by leveraging the quasi-Helmholtz projectors. Similar to what is done in the junction-less case, we define the projector to the nonsolenoidal subspace as

$$ \mathbf{P}^{\Sigma} = \mathbf{\Sigma}(\mathbf{\Sigma}^{T}\mathbf{\Sigma})^{+}\mathbf{\Sigma}^{T} \quad (9) $$

where “+” denotes the Moore-Penrose pseudoinverse [14], and the projector to the solenoidal subspace as

$$ \mathbf{P}^{\Lambda H} = \mathbf{I} - \mathbf{P}^{\Sigma}. \quad (10) $$

Finally, we use the quasi-Helmholtz decomposition provided by the projectors to compensate for the frequency scaling of the solenoidal and nonsolenoidal parts of the EFIE system when the frequency is decreasing [1], [3]. To this end, we define the preconditioning matrix $\mathbf{P} = (1/\sqrt{k})\mathbf{P}^{\Lambda H} + j\sqrt{k}\mathbf{P}^{\Sigma}$ to form the preconditioned EFIE

$$ \mathbf{P} \left( jk \mathbf{T}_s - \frac{1}{jk} \mathbf{T}_h \right) \mathbf{P} \mathbf{y} = \mathbf{P} \mathbf{e} \quad (11) $$

that exhibits a bounded condition number at low-frequency [1], and where $\mathbf{P} \mathbf{y} = j$.

In addition, to avoid loss of significant digits at low frequency, (11) needs to be expanded in order to impose $\mathbf{P}^{\Lambda H} \mathbf{T}_h = \mathbf{T}_h \mathbf{P}^{\Lambda H} = \mathbf{0}$ and the Loop contributions of the right-hand side and postprocessing field solutions need to be computed using specific techniques [1].

### A. Study of Orthogonality and Completeness

To ensure the proper applicability of the quasi-Helmholtz projectors to arbitrary problems with junctions, orthogonality between the Loop and Star bases and the completeness of the proposed quasi-Helmholtz decomposition must be verified.

The orthogonality between the Loop (or quasi-Harmonic-functions) and the Star transformation matrices follows from the fact that, across a triangle, a cycle has to enter through one edge and to exit through another one, whereas a Star associated with a given cell exits across all the branches made by the RWGs on that cell, even in the presence of junctions. It results in a null dot product between the columns of $\Lambda$ (or $\mathbf{H}$) and $\Sigma$, which finally gives $\mathbf{\Sigma}^T \Lambda = \mathbf{0}$ and $\mathbf{\Sigma}^T \mathbf{H} = \mathbf{0}$.

Next, we show the completeness of the quasi-Helmholtz decomposition, i.e., that the Loops, quasi-Harmonic functions, and Stars span the entire RWG space. We establish our result for connected geometries. The generalization to the nonconnected case is straightforward, by individually applying the following result to each connected subdomain. We first describe the elements of the dual graph: $S$ is the number of nodes that correspond to the cells of the mesh, $F$ is the number of faces that correspond to the local Loops and apertures, including the exterior boundary of open objects, and $H$ is the number of handles of the structure, present in both the mesh and its dual graph, which corresponds to $2H$ quasi-harmonic functions; $N$ denotes the number of RWGs and so the number of edges of the dual graph. We can then relate these quantities using the Euler-Poincaré formula [15] applied to the dual graph $N = (S-1) + (F-1) + 2H$. First, the number of independent cycles (Loops and quasi-Harmonic functions) is known to be $(F-1) + 2H$ [16] where the linear dependency between the nonharmonic cycles is taken into account. Next, the rank of $\Sigma$ is $S-1$ since it has a nullspace associated with

the constant vector (it is related to the fact that a div-conforming current has a neutral charge along the whole structure [3]). As a consequence, the rank of the subspace orthogonal to the range of $\Sigma$ is $N - (S-1) = F - 1 + 2H$, which corresponds to the total number of independent cycles. Because $\Lambda$ and $\mathbf{H}$ are both orthogonal to $\Sigma$, it results in the completeness between the Loop and Star bases and that $\mathbf{P}_{\Lambda H}$ is the projector on the solenoidal subspace (local Loops and global cycles).

### B. Spectral Analysis of Graph Laplacians With Junctions

The computation of the Star projector (9) requires solving systems arising from the pseudoinversion of $\mathbf{\Sigma}^T \mathbf{\Sigma}$. Thus, the spectral properties of this matrix directly impact the computational efficiency of this procedure. In particular, the definition (7) of $\Sigma$ is the one of the adjacency matrix of the dual graph described above, which in turn means that $\mathbf{\Sigma}^T \mathbf{\Sigma}$ is the graph Laplacian associated with this graph [17]. The spectral properties of this class of matrices have been thoroughly studied [17], [18], [19].

First, the graph Laplacian $\mathbf{\Sigma}^T \mathbf{\Sigma}$ has a nullspace associated with the constant functions restricted to each connected subdomain of $\Gamma$. Because the pseudoinverse of $\mathbf{\Sigma}^T \mathbf{\Sigma}$ in (9) will only be applied to vectors living in the range of $\mathbf{\Sigma}^T$, we will focus on the spectral properties of the graph Laplacian restricted to this subspace. In particular, we will consider the condition number “restricted” to that subspace that we define as the ratio between the largest singular value $\mu_{max}$ and the smallest nonzero singular value $\mu_{min}$ of $\mathbf{\Sigma}^T \mathbf{\Sigma}$, i.e., $\mathrm{cond}_r(\mathbf{\Sigma}^T \mathbf{\Sigma}) = \mu_{max}/\mu_{min}$.

In practice, the convergence of several iterative methods will be faster if this condition number is close to one, provided that the initial guess belongs to the range of $\mathbf{\Sigma}^T$ (see for instance the convergence results on the conjugate gradient method in solving singular symmetric systems [20]). The maximal singular value can be bounded as [18]

$$ \delta \le \mu_{max} \le 2\delta \quad (12) $$

where $\delta$ denotes the maximal degree of connectivity of the graph associated with $\mathbf{\Sigma}^T \mathbf{\Sigma}$. Regarding the minimal nonzero singular value, we have the upper bound [17]

$$ \mu_{\min} \le 1 - \frac{\sqrt{\delta - 1}}{\delta} \left(1 - \frac{2}{D}\right) + \frac{2}{D} \quad (13) $$

where $D$ denotes the diameter of the graph associated with $\mathbf{\Sigma}^T \mathbf{\Sigma}$. Let us first consider the case of $\delta$ growing. Combining (12) and (13) we get

$$ \mathrm{cond}_r(\mathbf{\Sigma}^T \mathbf{\Sigma}) \ge \frac{\delta}{1 - \frac{\sqrt{\delta - 1}}{\delta} \left(1 - \frac{2}{D}\right) + \frac{2}{D}}. \quad (14) $$

This inequality shows that the condition number increases asymptotically with the maximal degree of connectivity $\delta$. When $\delta$ is bounded and $D$ grows, however, (14) does not provide a useful bound, but the following [19]:

$$ 2 \left[ \sqrt{\mathrm{cond}_r(\mathbf{\Sigma}^T \mathbf{\Sigma})} \sqrt{\frac{\alpha^2 - 1}{4\alpha}} + 1 \right] \left[ \log_{\alpha} \frac{S}{2} \right] \ge D \quad (15) $$

can be used instead, where $\alpha$ is any arbitrary real number greater than 1, and $S$ is the dimension of $\mathbf{\Sigma}^T \mathbf{\Sigma}$. This shows that the condition number increases asymptotically with respect to the graph diameter $D$, provided that $\log_{\alpha}(S/2)$ increases slower than $D$. This requirement is met whenever $\delta$ is bounded, since in that case $D$ increases as $S^d$, with $d = 1$ for elongated (quasi-lineic) surfaces, $d = 1/2$ for surfaces of balanced length and width and $d = 1/3$ for
<figure><img src="image_5.png" /><figcaption>Fig. 4. Spectrum of the graph Laplacian $\Sigma^T \Sigma$ for structures with multiple sheets, each discretized with two triangles. The two connectivity configurations of Fig. 2 are tested and yield different spectra for the same structures.</figcaption></figure>

3-D lattice (quasi-volumic) structures. Equations (14) and (15) only imply an increase of the condition number in the asymptotic limit. However, these bounds are sharp enough that this increase has been observed numerically for small degrees of connectivity and small graph diameters, as shown in Fig. 4.

A direct consequence of (15), is that $cond_r(\Sigma^T \Sigma)$ increases asymptotically when refining the mesh, since $D$ will be increasing. The second consequence from the above analysis is that the condition number $cond_r(\Sigma^T \Sigma)$ also increases asymptotically when inserting additional sheets to the structure. To further clarify this, let us consider the two extreme configurations of RWG functions along the junctions, illustrated in Fig. 2, both implying different asymptotic behavior for the degrees of connectivity and graph diameters. In the first case (as in Fig. 2, left), the graph diameter is not impacted by the presence of junctions but the connectivity increases with the maximal number of sheets connected together, which results in an increase of the condition number. In the second case (as in Fig. 2, right), the connectivity is bounded but the graph diameter increases with the maximal number of sheets connected together, which also results in an increase of the condition number.

Because of this conditioning behavior, the computation of the quasi-Helmholtz projectors requires an efficient method to precondition and solve graph Laplacian systems. In Section IV, we will investigate numerically the effects of using multigrid approaches [6], [7], [8], [9] for this purpose.

## IV. NUMERICAL RESULTS

This section will first show the performance of AGMG in performing the graph Laplacian pseudoinversion required by quasi-Helmholtz projectors, then the overall preconditioning approach will be validated on simply and nonsimply connected structures containing junctions.

### A. Pseudoinversion of the Graph Laplacian

First, we investigate the use of AGMG, described in [6] and [7] and its robustness in preconditioning the graph Laplacian $\Sigma^T \Sigma$ in the presence of junctions. We study the number of iterations as a function of the number of unknowns when refining the mesh of structures having a bounded number of junctions and with the RWG configuration on junctions that maximizes the connectivity

<table><caption>TABLE I
NUMBER OF ITERATIONS ACHIEVED BY AGMG FOR INCREASING DISCRETIZATION WITH A RESIDUAL ERROR BELOW $10^{-8}$</caption><thead><tr><th>Structure</th><th>Unk.</th><th>Cond.<br>Iter.</th><th>Structure</th><th>Unk.</th><th>Cond.<br>Iter.</th></tr></thead><tbody><tr><td rowspan="6">Structure 1</td><td>7.2E3</td><td>1.1E4</td><td rowspan="6">Structure 2</td><td>2.9E3</td><td>2.5E3</td></tr><tr><td>21</td><td></td><td>18</td></tr><tr><td>6.9E4</td><td>1.1E5</td><td>2.8E5</td><td>2.3E5</td></tr><tr><td>22</td><td></td><td>24</td></tr><tr><td>7.0E5</td><td>1.0E6</td><td>2.7E6</td><td>2.3E6</td></tr><tr><td>22</td><td></td><td>22</td></tr><tr><td rowspan="6">Structure 3</td><td>1.2E3</td><td>1.5E3</td><td rowspan="6">Structure 4</td><td>3.2E3</td><td>8.4E3</td></tr><tr><td>18</td><td></td><td>19</td></tr><tr><td>1.2E5</td><td>1.3E5</td><td>3.0E5</td><td>7.5E5</td></tr><tr><td>22</td><td></td><td>22</td></tr><tr><td>1.2E6</td><td>1.3E6</td><td>3.0E6</td><td>7.5E6</td></tr><tr><td>23</td><td></td><td>23</td></tr><tr><td rowspan="6">Structure 4</td><td>9.6E3</td><td>5.1E3</td><td rowspan="6">Structure 5</td><td>4.0E3</td><td>7.5E3</td></tr><tr><td>19</td><td></td><td>19</td></tr><tr><td>2.4E5</td><td>1.2E5</td><td>4.0E5</td><td>7.1E5</td></tr><tr><td>20</td><td></td><td>22</td></tr><tr><td>2.7E6</td><td>1.2E6</td><td>1.0E6</td><td>1.8E6</td></tr><tr><td>24</td><td></td><td>24</td></tr></tbody></table>

(Fig. 2, left). Table I shows that the number of iterations of AGMG increases logarithmically with the number of unknowns for various structures while the restricted condition number increases linearly with the number of unknowns. This results in a quasi-linear number of operations since each iteration requires sparse matrix-vector products.

Next, we address the question of the AGMG performance when the number of connected sheets increases. Fig. 5 shows the time required by the setup and solving phases of AGMG as a function of the number of sheets, when the graph Laplacian is built using the two extremal RWG configurations (minimal and maximal connectivity). We observe that the setup time increases almost quadratically when the number of sheets and the connectivity both increase, while it remains linear when the number of sheets increases but the connectivity is bounded. In both cases, the solving time is quasi-linear with respect to the number of unknowns. Note that for a smaller number of junctions, we however observe better performance with the configuration having the maximal degree of connectivity. The computation time for the first experiment has not been reported above because the scenario is similar to the second configuration in Fig. 5, and numerical observations have shown similar behavior.

For this particular version of algebraic multigrid, our numerical observations thus point out that the connectivity plays an important role in the performance of the graph Laplacian inversion using AGMG and that the proper RWG configuration on junctions can be chosen to get a robust algorithm.

### B. Efficiency of the Quasi-Helmholtz Preconditioning for Structures With Junctions

We now numerically test the efficiency of the proposed method with application to the computation of the radar cross section (RCS) of a jet engine illuminated by an incident plane wave $E^{\text{inc}}(r) = \hat{z}e^{-jkr \cdot \hat{x}}$. The shell of the engine is a toroidal geometry that gives two quasi-harmonic functions. On the front, the blades of the turbine are approximated by metallic sheets and are connected to the rotor by junctions. This part of the structure is separated from the shell and generates an independent set of Loops and Stars (with respect to the Euler-Poincaré formula). In particular, $\Sigma$ has nullspace of dimension 2 generated by the constant functions restricted to each separated subdomain. On the back, a fence made of metallic sheets connects the shell to the nozzle forming junctions at both extremities of the sheets
<figure><img src="image_6.png" /><figcaption>Fig. 5. Computational time required by AGMG to solve graph-Laplacian systems for increasing number of junctions. Configuration (1) and (2) are, respectively, corresponding to the maximal and minimal degree of connectivity that we get with a given RWG configuration on junctions.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Surface density current calculated along the surface of the turbine with the quasi-Helmholtz projectors at 1 Hz.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Condition number of the EFIE with and without preconditioning in function of the frequency for the model of turbine.</figcaption></figure>

as well as apertures between them. The nozzle itself is an open surface with an aperture at its end. The object is discretized using 9191 RWGs that can be split into 3011 and 6180 Loops and Stars, respectively. The average number of iterations of AGMG required to solve a single graph Laplacian system with a desired accuracy of $10^{-6}$ is 19, which is repeated at every multiplication involving the projectors in the iterative inversion of the preconditioned EFIE. The surface current density calculated at a frequency of 1 Hz using the EFIE preconditioned with the quasi-Helmholtz projectors is given in Fig. 6.

Our work (11) is compared with the standard EFIE (without preconditioning) as well as the Loop-Star decomposed EFIE [1]. In Fig. 7, we observe that the two preconditioning techniques give rise to systems with a bounded condition number when the frequency

<figure><img src="image_9.png" /><figcaption>Fig. 8. Bistatic RCS (dBsm) in the (x̃, z̃) plane given the scattered far-field produced by the turbine with an incident plane wave at a frequency of 1 Hz.</figcaption></figure>

<table><caption>TABLE II
NUMBER OF ITERATIONS AND RUN TIME OF TFQMR TO SOLVE THE EFIE WITH AND WITHOUT PRECONDITIONING AT MEDIUM AND LOW FREQUENCIES FOR THE MODEL OF TURBINE WITH 9191 RWGS</caption><thead><tr><th rowspan="2">Freq.</th><th colspan="2">Standard EFIE</th><th colspan="2">Loop-Star</th><th colspan="2">This work</th></tr><tr><th>Iter.</th><th>Time (s)</th><th>Iter.</th><th>Time (s)</th><th>Iter.</th><th>Time (s)</th></tr></thead><tbody><tr><td>5 MHz</td><td>1632</td><td>176.7</td><td>4064</td><td>690.8</td><td>393</td><td>137.2</td></tr><tr><td>1 Hz</td><td colspan="2">No convergence</td><td>4416</td><td>779.0</td><td>419</td><td>142.3</td></tr></tbody></table>

decreases, while the condition number dramatically increases for the standard formulation until it saturates due to numerical effects. In addition, the condition numbers obtained with the quasi-Helmholtz projectors are significantly lower than the Loop-Star decomposed EFIE. Solving with the solver TFQMR at 1 Hz for a desired accuracy of $10^{-6}$, the iterative procedure requires fewer iterations when the system is preconditioned with the quasi-Helmholtz projectors than when preconditioned using a Loop-Star decomposition since only 419 iterations are required for the former instead of 4416 for the latter. Under these conditions, the EFIE does not converge without preconditioning. We also note that the two preconditioning methods give comparable results at low frequency, unlike the standard formulation (solved via direct factorization), as is shown in Fig. 8. Finally, in Table II we report the number of iterations and computation time required to solve the EFIE with and without preconditioning, showing better performance for the proposed scheme at a frequency of 5 MHz.

V. CONCLUSION

In this communication, we theoretically and numerically showed how the quasi-Helmholtz projectors method can be used to model low-frequency problems containing junctions with the EFIE. Expla- nations concerning the generalization of the Loops and Stars for such structures have been given. Finally, we have presented the successful application of the quasi-Helmholtz projectors to scenarios of industrial relevance.

REFERENCES

[1] S. B. Adrian, A. Dély, D. Consoli, A. Merlini, and F. P. Andriulli, "Electromagnetic integral equations: Insights in conditioning and preconditioning," *IEEE Open J. Antennas Propag.*, vol. 2, pp. 1143-1174, 2021.

[2] G. Vecchi, "Loop-star decomposition of basis functions in the discretization of the EFIE," *IEEE Trans. Antennas Propag.*, vol. 47, no. 2, pp. 339-346, Feb. 1999.
[3] J.-F. Lee, R. Lee, and R. J. Burkholder, "Loop star basis functions and a robust preconditioner for EFIE scattering problems," *IEEE Trans. Antennas Propag.*, vol. 51, no. 8, pp. 1855-1863, Aug. 2003.

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

[5] D. R. Wilton, "Topological consideration in surface patch and volume cell modeling of electromagnetic scatterers," in *Proc. URSI Int. Symp. Electromagn. Theory*, 1983, pp. 65-68.

[6] A. Napov and Y. Notay, "An algebraic multigrid method with guaranteed convergence rate," *SIAM J. Sci. Comput.*, vol. 34, no. 2, pp. A1079-A1109, Jan. 2012.

[7] Y. Notay, "An aggregation-based algebraic multigrid method," *Electron. Trans. Numer. Anal.*, vol. 37, pp. 123-146, Jan. 2010.

[8] M. Bolten, S. Friedhoff, F. Andreas, H. Matthis, and K. Karsten, "Algebraic multigrid methods for Laplacians of graphs," *Linear Algebra Its Appl. Linear Algebra Appl.*, vol. 434, no. 11, pp. 2225-2243, Jun. 2009.

[9] C. Iwamura, F. S. Costa, I. Sbarski, A. Easton, and N. Li, "An efficient algebraic multigrid preconditioned conjugate gradient solver," *Comput. Methods Appl. Mech. Eng.*, vol. 192, nos. 20-21, pp. 2299-2318, May 2003.

[10] J. Bourhis, A. Merlini, and F. P. Andriulli, "On the effectiveness of quasi-helmholtz projectors in preconditioning problems with junctions," in *Proc. IEEE Int. Symp. Antennas Propag. USNC-URSI Radio Sci. Meeting (AP-S/URSI)*, Jul. 2022, pp. 1944-1945.

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

[12] M. Carr, E. Topsakal, and J. L. Volakis, "A procedure for modeling material junctions in 3-D surface integral equation approaches," *IEEE Trans. Antennas Propag.*, vol. 52, no. 5, pp. 1374-1378, May 2004.

[13] P. Ylä-Oijala, M. Taskinen, and J. Sarvas, "Surface integral equation method for general composite metallic and dielectric structures with junctions," *Prog. Electromagn. Res.*, vol. 52, pp. 81-108, 2005.

[14] A. Ben-Israel and T. N. Greville, *Generalized Inverses: Theory and Applications*, vol. 15. New York, NY, USA: Springer, 2003.

[15] C. Dodson and P. Parker, *A User's Guide to Algebraic Topology*. New York, NY, USA: Springer, 1997.

[16] J. L. Gross and J. Yellen, *Graph Theory and Its Applications*. Boca Raton, FL, USA: CRC Press, 2005.

[17] N. M. M. de Abreu, "Old and new results on algebraic connectivity of graphs," *Linear Algebra Its Appl.*, vol. 423, no. 1, pp. 53-73, 2007.

[18] R. Merris, "Laplacian matrices of graphs: A survey," *Linear Algebra Its Appl.*, vol. 197, pp. 143-176, Jan. 1994.

[19] B. Mohar, "Eigenvalues, diameter, and mean distance in graphs," *Graphs Combinatorics*, vol. 7, no. 1, pp. 53-64, Mar. 1991.

[20] W. J. Kammerer and M. Z. Nashed, "On the convergence of the conjugate gradient method for singular linear operator equations," *SIAM J. Numer. Anal.*, vol. 9, no. 1, pp. 165-181, Mar. 1972.