# Calderón Strategies for the Convolution Quadrature Time Domain Electric Field Integral Equation

Pierrick Cordel, Alexandre Dély, Adrien Merlini, Francesco Andriulli

► To cite this version:

Pierrick Cordel, Alexandre Dély, Adrien Merlini, Francesco Andriulli. Calderón Strategies for the Convolution Quadrature Time Domain Electric Field Integral Equation. IEEE Open Journal of Antennas and Propagation, 2024, 1, pp.1-1. 10.1109/OJAP.2024.3354044. [hal-04405445](https://hal.science/hal-04405445v1)

HAL Id: hal-04405445

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

Submitted on 6 Feb 2024

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
<figure><img src="image_2.png" /></figure>

Received 21 October, 2023; revised 07 December, 2023; accepted 11 December, 2023; Date of publication XX Month, XXXX; date of current version 11 December, 2023.

Digital Object Identifier 10.1109/XXXXXXXXXXXX

# Calderón Strategies for the
Convolution Quadrature Time Domain
Electric Field Integral Equation

PIERRICK CORDEL*, Student Member, IEEE, ALEXANDRE DÉLY*, ADRIEN MERLINI†,
Senior Member, IEEE AND FRANCESCO P. ANDRIULLI*, Fellow, IEEE

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

This work was supported by the European Research Council (ERC) through the European Union's Horizon 2020 Research and Innovation Programme under Grant 724846 (Project 321) and by the H2020-MSCA-ITN-EID project COMPETE GA No 955476.

**ABSTRACT** In this work, we introduce new integral formulations based on the convolution quadrature method for the time-domain modeling of perfectly electrically conducting scatterers that overcome some of the most critical issues of the standard schemes based on the electric field integral equation (EFIE). The standard time-domain EFIE-based approaches typically yield matrices that become increasingly ill-conditioned as the time-step or the mesh discretization density increase and suffer from the well-known DC instability. This work presents solutions to these issues that are based both on new Calderón strategies and quasi-Helmholtz projectors regularizations. In addition, to ensure an efficient computation of the marching-on-in-time, the proposed schemes leverage properties of the Z-transform—involved in the convolution quadrature discretization scheme—when computing the stabilized operators. The two resulting formulations compare favorably with standard, well-established schemes. The properties and practical relevance of these new formulations will be showcased through relevant numerical examples that include canonical geometries and more complex structures.

**INDEX TERMS** Boundary element method, Calderón preconditioning, computational electromagnetic, convolution quadrature method, EFIE, time-domain integral equations

## I. INTRODUCTION

TIME domain boundary integral equations (TDIEs) are widely used in the simulation of transient electromagnetic fields scattered by perfectly electrically conducting (PEC) objects [1]–[4]. Like their frequency-domain counterparts, the spatial discretization of these equations is often performed via the boundary element method. The time discretization, however, can be tackled in different ways. A popular approach leverages time basis functions either within a Marching-On-in-Time (MOT) scheme [5]–[7] or within a Marching-On-In-Order procedure [8]. The convolution quadrature (CQ) approach [9], [10] is an attractive alternative to these methods in which only space basis functions are explicitly defined. The approach has been applied to several equations in elastodynamics and acoustics [11], [12] and then in electromagnetics [13]. It provides an efficient time-

stepping scheme with matrices derived from the space-discretized Laplace domain operators.

Another advantage of the CQ method is the use of implicit schemes (e.g. Runge Kutta methods [14]–[18]), which are generally more stable and typically allow for a better accuracy control of the solution over time [19], [20]. However, the CQ time stepping scheme is solved via a computationally expensive MOT algorithm. Nowadays, fast solvers can reach quasi-linear complexity in time and space [21], [22]. Usually, this fast technology uses iterative solvers, resulting in an overall computational cost that is proportional to the number of iterations which is low for well-conditioned systems. Working with well-conditioned matrices is therefore essential to reduce the computational cost of the solution process, in addition to being necessary to obtain accurate results [23].

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

<footnote n="2">²Microwave Department, IMT Atlantique, 29285 Brest, France</footnote>
Lamentably, however, the CQ discretized time domain electric field integral equation (EFIE) is plagued by several drawbacks. Indeed, the matrices resulting from the discretization of the EFIE are known to become ill-conditioned for large time steps or at dense mesh discretizations: the condition number of the MOT matrices grows quadratically with the time step and with the inverse of the average mesh edge length. These two phenomena are the CQ counterparts of what for standard MOT schemes are known as the large time step breakdown [24]–[27] and the dense discretization breakdown (or h-refinement breakdown) [28], [29]. Another challenge in handling the CQ EFIE is that it involves operators whose definitions include a time integration. To avoid dealing with this integral, the time-differentiated counterpart of this formulation is often used [13], [30], but this differentiation is subject to a source of instability in the form of spurious linear currents living in the nullspace of the operator that degrades the solution [28], [31]; this phenomenon is known as the direct current instability (DC instability).

In this work, we propose new Calderón-preconditioned and quasi-Helmholtz regularized formulations free from the limitations mentioned above. The Calderón identities they rely on are already a well-established preconditioning approach in both the frequency domain [32] and time domain discretized by the Galerkin method [28], [29], [33], [34] that is extended in this work to convolution quadrature discretizations and complemented with quasi-Helmholtz regularization. The contribution of this paper is twofold: (i) we present a first approach to tackle the regularization of the EFIE operator and to address the DC instability resulting in a new operator that presents no nullspace on simply connected geometries, thus stabilizing the solution, and (ii) we build upon this first regularized form of the EFIE to obtain an equation that, at the price of a higher number of matrix-vector products, is stable in the case of multiply connected geometries.

This article is structured as follows: the time domain formulations of interest are summarized in Section II along with the convolution quadrature method and the boundary element method for spatial discretization; in Section III, the new Calderón and projectors-based preconditioning strategies are presented; finally, Section IV presents the numerical studies that confirm the effectiveness of the different approaches before concluding. Preliminary studies pertaining to this work were presented in the conference contribution [35].

## II. Background and Notations

### A. Time Domain Integral Formulations

In this work, we consider the problem of time-domain scattering by a perfectly electrically conducting object that resides in free space. The object is illuminated by an electromagnetic field ($e^{\text{inc}}$, $h^{\text{inc}}$)($\mathbf{r}, t$) which induces a surface current density $j_{\Gamma}$ on its boundary $\Gamma$ that is the solution of the time-domain EFIE

$$ \eta_0 \mathcal{T}(j_{\Gamma}) (\mathbf{r}, t) = -\hat{\mathbf{n}}(\mathbf{r}) \times \mathbf{e}^{\text{inc}}(\mathbf{r}, t). \quad (1) $$

Here, $\hat{\mathbf{n}}$ is the outpointing normal to $\Gamma$ and $\eta_0$ is the characteristic impedance of the background. The electric field operator $\mathcal{T}$ includes the contributions of the vector and scalar potentials, respectively denoted $\mathcal{T}_s$ and $\mathcal{T}_h$ [30]

$$ \mathcal{T}(\mathbf{f}) (\mathbf{r}, t) = -\frac{1}{c_0} \frac{\partial}{\partial t} \mathcal{T}_s(\mathbf{f}) (\mathbf{r}, t) + c_0 \int_{-\infty}^{t} \mathcal{T}_h(\mathbf{f}) (\mathbf{r}, t') dt', \quad (2) $$

$$ \mathcal{T}_s(f)(\mathbf{r}, t) = \hat{\mathbf{n}}(\mathbf{r}) \times \iint_{\mathbf{r}' \in \Gamma} (\mathcal{G}_r *_{t} f)(\mathbf{r}', t) dS', \quad (3) $$

$$ \mathcal{T}_h(f)(\mathbf{r}, t) = \hat{\mathbf{n}}(\mathbf{r}) \times \nabla \iint_{\mathbf{r}' \in \Gamma} (\mathcal{G}_r *_{t} \nabla' \cdot f)(\mathbf{r}', t) dS', \quad (4) $$

where $c_0$ is the speed of light in the background. The temporal convolution product $*_t$ and the temporal Green function $\mathcal{G}$ are defined as

$$ (f *_{t} g)(t) = \int_{-\infty}^{\infty} f(\tau) g(t - \tau) d\tau, \quad (5) $$

$$ \mathcal{G}_r(\mathbf{r}', t) = \frac{\delta\left(t - \frac{|\mathbf{r} - \mathbf{r}'|}{c_0}\right)}{4\pi|\mathbf{r} - \mathbf{r}'|}, \quad (6) $$

with $\delta$ the time Dirac delta.

### B. Marching-On-In-Time with Convolution Quadratures

Let $\Theta$ be a placeholder for any of the integral operators previously presented and let $k(\mathbf{r}, t)$ be a causal function ($\forall t < 0, k(\mathbf{r}, t) = 0$). With these notations, most time domain integral equation take the form

$$ \Theta(f_c)(\mathbf{r}, t) = k(\mathbf{r}, t), \quad (7) $$

where $f_c$ is the solution to be solved for. The first step of the Marching-On-in-Time solution scheme with convolution quadratures is to apply the boundary element method [3], [36], [37] as spatial discretization. Assuming separability between the space and time variables, the unknown function $f_c$ is expanded as a linear combination of $N_s$ spatial basis functions such that [38]

$$ f_c(\mathbf{r}, t) \approx \sum_{n=1}^{N_s} [f_{\Gamma}]_n(t) f_n^{\text{src}}(\mathbf{r}), \quad (8) $$

where $\{f^{\text{src}}\}$ are the source spatial basis functions and their associated time coefficients are stored in the vector $f_{\Gamma}(t)$. Then, the equation (7) is tested by the spatial basis functions $\{f^{\text{tst}}\}$ leading to the time-dependent matrix system

$$ (\boldsymbol{\theta} * \mathbf{f})(t) = \mathbf{k}(t), \quad (9) $$

where for $n$ and $m$ in $[1, N_s]$, we have,

$$ (\boldsymbol{\theta} * \mathbf{f})_m(t) = \sum_{n=1}^{N_s} (\boldsymbol{\theta}_{m,n} *_{t} f_n)(t), \quad (10) $$

$$ [\boldsymbol{\theta}]_{m,n}(t) = \langle f_m^{\text{tst}}, \boldsymbol{\Theta} (\delta f_n^{\text{src}}) \rangle_{\Gamma}(t), $$

$$ [\mathbf{k}_{\Gamma}]_m(t) = \langle f_m^{\text{tst}}, \mathbf{k} \rangle_{\Gamma}(t), $$

with

$$ \langle \mathbf{f}, \mathbf{g} \rangle_{\Gamma} = \iint_{\mathbf{r} \in \Gamma} \mathbf{f}(\mathbf{r}) \cdot \mathbf{g}(\mathbf{r}) \, dS. \quad (11) $$
<figure><img src="image_3.png" /></figure>

The second step is the discretization in time with the con-
volution quadrature method [9], [10], [12]–[14], [38]. First,
the system (9) must be transformed in the Laplace domain
[39] and we denote θ_L, f_L and k_L the Laplace transform of
θ, f_R and k_R. The system (9) is then equivalent to

$$
\mathcal{L}\{t \mapsto (\boldsymbol{\theta} * \mathbf{f})(t)\}(s) = \boldsymbol{\theta}_{\mathcal{L}}(s) \mathbf{f}_{\mathcal{L}}(s) = \mathbf{k}_{\mathcal{L}}(s), \quad (12)
$$

where *s*, a complex number, is the Laplace parameter. Then,
a representation on the Z-domain discretizes the system (12).
The Laplace parameter *s*, in the operator *θ*<sub>*L*</sub>, is replaced by
the matrix-valued parameter

$$
\mathbf{s}_{\text{cq}}(z) = \frac{1}{\Delta t} \left( \mathbf{A} + \frac{\mathbf{1}_p}{z-1} \mathbf{b}^T \right)^{-1}, \quad (13)
$$

diagonalizable for the considered z values, with the following
eigenvalue decomposition s<sub>cq</sub>(z) = Q(z)Λ(z)Q<sup>-1</sup>(z). The
elements of the diagonal matrix Λ(z) are the eigenvalues
of s<sub>cq</sub>(z) and the columns of Q(z) are their associated
eigenvectors. The time step size is denoted Δt and the matrix
A and the vectors 1<sub>p</sub>, c, b of size p are determined by
the implicit scheme used [40]. The discretized Z-domain
operator θ<sub>Z</sub> is defined such that for k and l in [1, p] and
for m and n in [1, N<sub>s</sub>]

$$
[\theta_Z]_{\Phi_{m,k}, \Phi_{n,l}}(z) = [Q(z)\theta_{m,n}^{\Lambda}(z)Q^{-1}(z)]_{k,l}, \quad (14)
$$

where $\Phi_{\alpha,\beta} = p(\alpha-1)+\beta$ is an appropriate indexing function and $\theta_{m,n}^{\Lambda}(z)$ is a diagonal matrix defined as

$$
\left[ \theta_{m,n}^{\Lambda} \right]_{k,k}(z) = \left[ \theta_{\mathcal{L}}(\Lambda_{k,k}) \right]_{m,n}(z) . \quad (15)
$$

The vectors $\mathbf{f}_Z$ and $\mathbf{k}_Z$ are the Z-domain representation [40] of the respective time-discretized vectors

$$
\begin{align}
[k_i]_{\Phi_{m,k}} &= [k_{\Gamma}]_m (\Delta t (i + [c]_k)), \nonumber \\
[f_i]_{\Phi_{n,k}} &= [f_{\Gamma}]_n (\Delta t (i + [c]_k)), \tag{16}
\end{align}
$$

yielding to the following discretization of the system (12)

$$
\theta_Z(z) \mathbf{f}_Z(z) = \mathbf{k}_Z(z). \quad (17)
$$

Finally, the equivalent time discretized system of (9) is
obtained by applying the inverse Z-transform on (17)

$$
\mathcal{Z}^{-1}\{\theta_Z(z)\mathbf{f}_Z(z)\}_i = [\mathbf{Z}_{\theta}*_s \mathbf{f}]_i = \sum_{j=0}^{i} \mathbf{Z}_{\theta,j} \mathbf{f}_{i-j} = \mathbf{k}_i, \quad (18)
$$

where $Z_{\theta,i} = Z^{-1}\{z \to \theta_Z(z)\}_i$ are the time domain inter-
action matrices and $*_s$ is the sequence convolution product.
The system sequence (18) is rewritten in the following
Marching-On-In-Time that can be solved for $\mathbf{f}_i$

$$
Z_{\theta,0} \mathbf{f}_i = \mathbf{k}_i - \sum_{j=1}^{i} Z_{\theta,j} \mathbf{f}_{i-j}. \quad (19)
$$

C. Classic Integral Marching-On-In-Times

In this subsection, the discretization scheme described above
will be applied to the specific case of the EFIE. The Rao-
Wilton-Glisson (RWG) basis functions {$f_n^{\text{rwg}}$}$_{N_s}$ [37], [41]
are used to expand the current density as

$$
j_{\Gamma}(\mathbf{r}, t) \approx \sum_{n=1}^{N_s} [j_{\Gamma}]_n (t) f_n^{\text{rwg}}(\mathbf{r}), \quad (20)
$$

where the current coefficients are gathered in an unknown
vector function of time j<sub>Γ</sub>(t). The EFIE is then tested with
rotated RWG basis functions {$\hat{n} \times f_n^{\text{rwg}}$}$_{N_s}$, leading to the
following Marching-On-In-Time

$$
Z_{\tau,0} j_i = -\eta_0^{-1} e_i^{\text{inc}} - \sum_{j=1}^{i} Z_{\tau,j} j_{i-j}, \quad (21)
$$

where the vector sequences $j_i$ and $e_i^{\text{inc}}$, and the time domain interaction matrices $Z_{\tau,i}$ are respectively generated by the convolution quadrature method described in Subsection B of $j_{\Gamma}(t)$ and the following space-discretized vector and matrix

$$
\left[ \mathbf{e}_{\Gamma}^{\text{inc}} \right]_{m}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_{m}^{\text{rwg}}, \hat{\mathbf{n}} \times \mathbf{e}^{\text{inc}} \rangle_{\Gamma}(t), \quad (22)
$$

$$
[\mathcal{T}]_{m,n}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{rwg}}, \mathcal{T}(\delta \mathbf{f}_n^{\text{rwg}}) \rangle_{\Gamma}(t) . \quad (23)
$$

However, the time integral contribution of this operator $\mathcal{T}$ involves an unbounded number of non-vanishing matrices $Z_{\tau,i}$ (21), leading to a prohibitive quadratic complexity with the number of time steps [42]. Historically, the time differentiated formulation is preferred because it is not afflicted by this drawback [13], [30], and leads to the following MOT [30]

$$
Z_{\tau,0} j_i = -\eta_0^{-1} \dot{e}_i^{\text{inc}} - \sum_{j=1}^{N_{\text{conv}}} Z_{\tau,j} j_{i-j}, \quad (24)
$$

where $\dot{e}_n^{\text{inc}}$ and $\dot{Z}_n$ are respectively the time domain vec-
tors and interaction matrices generated by the convolution
quadrature method described in Subsection B of

$$
\left[ \dot{\mathbf{e}}_{\Gamma}^{\text{inc}} \right]_{m}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_{m}^{\text{rwg}}, \hat{\mathbf{n}} \times \frac{\partial}{\partial t} \mathbf{e}^{\text{inc}} \rangle_{\Gamma}, \quad (25)
$$

$$
\left[ \dot{\mathcal{T}} \right]_{m,n}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{rwg}}, \frac{\partial}{\partial t} \mathcal{T} \left( \delta \mathbf{f}_n^{\text{rwg}} \right) \rangle_{\Gamma}. \quad (26)
$$

D. EFIE DC instability

The electric field integral operator suffers from the DC instability: since for all constant-in-time solenoidal current $j_{cs}$ we have $\nabla \cdot j_{cs} = 0$ and $\frac{\partial}{\partial t} j_{cs} = 0$, we can conclude that

$$
\mathcal{T}(j_{cs}) = \mathbf{0}. \tag{27}
$$

Therefore, the EFIE solution is only determined up to
a constant solenoidal current [30]. Its time differentiated
counterpart inherits these drawbacks and amplifies the DC
instability by further adding linear in time solenoidal cur-
rents to the nullspace. This latter deteriorates the late time
simulation in which spurious currents grow exponentially in
the operator nullspace [43]. This behaviour is predicted by
the polynomial eigenvalues analysis of the MOT: a stable
MOT has all its eigenvalues inside the unit circle in the complex plane while a MOT that suffers from the DC instability has some eigenvalues that cluster around 1 [44]. The eigenvalue distribution of the time differentiated EFIE MOT is represented in Figure 1 in which such a cluster is clearly visible around 1.

<figure><img src="image_4.png" /><figcaption>FIGURE 1. Polynomial eigenvalues of the TD-EFIE MOT scheme on the unit sphere, Ns = 270 and Δt = 3 ns</figcaption></figure>

E. Quasi-Helmholtz projectors

Previous works show that the electric field integral equation discretized in space using RWG basis functions can be sta-
bilized by the quasi-Helmholtz projectors [30], [34]. These
projectors are formed from the star-to-rwg transformation
matrix, denoted Σ and defined in [45], which maps the
discretized current into the non-solenoidal contributions [46],
[47]. The quasi-Helmholtz projectors on the non-solenoidal
space and its complementary (the one on solenoidal/quasi-
harmonic space) are respectively

$$
P^{\Sigma} = \left( \Sigma \left( \Sigma^T \Sigma \right)^{+} \Sigma^T \right) \quad \text{and} \quad P^{\Lambda H} = I - P^{\Sigma}, \quad (28)
$$

where + denotes the Moore–Penrose pseudoinverse [45].

III. Calderón preconditioning of the EFIE

EFIE formulations based on the quasi-Helmholtz projectors
cure the DC instability and the conditioning at large time
steps. However, these formulations still suffer from a dense
discretization breakdown. One appealing strategy could be to
apply standard preconditioning schemes to the Marching-On-
In-Time matrices directly to cure the matrix conditioning is-
sues. However, the solution currents would remain unaltered
and subject to DC instabilities as the original scheme. This is
why the preconditioning has to be performed on the continu-
ous equations to build a new operator without nullspace and
then discretize the formulation to obtain a well-conditioned
scheme. In this part, Calderón preconditioning strategies are
proposed to cure the DC instability and the conditioning
breakdowns.

**A. A Convolution Quadrature Calderón time-domain EFIE**

Calderón preconditioners are based on the Calderón identity
[29], [48]

$$
\mathcal{T}^2 = \mathcal{T} \circ \mathcal{T} = -\mathcal{I}/4 + \mathcal{K}^2, \quad (29)
$$

where ◦ is the composition operator, I is the identity
operator and K is defined as

$$
\mathcal{K}(f)(\mathbf{r}, t) = \hat{\mathbf{n}}(\mathbf{r}) \times \nabla \times \iint_{\mathbf{r}' \in \Gamma} (\mathcal{G}_{\mathbf{r}} *_{\mathbf{t}} f)(\mathbf{r}', t) dS'. \quad (30)
$$

The operator $-I/4 + \mathcal{K}^2$ is a well-behaved operator for increasing discretization densities. As a consequence, with a proper discretization, $\mathcal{T}^2$ is well-conditioned for large time steps and dense meshes for simply connected structures [45]. In practice, a discretization of $\mathcal{T}^2$ is used in which the right EFIE operator is discretized with RWG basis functions and the left preconditioner is discretized with Buffa-Christiansen (BC) basis functions ($f_n^{bc}$) [49]–[51]

$$
[\mathbb{T}]_{m,n}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{bc}}, \mathcal{T}(\delta \mathbf{f}_n^{\text{bc}}) \rangle_{\Gamma}(t). \quad (31)
$$

The preconditioning leads to the following space-discretized
formulation

$$
\left( \mathbb{T} \mathbf{G}_m^{-1} * \left( \mathbb{T} * j_{\Gamma} \right) \right) (t) = - \eta_0^{-1} \left( \mathbb{T} \mathbf{G}_m^{-1} * e_{\Gamma}^{\text{inc}} \right) (t), \quad (32)
$$

where the matrix $\mathbf{G}_m$ is the mixed gram matrix linking the two discretizations

$$
[\mathbf{G}_m]_{m,n} = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{rwg}}, \mathbf{f}_n^{\text{bc}} \rangle_{\Gamma}. \quad (33)
$$

Then, the convolution quadrature leads to the MOT scheme

$$
\begin{equation}
\begin{split}
\left[ Z_{\mathbb{T}} \tilde{G}_m^{-1} *_{s} Z_{\mathbb{T}} \right]_{0} j_i = & -\eta_0^{-1} \left[ Z_{\mathbb{T}} \tilde{G}_m^{-1} *_{s} e^{\text{inc}} \right]_{i} \\
& - \sum_{j=1}^{i} \left[ Z_{\mathbb{T}} \tilde{G}_m^{-1} *_{s} Z_{\mathbb{T}} \right]_{j} j_{i-j},
\end{split}
\tag{34}
\end{equation}
$$

where $Z_{\mathbb{T},i}$ are the time domain interaction matrices of the space-discretized operator $\mathbb{T}(t)$ (31) generated by the convolution quadrature method described in Subsection B, the sequence convolution quadrature product $*_s$ is the discretization of the space-discretized temporal convolution product * and

$$
\tilde{G}_m = G_m \otimes I_p . \qquad (35)
$$

The Kronecker product ⊗Iₚ is required to match with the
convolution quadrature method where Iₚ is the identity
matrix of size p. Unfortunately, the MOT in (34), involves
operators with temporal integrations leading to a time con-
suming MOT. A more favorable scheme can be obtained by
noticing the following commutative properties

$$
\frac{\partial}{\partial t} \mathcal{T}_{s/h}(f)(\mathbf{r}, t) = \mathcal{T}_{s/h}\left(\frac{\partial}{\partial t} f\right)(\mathbf{r}, t), \quad (36)
$$

$$
\int_{-\infty}^{t} \mathcal{T}_{s/h}(f)(\mathbf{r}, t') dt' = \mathcal{T}_{s/h}\left(\int_{-\infty}^{t} f(\mathbf{r}, t') dt'\right), \quad (37)
$$
<figure><img src="image_5.png" /></figure>

and the cancellation property $\mathcal{T}_h^2 = \mathbf{0}$ [32], we have

$$ \begin{aligned} \mathcal{T}^2 &= \left(-\frac{1}{c_0} \frac{\partial}{\partial t} \mathcal{T}_s + c_0 \int_{-\infty} \dot{\mathcal{T}}_h \right) \circ \left(-\frac{1}{c_0} \frac{\partial}{\partial t} \mathcal{T}_s + c_0 \int_{-\infty} \dot{\mathcal{T}}_h \right) \quad (38) \\ &= c_0^{-2} \frac{\partial^2}{\partial t^2} \mathcal{T}_s^2 - \mathcal{T}_s \mathcal{T}_h - \mathcal{T}_h \mathcal{T}_s. \end{aligned} $$

This is advantageous since, besides not involving any time integration contribution, the operator $c_0^{-2} \frac{\partial^2}{\partial t^2} \mathcal{T}_s^2 - \mathcal{T}_s \mathcal{T}_h - \mathcal{T}_h \mathcal{T}_s$ has no nullspace for simply connected geometries leading to a DC-stable discretization (“dottrick TDEFIE”) [52]. By extending the previous notations on $\mathcal{T}_s$ and $\mathcal{T}_h$

$$ [\mathcal{T}_{s/h}]_{m,n}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{rwg}}, \mathcal{T}_{s/h} (\delta \mathbf{f}_n^{\text{rwg}}) \rangle_{\Gamma}(t), \quad (39) $$

$$ [\mathcal{T}_{s/h}]_{m,n}(t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{bc}}, \mathcal{T}_{s/h} (\delta \mathbf{f}_n^{\text{bc}}) \rangle_{\Gamma}(t), \quad (40) $$

the proposed space-discretized operator is denoted

$$ T_c = c_0^{-2} \frac{\partial^2}{\partial t^2} \mathcal{T}_s \mathbf{G}_m^{-1} * \mathcal{T}_s - \mathcal{T}_s \mathbf{G}_m^{-1} * \mathcal{T}_h - \mathcal{T}_h \mathbf{G}_m^{-1} * \mathcal{T}_s, \quad (40) $$

yielding to the following space-discretized formulation

$$ (T_c * j\Gamma)(t) = -\eta_0^{-1} (\mathcal{T}\mathbf{G}_m^{-1} * \mathbf{e}_\Gamma^{\text{inc}})(t). \quad (41) $$

The right-hand side operator still involves a temporal integration in (41). However, given the commutative properties (37), the temporal integral on the scalar potential $\mathcal{T}_h$ is evaluated with the incident field

$$ c_0 \left( \int_{-\infty} \dot{\mathcal{T}}_h \mathbf{G}_m^{-1} * \mathbf{e}_\Gamma^{\text{inc}} \right) (t) = c_0 \left( \mathcal{T}_h \mathbf{G}_m^{-1} * \mathbf{e}_\Gamma^{\text{prim}} \right) (t), \quad (42) $$

where

$$ \left[ \mathbf{e}_{\Gamma}^{\text{prim}} \right]_m (t) = \langle \hat{\mathbf{n}} \times \mathbf{f}_m^{\text{rwg}}, \int_{-\infty}^t \hat{\mathbf{n}} \times \mathbf{e}^{\text{inc}}(t') dt' \rangle_{\Gamma}, \quad (43) $$

$$ \left[ \mathbf{e}_i^{\text{prim}} \right]_{\Phi_{m,k}} = \left[ \mathbf{e}_{\Gamma}^{\text{prim}} \right]_m \left( \Delta t (i + [c]_k) \right). \quad (44) $$

Therefore, the previous MOT is rewritten as

$$ Z_{T_c, 0j} = -\eta_0^{-1} \sum_{j=0}^{N_{\text{conv}}} \left( Z_{T_s, j} \mathbf{e}_{i-j}^{\text{inc}} + Z_{T_h, j} \mathbf{e}_{i-j}^{\text{prim}} \right) - \sum_{j=1}^{N_{\text{conv}}} Z_{T_c, j} \dot{j}_{i-j}, \quad (44) $$

where the time domain interaction matrices $Z_{T_c,i}$, $Z_{T_s,i}$ and $Z_{T_h,i}$ are respectively generated by the convolution quadrature method described in Subsection B of the space-discretized operators $T_c$, $\dot{T}_s = -c_0^{-1} \frac{\partial}{\partial t} \mathcal{T}_s \mathbf{G}_m^{-1}$ and $T_h = c_0 \mathcal{T}_h \mathbf{G}_m^{-1}$. As in (34), the interaction matrix sequence $Z_{T_c,i}$ involves computationally expensive sequence convolution products $*_s$, however, the convolution quadrature method allows the substitution of the sequence convolution products $*_s$ by matrix multiplications in the Z-domain, that can be evaluated at a lesser cost. By extending the notations of the convolution quadrature described in Subsection B on the space-discretized operators $T_{s/h}$ and $T_{s/h}$ and by using the Z-domain properties, the matrix sequence $Z_{T_c,i}$ is equal to

$$ Z_{T_c,i} = c_0^{-2} \mathcal{Z}^{-1} \left\{ \bar{s}_{\text{cq}}^2 \mathcal{T}_s, \mathcal{Z} \tilde{G}_m \mathcal{T}_s, \mathcal{Z} \right\} \\ - \mathcal{Z}^{-1} \left\{ \mathcal{T}_h, \mathcal{Z} \tilde{G}_m \mathcal{T}_s, \mathcal{Z} + \mathcal{T}_s, \mathcal{Z} \tilde{G}_m \mathcal{T}_h, \mathcal{Z} \right\}, \quad (45) $$

where the matrix $\bar{s}_{\text{cq}}(z) = I_{N_s} \otimes s_{\text{cq}}(z)$ is the Z-discretization of the time derivative and $I_{N_s}$ is the identity matrix of size $N_s$. The formulation (44) is a good candidate to obtain a stable current solution, however, the proposed operator $\mathcal{T}^2$ has static nullspaces for multiply connected geometries [33]. As such, (40) is still subject to DC-instabilities for multiply connected geometries. The polynomial eigenvalue analysis on a sphere and on a torus (respectively Figure 2 and Figure 3) illustrate this phenomenon. While all the eigenvalues cluster in 0 in the spherical case, an analysis on a torus highlights four eigenvalues of this MOT clustered around 1, corresponding to the four constant regime solutions [33].

<figure><img src="image_6.png" /><figcaption>FIGURE 2. Polynomial eigenvalues of the Calderón EFIE MOT scheme on the unit sphere, $N_s = 270$ and $\Delta t = 3$ ns.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 3. Polynomial eigenvalues of the Calderón EFIE MOT scheme on a torus with an inner and outer radii respectively equal to 0.2 m and 0.5 m, $N_s = 387$ and $\Delta t = 3$ ns. Near 1, four eigenvalues are clustered, superposed two by two on this figure.</figcaption></figure>

## B. A Convolution Quadrature Calderón time-domain EFIE regularized with quasi-Helmholtz projectors

The previous Calderón formulation is perfectly adapted to simply connected geometries, ensuring that the new operator
has no nullspace. However, on multiply-connected geome-
tries, the harmonic subspace is non-empty, thus enlarging the
nullspace of $\mathcal{T}$ which is a new source of DC instability in
(44)[33]. The discretized EFIE operators can be regularized
using the quasi-Helmholtz projectors to address this issue.
Because the regularization is based on projectors, it does
not compromise the $h$-refinement regularizing effect of the
original Calderón scheme. The regularized EFIE space-
discretized operators are

$$
\mathcal{T}^{\text{reg}} = \left( \frac{c_0}{a} \int_{-\infty}^{\cdot} \mathbf{P}^{\Lambda H} + \mathbf{P}^{\Sigma} \right) * \mathcal{T} * \left( \mathbf{P}^{\Lambda H} + \frac{a}{c_0} \frac{\partial}{\partial t} \mathbf{P}^{\Sigma} \right), \quad (46)
$$

$$
\mathcal{T}^{\text{reg}} = \left( \frac{c_0}{a} \int_{-\infty}^{\cdot} \mathbf{P}^{\Sigma H} + \mathbf{P}^{\Lambda} \right) * \mathcal{T} * \left( \mathbf{P}^{\Sigma H} + \frac{a}{c_0} \frac{\partial}{\partial t} \mathbf{P}^{\Lambda} \right), \quad (47)
$$

where the BC quasi-Helmholtz projectors are defined with
the loop-to-RWG transformation matrix $\Lambda$ [45] such that

$$
\mathbf{P}^{\Lambda} = (\mathbf{\Lambda}(\mathbf{\Lambda}^{T}\mathbf{\Lambda})^{+} \mathbf{\Lambda}^{T}) \quad \text{and} \quad \mathbf{P}^{\Sigma H} = \mathbf{I} - \mathbf{P}^{\Lambda}, \qquad (48)
$$

and where the scaling $c_0/a$ with $a$ defined as the maximal
diameter of the scatterer, ensures consistent dimensionality
and helps reduce the conditioning further. This application
of the projectors is equivalent to differentiating the non-
solenoidal contributions on the left and in time integrating
the solenoidal contributions on the right of each EFIE op-
erator [30]. The regularized Calderón operator in the space-
discretized time domain is

$$
\mathbf{T}_{\mathrm{c}}^{\mathrm{reg}} = \mathbf{T}^{\mathrm{reg}} \mathbf{G}_{\mathrm{m}}^{-1} * \mathbf{T}^{\mathrm{reg}}. \quad (49)
$$

At first sight, (46) and (47) seem to involve unpractical temporal integrals. However, the problematic contributions in the regularized EFIE operator $\mathbf{T}^{\text{reg}}$ will vanish, since $\mathbf{P}^{\Lambda H} \mathbf{T}_{\text{h}} = \mathbf{T}_{\text{h}} \mathbf{P}^{\Lambda H} = \mathbf{0}$ and $\mathbf{P}^{\Sigma} \mathbf{T}_{\text{h}} \mathbf{P}^{\Sigma} = \mathbf{T}_{\text{h}}$, and we have

$$
\begin{equation}
\begin{split}
\mathbf{T}^{\text{reg}} = {}& -a^{-1} \mathbf{P}^{\Lambda H} \mathbf{T}_s \mathbf{P}^{\Lambda H} - c_0^{-1} \mathbf{P}^{\Lambda H} \frac{\partial}{\partial t} \mathbf{T}_s \mathbf{P}^{\Sigma} \\
& - c_0^{-1} \mathbf{P}^{\Sigma} \frac{\partial}{\partial t} \mathbf{T}_s \mathbf{P}^{\Lambda H} - \frac{a}{c_0^2} \mathbf{P}^{\Sigma} \frac{\partial^2}{\partial t^2} \mathbf{T}_s \mathbf{P}^{\Sigma} + a \mathbf{T}_h.
\end{split}
\tag{50}
\end{equation}
$$

Similarly, the dual EFIE operator simplifies as

$$
\begin{equation}
\begin{split}
\mathbf{T}^{\text{reg}} = {}& -a^{-1} \mathbf{P}^{\Sigma H} \mathbf{T}_s \mathbf{P}^{\Sigma H} - c_0^{-1} \mathbf{P}^{\Sigma H} \frac{\partial}{\partial t} \mathbf{T}_s \mathbf{P}^{\Lambda} \\
& - c_0^{-1} \mathbf{P}^{\Lambda} \frac{\partial}{\partial t} \mathbf{T}_s \mathbf{P}^{\Sigma H} - \frac{a}{c_0^2} \mathbf{P}^{\Lambda} \frac{\partial^2}{\partial t^2} \mathbf{T}_s \mathbf{P}^{\Lambda} + a \mathbf{T}_h.
\end{split}
\tag{51}
\end{equation}
$$

The space-discretized formulation of the regularized
Calderón EFIE is

$$
\mathbf{T}_{\mathrm{c}}^{\mathrm{reg}} * \mathbf{y}_{\Gamma} = -\eta_{0}^{-1} \mathbf{T}^{\mathrm{reg}} \mathbf{G}_{\mathrm{m}}^{-1} * \mathbf{R} * \mathbf{e}_{\Gamma}^{\mathrm{inc}}, \quad (52)
$$

where

$$
\mathbf{R} = \frac{c_0}{a} \int_{-\infty}^{\cdot} \mathbf{P}^{\Lambda H} + \mathbf{P}^{\Sigma} \quad \text{and} \quad j_{\Gamma} = \left( \mathbf{P}^{\Lambda H} + \frac{a}{c_0} \frac{\partial}{\partial t} \mathbf{P}^{\Sigma} \right) y_{\Gamma}. \qquad (53)
$$

The right-hand side of (52) has a temporal integral which
is directly evaluated on the incident field to avoid quadratic

complexity with the number of time steps. This leads to the
following MOT

$$
\begin{equation}
\begin{split}
Z_{\mathbf{T}^{\text{reg}},0} y_i = -\eta_0^{-1} \sum_{j=0}^{N_{\text{conv}}} Z_{\mathbf{T}^{\text{reg}},j} \left( \tilde{\mathbf{P}}^{\Sigma} e_{i-j}^{\text{inc}} + \frac{c_0}{a} \tilde{\mathbf{P}}^{\Lambda H} e_{i-j}^{\text{prim}} \right) \\
\qquad - \sum_{j=1}^{N_{\text{conv}}} Z_{\mathbf{T}^{\text{reg}},j} y_{i-j},
\end{split}
\tag{54}
\end{equation}
$$

where $\tilde{\mathbf{P}}^{\Lambda H} = \mathbf{P}^{\Lambda H} \otimes I_p$, $\tilde{\mathbf{P}}^{\Sigma} = \mathbf{P}^{\Sigma} \otimes I_p$, the vector sequence $\mathbf{y}_n$ is the time discretization of $\mathbf{y}_\Gamma(t)$, and the time domain interaction matrices $Z_{\mathbf{T}^{\text{reg}},i}$ and $Z_{\mathbf{T}^{\text{reg}},i}$ are respectively generated by the convolution quadrature method described in Subsection B of the space-discretized operators $\mathbf{T}_c^{\text{reg}}$ and $\mathbf{T}^{\text{reg}} = \mathbf{T}^{\text{reg}}\mathbf{G}_m^{-1}$. Once the computation of $\mathbf{y}_n$ is done, the current $j$ still has to be evaluated. The convolution quadrature discretization of the time derivative is

$$
\begin{align}
Z_{\frac{\partial}{\partial t},i} &= Z^{-1} \{z \to \bar{s}_{\text{cq}}(z)\}_i \nonumber \\
&= \Delta t^{-1} \overline{\mathbf{A}^{-1}} \delta_{i,0} - \Delta t^{-1} \overline{\mathbf{A}^{-1} \mathbf{1}_p b^T \mathbf{A}^{-1}} \delta_{i,1}, \tag{55}
\end{align}
$$

where $\delta_{i,0}$ is the Kronecker delta, $\overline{\mathbf{A}^{-1}} = I_{N_s} \otimes \mathbf{A}^{-1}$ and $\overline{\mathbf{A}^{-1} \mathbf{1}_p b^T \mathbf{A}^{-1}} = I_{N_s} \otimes \mathbf{A}^{-1} \mathbf{1}_p b^T \mathbf{A}^{-1}$ [30]. Therefore, the current solution is obtained as

$$
\begin{align}
j_i &= \left[ \left( \tilde{\mathbf{P}}^{\Lambda H} + \frac{a}{c_0} Z_{\frac{\partial}{\partial t}} \tilde{\mathbf{P}}^{\Sigma} \right) *_s y_i \right]_i \tag{56} \\
&= \tilde{\mathbf{P}}^{\Lambda H} y_i + \frac{a}{\Delta t c_0} \tilde{\mathbf{P}}^{\Sigma} \left[ \overline{\mathbf{A}^{-1} y_i} - \overline{\mathbf{A}^{-1} \mathbf{1}_p b^T \mathbf{A}^{-1} y_i} \right]_{i-1}. \nonumber
\end{align}
$$

IV. Results

To test the effectiveness of the proposed schemes, simula-
tions have been realized with different geometries, excited
by a Gaussian pulse plane wave

$$
e^{\text{inc}}(\mathbf{r}, t) = A_0 \exp\left(-\frac{(t - \frac{\hat{k} \cdot \mathbf{r}}{c})^2}{2\sigma^2}\right) \hat{\mathbf{p}}, \quad (57)
$$

$$
\mathbf{h}^{\text{inc}}(\mathbf{r}, t) = \frac{1}{\eta_0} \hat{\mathbf{k}} \times \mathbf{e}^{\text{inc}}(\mathbf{r}, t), \quad (58)
$$

where $\sigma = 6/(2\pi f_{\text{bw}})$, $\hat{\mathbf{p}} = \hat{\mathbf{x}}$, $\hat{\mathbf{k}} = -\hat{\mathbf{z}}$, $A_0 = 1 \text{ V m}^{-1}$ and $f_{\text{bw}}$ is the frequency bandwidth. Notice that this frequency bandwidth is proportional to the maximal frequencies excited by the pulse Gaussian. In this work, the Runge-Kutta Radau IIA method of stage 2 is used for all simulations [17], [18]. The time step size $\Delta t$ has been chosen equal to $(\psi f_{\text{max}})^{-1}$ where $f_{\text{max}}$ is the upper frequency of the excitation and $\psi = 3$ is an oversampling parameter.

**A. Canonical geometries**

To illustrate the key properties of the newly proposed
schemes, namely the Calderón preconditioned MOT (44) and
the Calderón preconditioned MOT regularized by the quasi-
Helmholtz-projectors (54), they are compared in the case
of modelling of canonical scatterers to other formulations
present on the literature: the EFIE MOT schemes (MOT
EFIE) (21), the time-differentiated one (MOT TD-EFIE)
(24), the MOT regularized by the quasi-Helmholtz-projectors
<figure><img src="image_8.png" /></figure>

(MOT qH-EFIE) [30]. In this subsection, the excitation parameters have been chosen not to excite the first resonant mode of the geometries.

The first set of numerical tests were performed on the unit sphere, discretized with 270 RWG functions. The intensity of the resulting currents at one point of the geometry are shown in Figure 4. As expected, the time differentiated and the non-differentiated EFIEs are the only formulations suffering from DC-instabilities on this simply connected scenario. In addition, the condition number of the matrices to invert for each MOT are presented in Figure 5 and Figure 6 with respect to the time step size $\Delta t$ and the mesh density $h^{-1}$. The standard EFIE formulation and its time-differentiated counterpart suffer from ill-conditioning at large time steps while the stabilized ones remain well-conditioned. Instead, only the Calderón preconditioned formulations presented in this work remain well-conditioned for dense discretizations (Figure 6).

<figure><img src="image_9.png" /><figcaption>**FIGURE 4.** Evolution in time of the current intensity at $r_0 = (-0.36, 0.89, 0.11)$ m of the sphere where $f_{bw} = 25$ kHz with parameters $N_s = 270$ and $\Delta t = 4.5 \times 10^4$ ns.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**FIGURE 5.** Condition number of the EFIE MOT matrices with respect to the time step size $\Delta t$ ($N_s = 270$) on a sphere.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 6.** Condition number of the EFIE MOT matrices with respect to the mesh size $h$ ($\Delta t = 45$ ns) on a sphere.</figcaption></figure>

The second set of tests focused on the stability of the different formulations when modelling multiply connected scatters, here a torus with inner radius of 0.2 m and outer radius of 0.5 m. The current densities at the probe point are shown in Figure 7 and the conditioning studies are represented in Figure 8 and Figure 9. In line with the polynomial eigenvalue analysis of the non-regularized Calderón EFIE formulation (Figure 3), the formulation (44) suffers from DC instability. Moreover, the static nullspace of the continuous operator deteriorates the condition number of the matrix to invert for large time steps (Figure 9). However, the newly proposed regularized Calderón formulation (54) is stable and remains well-conditioned at large time steps and dense meshes for this geometry.

<figure><img src="image_12.png" /><figcaption>**FIGURE 7.** Evolution in time of the current intensity at $r_0 = (0.62, -0.13, 0.11)$ m of the torus where $f_{bw} = 25$ kHz with parameters $N_s = 387$ and $\Delta t = 4.5 \times 10^4$ ns.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 8. Condition number with respect to the mesh size $h$ ($\Delta t$ = 4.5 ns) on a torus of the EFIE MOT schemes.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 9. Condition number with respect to the time step size $\Delta t$ ($N_s$ = 387) on a torus of the EFIE MOT schemes.</figcaption></figure>

## B. Non-canonical geometries

The final set of numerical tests is dedicated to more complex test structures (Figure 10). In addition, instead of direct solver we rely on the iterative solver GMRES with different relative target tolerances $\epsilon$ [53]. For practical reason, the maximum number of iteration has been limited to 200 without restart.

All the structures have been illuminated by a pulse Gaussian plane wave with $f_{bw} = 1.6$ MHz. Table 1 shows the condition number and the number of iterations needed. The Calderón preconditioned formulation (CP-EFIE) (44) and the Calderón preconditioned formulation regularized by the quasi-Helmholtz-projectors (qH-CP-EFIE) (54) are the only formulations requiring less than 200 iterations to converge at each time steps. As expected, the condition number of the CP-EFIE is high on non-simply connected geometries because of the presence of the operator DC nullspace on these structures. Although, in this case, the number of iterations remains low because the nullspace results in isolated elements in the spectrum, the solution is corrupted by DC instability arising from this nullspace. This phenomenon is absent for the qH-CP-EFIE formulation which is free from high conditioning or DC instability and yields stable solutions up to the target precision of the iterative solver.

Finally, the setup time of the different formulations for the cube with multiple holes scatterer and the corresponding times required to compute a single time step are listed

<figure><img src="image_15.png" /><figcaption>FIGURE 10. Test structures. (a) space shuttle; (b) compression spring; (c) cube with multiple holes.</figcaption></figure>

<table><caption>Space shuttle: $N_s$ = 4311</caption><thead><tr><th rowspan="2">Tolerance $\epsilon$</th><th>Cond(Z<sub>0</sub>)</th><th colspan="3">N<sub>iter</sub></th></tr><tr><th>Ø</th><th>10<sup>-3</sup></th><th>10<sup>-6</sup></th><th>10<sup>-10</sup></th></tr></thead><tbody><tr><td>MOT TD-EFIE</td><td>4.3 · 10<sup>6</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>MOT EFIE</td><td>6.5 · 10<sup>6</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>MOT qH-EFIE</td><td>5.3 · 10<sup>4</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>This work (44)</td><td>38</td><td>15</td><td>28</td><td>41</td></tr><tr><td>This work (54)</td><td>41</td><td>4</td><td>18</td><td>33</td></tr></tbody></table>

<table><caption>Compression spring: $N_s$ = 3906</caption><thead><tr><th rowspan="2">Tolerance $\epsilon$</th><th>Cond(Z<sub>0</sub>)</th><th colspan="3">N<sub>iter</sub></th></tr><tr><th>Ø</th><th>10<sup>-3</sup></th><th>10<sup>-6</sup></th><th>10<sup>-10</sup></th></tr></thead><tbody><tr><td>MOT TD-EFIE</td><td>2.9 · 10<sup>6</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>MOT EFIE</td><td>4.3 · 10<sup>6</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>MOT qH-EFIE</td><td>1.6 · 10<sup>3</sup></td><td>57</td><td>154</td><td>&gt; 200</td></tr><tr><td>This work (44)</td><td>158</td><td>8</td><td>31</td><td>50</td></tr><tr><td>This work (54)</td><td>133</td><td>2</td><td>19</td><td>39</td></tr></tbody></table>

<table><caption>TABLE 1. Condition number and maximum number of iterations of the different formulations with different relative tolerances.</caption><caption>Cube with multiple holes: $N_s$ = 3267</caption><thead><tr><th rowspan="2">Tolerance $\epsilon$</th><th>Cond(Z<sub>0</sub>)</th><th colspan="3">N<sub>iter</sub></th></tr><tr><th>Ø</th><th>10<sup>-3</sup></th><th>10<sup>-6</sup></th><th>10<sup>-10</sup></th></tr></thead><tbody><tr><td>MOT TD-EFIE</td><td>1.2 · 10<sup>9</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>MOT EFIE</td><td>1.8 · 10<sup>9</sup></td><td>&gt; 200</td><td>&gt; 200</td><td>&gt; 200</td></tr><tr><td>MOT qH-EFIE</td><td>6.6 · 10<sup>3</sup></td><td>50</td><td>196</td><td>&gt; 200</td></tr><tr><td>This work (44)</td><td>9.5 · 10<sup>6</sup></td><td>5</td><td>20</td><td>82</td></tr><tr><td>This work (54)</td><td>50</td><td>1</td><td>11</td><td>25</td></tr></tbody></table>
<figure><img src="image_16.png" /></figure>

<table><caption>TABLE 2. Computation times of the different formulations with $\epsilon = 10^{-6}$.</caption><thead><tr><th colspan="4">Cube with holes: N<sub>s</sub> = 3267</th></tr><tr><th></th><th>Setup time</th><th>Computation time for<br>a single time step</th><th>DC stable</th></tr></thead><tbody><tr><td>MOT TD-EFIE</td><td>10 min</td><td>69 s</td><td>No</td></tr><tr><td>MOT EFIE</td><td>10 min</td><td>78 s</td><td>No</td></tr><tr><td>MOT qH-EFIE</td><td>19 min</td><td>3.4 s</td><td>Yes</td></tr><tr><td>This work (44)</td><td>134 min</td><td>0.84 s</td><td>No</td></tr><tr><td>This work (54)</td><td>138 min</td><td>0.40 s</td><td>Yes</td></tr></tbody></table>

in Table 2. As expected, the preconditioned formulations require a longer setup time, which is rapidly compensated by a faster computation of the time steps. The setup time of these formulations can be significantly reduced by leveraging fast solvers [54]; this, however, goes beyond the scope of this work.

## Conclusion

In this paper, novel Calderón preconditioned techniques have been presented for the time domain Electric Field Integral Equations solved with Marching-On-In-Time with convolution quadratures. These formulations eliminate the DC-instability for simply and multiply connected geometries. In addition, they cure the *h*-refinement and large time step breakdowns and generate well-conditioned Marching-On-In Time. Finally, numerical results on complex geometries showcased the effectiveness of the proposed schemes.

## References

[1] A. F. Peterson, S. L. Ray, R. Mittra, I. of Electrical, and E. Engineers, *Computational methods for electromagnetics*. IEEE press New York, 1998, vol. 351.

[2] J.-M. Jin, *Theory and computation of electromagnetic fields*. John Wiley & Sons, 2011.

[3] W. C. Gibson, *The method of moments in electromagnetics*. Chapman and Hall/CRC, 2021.

[4] D. Colton and R. Kress, *Integral equation methods in scattering theory*. SIAM, 2013.

[5] C. L. Bennett Jr, *A technique for computing approximate electromagnetic impulse response of conducting bodies*. Purdue University, 1968.

[6] S. M. Rao and D. R. Wilton, "Transient scattering by conducting surfaces of arbitrary shape," *IEEE Transactions on Antennas and Propagation*, vol. 39, no. 1, pp. 56–61, 1991.

[7] F.-J. Sayas, *Retarded potentials and time domain boundary integral equations: A road map*. Springer, 2016, vol. 50.

[8] Y.-S. Chung, T. K. Sarkar, B. H. Jung, et al., "Solution of time domain electric field integral equation using the laguerre polynomials," *IEEE Transactions on Antennas and Propagation*, vol. 52, no. 9, pp. 2319–2328, 2004.

[9] C. Lubich, "Convolution quadrature and discretized operational calculus. i," *Numerische Mathematik*, vol. 52, no. 2, pp. 129–145, 1988.

[10] C. Lubich, "Convolution quadrature and discretized operational calculus. ii," *Numerische Mathematik*, vol. 52, no. 4, pp. 413–425, 1988.

[11] M. Schanz and H. Antes, "A new visco-and elastodynamic time domain boundary element formulation," *Computational Mechanics*, vol. 20, no. 5, pp. 452–459, 1997.

[12] L. Banjai, M. Messner, and M. Schanz, "Runge-kutta convolution quadrature for the boundary element method," *Computer methods in applied mechanics and engineering*, vol. 245, pp. 90–101, 2012.

[13] X. Wang and D. S. Weile, "Implicit runge-kutta methods for the discretization of time domain integral equations," *IEEE transactions on antennas and propagation*, vol. 59, no. 12, pp. 4651–4663, 2011.

[14] C. Lubich and A. Ostermann, "Runge-kutta methods for parabolic equations and convolution quadrature," *mathematics of computation*, vol. 60, no. 201, pp. 105–131, 1993.

[15] A. V. Oppenheim, J. R. Buck, and R. W. Schafer, *Discrete-time signal processing*. Vol. 2. Upper Saddle River, NJ: Prentice Hall, 2001.

[16] F. Chipman, "A-stable runge-kutta processes," *BIT Numerical Mathematics*, vol. 11, no. 4, pp. 384–388, 1971.

[17] C. Runge, "Über die numerische auflösung von differentialgleichungen," *Mathematische Annalen*, vol. 46, no. 2, pp. 167–178, 1895.

[18] W. Kutta, "Beitrag zur naherungsweisen integration totaler differentialgleichungen," *Z. Math. Phys.*, vol. 46, pp. 435–453, 1901.

[19] J. C. Butcher, *Numerical methods for ordinary differential equations*. John Wiley & Sons, 2016.

[20] L. Skvortsov, "Accuracy of runge-kutta methods applied to stiff problems," *Computational mathematics and mathematical physics*, vol. 43, no. 9, pp. 1320–1330, 2003.

[21] L. Banjai and M. Kachanovska, "Fast convolution quadrature for the wave equation in three dimensions," *Journal of Computational Physics*, vol. 279, pp. 103–126, 2014.

[22] T. Maruyama, T. Saitoh, T. Bui, and S. Hirose, "Transient elastic wave analysis of 3-d large-scale cavities by fast multipole bem using implicit runge-kutta convolution quadrature," *Computer Methods in Applied Mechanics and Engineering*, vol. 303, pp. 231–259, 2016.

[23] G. H. Golub and C. F. Van Loan, *Matrix computations*. JHU press, 2013.

[24] F. P. Andriulli, H. Bagci, F. Vipiana, G. Vecchi, and E. Michielssen, "Analysis and regularization of the td-efie low-frequency breakdown," *IEEE transactions on antennas and propagation*, vol. 57, no. 7, pp. 2034–2046, 2009.

[25] N.-W. Chen, K. Aygün, and E. Michielssen, "Integral-equation-based analysis of transient scattering and radiation from conducting bodies at very low frequencies," *IEE Proceedings-Microwaves, Antennas and Propagation*, vol. 148, no. 6, pp. 381–387, 2001.

[26] F. P. Andriulli, K. Cools, I. Bogaert, and E. Michielssen, "On a well-conditioned electric field integral operator for multiply connected geometries," *IEEE transactions on antennas and propagation*, vol. 61, no. 4, pp. 2077–2087, 2012.

[27] I. Bogaert, K. Cools, F. P. Andriulli, and H. Bağcı, "Low-frequency scaling of the standard and mixed magnetic field and müller integral equations," *IEEE Transactions on Antennas and Propagation*, vol. 62, no. 2, pp. 822–831, 2013.

[28] F. P. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, "Time domain calderón identities and their application to the integral equation analysis of scattering by pec objects part ii: Stability," *IEEE transactions on antennas and propagation*, vol. 57, no. 8, pp. 2365–2375, 2009.
[29] K. Cools, F. P. Andriulli, F. Olyslager, and E. Michielssen, “Time domain calderón identities and their application to the integral equation analysis of scattering by pec objects part i: Preconditioning,” *IEEE Transactions on Antennas and Propagation*, vol. 57, no. 8, pp. 2352–2364, 2009.

[30] A. Dély, F. P. Andriulli, and K. Cools, “Large time step and dc stable td-efie discretized with implicit runge-kutta methods,” *IEEE Transactions on Antennas and Propagation*, vol. 68, no. 2, pp. 976–985, 2019.

[31] Y. Shi, H. Bağcı, and M. Lu, “On the static loop modes in the marching-on-in-time solution of the time-domain electric field integral equation,” *IEEE Antennas and Wireless Propagation Letters*, vol. 13, pp. 317–320, 2014.

[32] F. P. Andriulli, K. Cools, H. Bagci, et al., “A multiplicative calderon preconditioner for the electric field integral equation,” *IEEE Transactions on Antennas and Propagation*, vol. 56, no. 8, pp. 2398–2412, 2008.

[33] K. Cools, F. P. Andriulli, F. Olyslager, and E. Michielssen, “Nullspaces of mfie and calderón preconditioned efie operators applied to toroidal surfaces,” *IEEE Transactions on Antennas and Propagation*, vol. 57, no. 10, pp. 3205–3215, 2009.

[34] Y. Beghein, K. Cools, and F. P. Andriulli, “A dc-stable, well-balanced, calderón preconditioned time domain electric field integral equation,” *IEEE Transactions on Antennas and Propagation*, vol. 63, no. 12, pp. 5650–5660, 2015.

[35] P. Cordel, A. Dély, A. Merlini, and F. P. Andriulli, “Calderón preconditioners for the td-efie discretized with convolution quadratures,” in *2022 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting (AP-S/URSI)*, IEEE, 2022, pp. 1672–1673.

[36] R. F. Harrington and J. L. Harrington, *Field computation by moment methods*. Oxford University Press, Inc., 1996.

[37] S. Rao, D. Wilton, and A. Glisson, “Electromagnetic scattering by surfaces of arbitrary shape,” *IEEE Transactions on antennas and propagation*, vol. 30, no. 3, pp. 409–418, 1982.

[38] L. Banjai and F. Sayas, “Integral equation methods for evolutionary pde,” *Springer Ser. Comput. Math.*, Springer, Cham, 2022.

[39] L. Schwartz, *Mathematics for the physical sciences*. Courier Dover Publications, 2008.

[40] A. Dély, A. Merlini, K. Cools, and F. P. Andriulli, “Convolution quadrature time domain integral equation methods for electromagnetic scattering,” *Advances in Time-Domain Computational Electromagnetic Methods*, pp. 321–359, 2022.

[41] P.-A. Raviart and J.-M. Thomas, “A mixed finite element method for 2-nd order elliptic problems,” in *Mathematical aspects of finite element methods*, Springer, 1977, pp. 292–315.

[42] Y. Beghein, K. Cools, and D. De Zutter, “A temporal galerkin discretization of the charge-current continuity equation,” in *2013 International Conference on Electromagnetics in Advanced Applications (ICEAA)*, IEEE, 2013, pp. 628–631.

[43] P. W. van Diepen, R. J. Dilz, A. P. Zwamborn, and M. C. van Beurden, “The role of jordan blocks in the mot-scheme time domain efie linear-in-time solution instability,” *Progress In Electromagnetics Research B*, vol. 95, no. 95, p. 123, 2022.

[44] S. Walker, M. Bluck, and I. Chatzis, “The stability of integral equation time-domain scattering computations for three-dimensional scattering; similarities and differences between electrodynamic and elastodynamic computations,” *International Journal of Numerical Modelling: Electronic Networks, Devices and Fields*, vol. 15, no. 5-6, pp. 459–474, 2002.

[45] S. B. Adrian, A. Dély, D. Consoli, A. Merlini, and F. P. Andriulli, “Electromagnetic integral equations: Insights in conditioning and preconditioning,” *IEEE Open Journal of Antennas and Propagation*, 2021.

[46] S. B. Adrian, A. Dély, D. Consoli, A. Merlini, and F. P. Andriulli, “Electromagnetic integral equations: Insights in conditioning and preconditioning,” *IEEE Open Journal of Antennas and Propagation*, vol. 2, pp. 1143–1174, 2021. DOI: 10.1109/OJAP.2021.3121097.

[47] J.-C. Nédélec, *Acoustic and electromagnetic equations: integral representations for harmonic problems*. Springer, 2001, vol. 144.

[48] G. C. Hsiao and R. E. Kleinman, “Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics,” *IEEE transactions on Antennas and Propagation*, vol. 45, no. 3, pp. 316–328, 1997.

[49] F. P. Andriulli, K. Cools, I. Bogaert, H. Bagci, P. Yla-Ojiala, and E. Michielssen, “Analysis and discretization of the yukawa-calderon preconditioned cfie,” in *28th Annual Review of Progress in Applied Computational Electromagnetics (ACES-2012)*, Curran Associates, Inc., 2012, pp. 454–463.

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

[51] K. Cools, F. Andriulli, D. De Zutter, and E. Michielssen, “Accurate and conforming mixed discretization of the mfie,” *IEEE antennas and wireless propagation letters*, vol. 10, pp. 528–531, 2011.

[52] F. P. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, “The “dottrick tdefie”: A dc stable integral equation for analyzing transient scattering from pec bodies,” in *2008 IEEE Antennas and Propagation Society International Symposium*, IEEE, 2008, pp. 1–4.

[53] Y. Saad and M. H. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems,” *SIAM Journal on scientific and statistical computing*, vol. 7, no. 3, pp. 856–869, 1986.

[54] H. Bagci, F. P. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, “A calderón multiplicative preconditioner for the combined field integral equation,” *IEEE Transactions on Antennas and Propagation*, vol. 57, no. 10, pp. 3387–3392, 2009. DOI: 10.1109/TAP.2009.2029389.

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

**PIERRICK CORDEL** (Student Member, IEEE) received the M.Sc. Eng. Degree from the Ecole Nationale Supérieure des Mines de Nancy, France, in 2021 as well as a M.Sc. in Industrial Mathematics from the University of Luxembourg, Luxembourg, the same year. Currently, he is doing a PhD thesis at the institute Politecnico di Torino, in Italy.

His research focuses on time domain integral equations discretized with the convolution quadra-ture method.
<figure><img src="image_18.png" /></figure>

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

**ALEXANDRE DÉLY** received the M.Sc. Eng. degree from the École Nationale Supérieure des Télécommunications de Bretagne (Télécom Bretagne), France, in 2015. He received the Ph.D. degree from the École Nationale Supérieure Mines-Télécom Atlantique (IMT Atlantique), France, and from the University of Nottingham, United Kingdom, in 2019.

His research focuses on preconditioned and fast solution of boundary element methods frequency domain and time domain integral equations. He is currently working in Thales, Elancourt, France, on electromagnetic modeling and numerical simulations

Chief of the IEEE Antennas and Propagation Magazine and he serves as a Track Editor for the IEEE Transactions on Antennas and Propagation, and as an Associate Editor for URSI Radio Science Letters. He served as an Associate Editor for the IEEE Transactions on Antennas and Propagation, IEEE Antennas and Wireless Propagation Letters, IEEE Access and IET-MAP.

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

**ADRIEN MERLINI** (S'16–M'19) received the M.Sc.Eng. degree from the École Nationale Supérieure des Télécommunications de Bretagne (Télécom Bretagne), France, in 2015 and received the Ph.D. degree from the École Nationale Supérieure Mines Télécom Atlantique (IMT Atlantique), France, in 2019. From 2018 to 2019, he was a visiting Ph.D. student at the Politecnico di Torino, Italy, which he then joined as a Research Associate. Since 2019, he has been an Associate Professor with the Microwave Department, IMT

Atlantique. His research interests include preconditioning and acceleration of integral equation solvers for electromagnetic simulations and their application in brain imaging.

Dr. Merlini received 2 Young Scientist Awards at the URSI GASS 2020 and the EMTS 2023 meetings. In addition, he has co-authored a paper that received the 2022 ICEAA-IEEE APWC Best Paper Award, 5 that received honorable mentions (URSI/IEEE-APS 2021, 2022, and 2023) and 3 best paper finalists (URSI GASS 2020, URSI/IEEE-APS 2021 and 2022). He is a member of IEEE-HKN, the IEEE Antennas and Propagation Society, URSI France, and of the Lab-STICC laboratory. He is currently serving as Associate Editor for the Antenna and Propagation Magazine.

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

**FRANCESCO P. ANDRIULLI** (S'05–M'09–SM'11–F'23) received the Laurea in electrical engineering from the Politecnico di Torino, Italy, in 2004, the MSc in electrical engineering and computer science from the University of Illinois at Chicago in 2004, and the PhD in electrical engineering from the University of Michigan at Ann Arbor in 2008. From 2008 to 2010 he was a Research Associate with the Politecnico di Torino. From 2010 to 2017 he was an Associate Professor (2010-2014) and then Full Professor with the

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

Prof. Andriulli received several best paper awards at conferences and symposia (URSI NA 2007, IEEE AP-S 2008, ICEAA IEEE-APWC 2015) also in co-authorship with his students and collaborators (ICEAA IEEE-APWC 2021, EMTS 2016, URSI-DE Meeting 2014, ICEAA 2009) with whom received also a second prize conference paper (URSI GASS 2014), a third prize conference paper (IEEE-APS 2018), seven honorable mention conference papers (ICEAA 2011, URSI/IEEE-APS 2013, 4 in URSI/IEEE-APS 2022, URSI/IEEE-APS 2023) and other three finalist conference papers (URSI/IEEE-APS 2012, URSI/IEEE-APS 2007, URSI/IEEE-APS 2006, URSI/IEEE-APS 2022). A Fellow of the IEEE, he is also the recipient of the 2014 IEEE AP-S Donald G. Dudley Jr. Undergraduate Teaching Award, of the triennium 2014-2016 URSI Issac Koga Gold Medal, and of the 2015 L. B. Felsen Award for Excellence in Electrodynamics.

Prof. Andriulli is a member of Eta Kappa Nu, Tau Beta Pi, Phi Kappa Phi, and of the International Union of Radio Science (URSI). He is the Editor-in-