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

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.


I. INTRODUCTION
T IME 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 spacediscretized 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].
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 ℎ-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 matrixvector 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].

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 ( inc ,  inc )(, ) which induces a surface current density   on its boundary  that is the solution of the time-domain EFIE Here, n is the outpointing normal to  and  0 is the characteristic impedance of the background.The electric field operator T includes the contributions of the vector and scalar potentials, respectively denoted T s and T h [30] T  (, ) = − 1 where  0 is the speed of light in the background.The temporal convolution product *  and the temporal Green function G are defined as with  the time Dirac delta.

B. Marching-On-In-Time with Convolution Quadratures
Let  be a placeholder for any of the integral operators previously presented and let  (, ) be a causal function (∀ < 0,  (, ) = 0).With these notations, most time domain integral equation take the form where  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  c is expanded as a linear combination of   spatial basis functions such that [38]  c (, ) ≈ where  src are the source spatial basis functions and their associated time coefficients are stored in the vector f  ().
Then, the equation ( 7) is tested by the spatial basis functions  tst leading to the time-dependent matrix system where for  and  in ⟦1,   ⟧, we have, with The second step is the discretization in time with the convolution 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  and k  .The system (9) is then equivalent to Then, a representation on the Z-domain discretizes the system (12).The Laplace parameter, in the operator θ L , is replaced by the matrix-valued parameter diagonalizable for the considered  values, with the following eigenvalue decomposition s cq () = Q()()Q −1 ().The elements of the diagonal matrix () are the eigenvalues of s cq () and the columns of Q() are their associated eigenvectors.The time step size is denoted  and the matrix A and the vectors 1  , c, b of size  are determined by the implicit scheme used [40].The discretized Z-domain operator θ Z is defined such that for any for  and  in ⟦1, ⟧ and for  and  in ⟦1, where  , = (−1)+ is an appropriate indexing function and θ  , () is a diagonal matrix defined as The vectors f Z and k Z are the Z-domain representation [40] of the respective time-discretized vectors yielding to the following discretization of the system (12) Finally, the equivalent time discretized system of ( 9) is obtained by applying the inverse Z-transform on ( 17) where Z θ, = Z −1  → θ Z ()  are the time domain interaction matrices and *  is the sequence convolution product.The system sequence ( 18) is rewritten in the following Marching-On-In-Time that can be solved for

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  rwg [37], [41] are used to expand the current density as where the current coefficients are gathered in an unknown vector function of time j  ().The EFIE is then tested with rotated RWG basis functions n ×  rwg    , leading to the following Marching-On-In-Time where the vector sequences j  and e inc  , and the time domain interaction matrices Z T , are respectively generated by the convolution quadrature method described in Subsection B of j  () and the following space-discretized vector and matrix However, the time integral contribution of this operator T involves an unbounded number of non-vanishing matrices Z T , (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 T ,0 j  = − −1 0 e inc  − where e inc  and Z  are respectively the time domain vectors and interaction matrices generated by the convolution quadrature method described in Subsection B of

D. EFIE DC instability
The electric field integral operator suffers from the DC instability: since for all constant-in-time solenoidal current  cs we have ∇ • j cs = 0 and    cs = 0, we can conclude that T  cs = 0 .
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 currents 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.

E. Quasi-Helmholtz projectors
Previous works show that the electric field integral equation discretized in space using RWG basis functions can be stabilized 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/quasiharmonic space) are respectively 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 issues.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 continuous 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] where • is the composition operator, I is the identity operator and K is defined as The operator −I/4 + K 2 is a well-behaved operator for increasing discretization densities.As a consequence, with a proper discretization, T 2 is well-conditioned for large time steps and dense meshes for simply connected structures [45].
In practice, a discretization of 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  bc The preconditioning leads to the following space-discretized formulation where the matrix G  is the mixed gram matrix linking the the two discretizations Then, the convolution quadrature leads to the MOT scheme where Z T, are the time domain interaction matrices of the space-discretized operator T () (31) generated by the convolution quadrature method described in Subsection B, the sequence convolution quadrature product *  is the discretization of the space-discretized temporal convolution product * and The Kronecker product ⊗I  is required to match with the convolution quadrature method where I  is the identity matrix of size .Unfortunately, the MOT in (34), involves operators with temporal integrations leading to a time consuming MOT.A more favorable scheme can be obtained by noticing the following commutative properties and the cancellation property T 2 h = 0 [32], we have This is advantageous since, besides not involving any time integration contribution, the operator  −2 0  2  2 T 2 s − T s T h − T h T s has no nullspace for simply connected geometries leading to a DC-stable discretization ("dottrick TDEFIE") [52].By extending the previous notations on T s and T h the proposed space-discretized operator is denoted ) yielding to the following space-discretized formulation 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 T ℎ is evaluated with the incident field where Therefore, the previous MOT is rewritten as where the time domain interaction matrices Z T c , , Z T s , and Z T h , are respectively generated by the convolution quadrature method described in Subsection B of the space- As in (34), the interaction matrix sequence Z T c , involves computationally expensive sequence convolution products *  , however, the convolution quadrature method allows the substitution of the sequence convolution products *  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 , is equal to where the matrix s cq () = I   ⊗ s cq () is the Z-discretization of the time derivative and I   is the identity matrix of size   .The formulation ( 44) is a good candidate to obtain a stable current solution, however, the proposed operator 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].

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 geometries, the harmonic subspace is non-empty, thus enlarging the nullspace of 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 ℎ-refinement regularizing effect of the original Calderón scheme.The regularized EFIE spacediscretized operators are where the BC quasi-Helmholtz projectors are defined with the loop-to-RWG transformation matrix Λ [45] such that and where the scaling  0  with  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 nonsolenoidal contributions on the left and in time integrating the solenoidal contributions on the right of each EFIE operator [30].The regularized Calderón operator in the spacediscretized time domain is At first sight, ( 46) and ( 47) seem to involve unpractical temporal integrals.However, the problematic contributions in the regularized EFIE operator T reg will vanish, since P  T h = T h P  = 0 and P  T h P  = T h , and we have Similarly, the dual EFIE operator simplifies as The space-discretized formulation of the regularized Calderón EFIE is where 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 where P  = P  ⊗   , P  = P  ⊗   , the vector sequence y  is the time discretization of y  (), and the time domain interaction matrices Z T reg c , and Z T reg , are respectively generated by the convolution quadrature method described in Subsection B of the space-discretized operators T reg c and Once the computation of y  is done, the current j still has to be evaluated.The convolution quadrature discretization of the time derivative is where  ,0 is the Kronecker delta, [30].Therefore, the current solution is obtained as

IV. Results
To test the effectiveness of the proposed schemes, simulations have been realized with different geometries, excited by a Gaussian pulse plane wave where  = 6/(2  bw ), p = x, k = − ẑ,  0 = 1 V m −1 and  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  has been chosen equal to (  max ) −1 where  max is the upper frequency of the excitation and  = 3 is an oversampling parameter.

A. Canonical geometries
To illustrate the key properties of the newly proposed schemes, namely the Calderón preconditioned formulation (44) and the and the Calderón preconditioned formulation 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 formulation regularized by the quasi-Helmholtz-projectors (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  and the mesh density ℎ −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).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 ( 52) is stable and remains well-conditioned at large time steps and dense meshes for this geometry.

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  [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  bw = 1.6MHz.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) (52) 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.Even if, in this case, the number of iterations remains low, 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.

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 DCinstability for simply and multiply connected geometries.In addition, they cure the ℎ-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.

FIGURE 3 .
FIGURE 3. Polynomial eigenvalues of the Calder ón EFIE MOT scheme on a torus with a inner and outer radii respectively equal to 0.2 m and 0.5 m,   = 387 and  = 3 ns.Near 1, four eigenvalues are clustered, superposed two by two on this figure.

FIGURE 8 .
FIGURE 8. Condition number with respect to the mesh size ℎ ( = 4.5 ns) on a torus of the EFIE MOT schemes.