A Coupled Hybridizable Discontinuous Galerkin and Boundary Integral Method for Analyzing Electromagnetic Scattering

A coupled hybridizable discontinuous Galerkin (HDG) and boundary integral (BI) method is proposed to efficiently analyze electromagnetic scattering from inhomogeneous/composite objects. The coupling between the HDG and the BI equations is realized using the numerical flux operating on the equivalent current and the global unknown of the HDG. This approach yields sparse coupling matrices upon discretization. Inclusion of the BI equation ensures that the only error in enforcing the radiation conditions is the discretization. However, the discretization of this equation yields a dense matrix, which prohibits the use of a direct matrix solver on the overall coupled system as often done with traditional HDG schemes. To overcome this bottleneck, a"hybrid"method is developed. This method uses an iterative scheme to solve the overall coupled system but within the matrix-vector multiplication subroutine of the iterations, the inverse of the HDG matrix is efficiently accounted for using a sparse direct matrix solver. The same subroutine also uses the multilevel fast multipole algorithm to accelerate the multiplication of the guess vector with the dense BI matrix. The numerical results demonstrate the accuracy, the efficiency, and the applicability of the proposed HDG-BI solver.


Introduction
In the past two decades, the discontinuous Galerkin (DG) method [1][2][3][4][5][6][7][8][9][10][11][12][13][14] has attracted significant attention in the computational electromagnetics research community because, compared the traditional finite element method (FEM) [15][16][17][18][19], it offers a higher-level of flexibility in discretization which allows for non-conformal meshes and an easier implementation of hand/or p-adaptivity.In addition, in time domain, when combined with an explicit time integration scheme, the DG method produces a very compact, fast, and easy-to-parallelize solver since the DG's block diagonal mass matrix is inverted once and very efficiently before the time marching is started.However, this increased efficiency does not carry over to the frequency domain.Due to doubling of the unknowns at the element boundaries and the fact that a sparse matrix system must still be solved, frequency-domain DG schemes usually require more computational resources than the traditional frequency-domain FEM.
Recently, this drawback has been alleviated with the introduction of the hybridizable discontinuous Galerkin (HDG) method [20].HDG introduces single-valued hybrid variables on the skeleton of the mesh (namely, a mesh that consists of only the faces of the elements) [20] and converts the local/elemental DG matrix systems into a coupled global matrix system, where these hybrid variables are the unknowns to be solved for.The computational requirements of HDG are lower than DG since the total degrees of freedom is now reduced [20].Indeed, HDG is competitive to FEM in terms of computational requirements when both methods use a high-order discretization, and at the same time, it maintains the advantages of the traditional DG over FEM [21,22].In addition, thanks to the local post-processing used after the global matrix system solution, HDG achieves an accuracy convergence of order p + 2 (superconvergence), where p is the order of polynomial basis functions used to expand the local/elemental field variables [23].Because of these benefits, HDG has been used to solve various equations of physics, such as convection-diffusion equations [24], Poisson equa-tion [25], and elastic/acoustic wave equations [26].For electromagnetics, HDG was first used to solve the two-dimensional (2D) Maxwell's equations [27].Since then, it has been extended to solve the three-dimensional (3D) Maxwell's equations and used in conjunction with a Schwarz-type domain decomposition method to analyze electromagnetic scattering from large objects [28,29].In addition, a hybridizable discontinuous Galerkin time-domain method (HDGTD) has been proposed to solve the time-dependent Maxwell's equations.This method combines an implicit and explicit time integration scheme and HDG for time marching and spatial discretization, respectively.[30][31][32].HDG has also been used in simulations of multiphysics problems: In [33,34], the coupled system of the Maxwell's equations and the hydrodynamic equation has been solved using HDG to simulate the non-local optical response of nanostructures.
Most of the HDG methods, which have been developed to simulate wave interactions, use approximate absorbing boundary conditions (ABCs) to truncate the computation domain [35][36][37][38][39].Although these boundary conditions yield sparse matrices upon discretization, their accuracy is limited and therefore they restrict the high-order convergence of the solution unless a very large computation domain is used.One can also use the method of perfectly matched layer (PML) to truncate the HDG computation domain [40][41][42][43].Indeed HDG with PML has recently been used in the simulation of waveguide transmission problems [44].However, to increase the "absorption" of PML (i.e., to increase its accuracy), one has to increase thickness of the layer or the value of the conductivity.The first option increases the size of the computation domain while the second option has to be done carefully since large values of conductivity often result in numerical reflection from the PML-computation domain interface and decrease the accuracy of the solution [14,45].
On the other hand, boundary integral (BI)-based approaches to truncating computation domains do not suffer from these bottlenecks [46][47][48][49][50][51][52][53][54].In this work, HDG is used together with a BI formulation to efficiently and accurately simulate electromagnetic scattering from However, the discretization of the BI equation yields a dense matrix, which prohibits the use of a direct matrix solver on the overall coupled system as often done with traditional HDG schemes [28,29].To overcome this bottleneck, in this work, a "hybrid" method is developed.This method uses an iterative scheme to solve the overall coupled system but within the matrix-vector multiplication subroutine of the iterations, the inverse of the HDG matrix is efficiently accounted for using a sparse direct matrix solver.The same subroutine also uses the multilevel fast multipole algorithm (MLFMA) [55][56][57][58][59][60][61][62][63] to accelerate the multiplication of the guess vector with the dense BI matrix.Another contribution of this work is that it describes in detail the first use of vector basis functions [64] within the HDG framework.
The rest of this paper is organized as follows.Section II first describes the electromagnetic scattering problem and introduces the mesh used to discretize the computation domain.This is followed by the formulation of the coupled HDG and BI equations and the description of the matrix system that is obtained by discretizing them.Finally, Section II introduces the hybrid scheme developed to efficiently solve this matrix system.Section III provides several numerical examples to demonstrate the computational benefits of the proposed HDG-BI solver.In Section IV, a short summary of the work is provided and several future research directions are briefly described.

Problem Description
Consider the electromagnetic scattering problem involving a dielectric object that resides in an unbounded background medium with permittivity ε 0 and permeability µ 0 (Fig. 1).The unbounded background medium is truncated into a finite computation domain that encloses the dielectric object.Let Ω and Γ denote this computation domain and its boundary.In Ω, the permittivity is given by ε 0 ϵ r (r) and the permittivity is given by µ 0 µ r (r).Note that ϵ r (r) = 1 and µ r (r) = 1 in the background medium enclosed in Ω and ϵ r (r) ̸ = 1 and µ r (r) ̸ = 1 inside the scatterer.The speed of light in the background medium is given by c 0 = 1/ √ ε 0 µ 0 .
The electric and magnetic fields incident on the object are represented by E inc (r) and H inc (r), respectively.It is assumed that the incident fields and all fields and currents generated as a result of this excitation are time-harmonic with time dependence e jωt , where t is the time and ω is the frequency of excitation.Let E sca (r) and H sca (r) denote the electric and magnetic fields scattered from the object, respectively.Then, one can express the total electric and magnetic fields as E(r) = E inc (r) + E sca (r) and H(r) = H inc (r) + H sca (r).On the computation domain boundary Γ, equivalent electric and magnetic currents are defined as J(r) = n0 (r) × H(r) and M(r) = −n 0 (r) × E(r).Here, n0 (r) is the outward-pointing unit normal vector on Γ.Note that the formulation presented in the rest of this section is derived for normalized electric fields E inc (r), E sca (r), and E(r), and the normalization factor is ϵ 0 /µ 0 .The wavenumber in the background medium is given by k The formulation presented in the rest of this section heavily uses two trace operators: (i) on S. Note that n(r) is the outward-pointing unit normal vector on S.
Furthermore, to keep the formulation concise, the inner products used by the Galerkin scheme are not written explicitly.The notation and the definition of the inner products between two vectors u(r) and v(r) in volume V and surface S are given by respectively.
In the rest of the formulation, the dependence of the variables and the operators on r is dropped for the sake of simplicity in the notation unless a new variable or an operator is introduced.

Computation Domain Discretization
The computation domain Ω is discretized into a mesh of non-overlapping tetrahedrons rep-  Ω h that have all their three corners on Γ. J and M on Γ are approximated by J h and M h that are expanded on pairs of Γ i of Γ h .
The traditional HDG method uses a global vector field, which is denoted by Λ h , to "connect" local solutions E h and H h on Ω i of Ω h [27,29].The HDG-BI solver proposed in this work uses Λ h to also "connect" E h and H h in Ω h to J h and M h on Γ h .This global unknown Λ h is defined on the "shared" triangular surfaces of Ω h and the triangular surfaces of Γ h : L h = S h Γ h , where S h = l S l , S l is the triangular surface shared by two tetrahedrons Ω i and Ω j , i.e., S l = ∂Ω i ∩ ∂Ω j (Fig. 2), and Γ h = i Γ i (as already defined above).

HDG-BI Formulation 2.3.1 HDG
In computation domain Ω, the electric field E and the magnetic field H satisfy the timeharmonic Maxwell's equations: Similar to the traditional DG schemes and FEM, HDG seeks E h and H h that are approximate solutions of (3) and ( 4) defined on mesh Ω h discretizing Ω.This is achieved via weak Galerkin formulation of (3) and (4).Let e and h represent the testing functions corresponding to E and H respectively.Then, in a given tetrahedron Ω i , one can express the weak form of ( 3) and ( 4) as Using the mathematical identity for the divergence of the cross product of two vectors on the second terms of the inner products and applying the divergence theorem to the resulting expressions [15], one can convert ( 5) and ( 6) into Here, E t h and H t h are the tangential components of E h and H h on ∂Ω i and are expressed as {H h }, respectively.To "couple" the local system of equations associated with tetrahedron Ω i in ( 7) and ( 8) to the global system equations, numerical fluxes Ĥt h and Êt h are introduced as [27,29]: where τ is a local stabilization parameter and it is set to 1.0 in the rest of formulation and the code that implements this formulation.Unlike the traditional DG schemes, where the local fields E h and H h of a given tetrahedron are coupled to those of its neighboring tetrahedrons via numerical fluxes that rely on mean and jump of the field values [1], the numerical fluxes used by HDG as described in ( 9) and ( 10) couple the local fields E h and H h and the global unknown Λ h .
Next, expressions of Êt h and Ĥt h in ( 9) and ( 10) are used to replace H t h and E t h in (7) and (8), respectively.This yields Note that the inner product (h, ∇ × E h ) Ω i in ( 12) is obtained after applying the divergence {h}, E t h ⟩ ∂Ω i and using the mathematical identity for the divergence of the cross product of two vectors on the resulting expression.
To ensure the continuity between local and global unknowns, one needs to enforce the field continuity condition on triangular surfaces of Ω h , namely S h = l S l and Γ h = j Γ j .
On a given shared/inner triangular surface S l = ∂Ω i ∩ ∂Ω j , the continuity of the fields in Ω i and Ω j is enforced using the numerical flux [28,29] π On a given boundary triangular surface Γ j = ∂Ω i ∩ Γ j , the continuity of the fields in Ω i across the computation domain boundary is enforced using the numerical flux Let η represent the testing function corresponding to the global unknown Λ h , Then, one can express the weak form of ( 13) and ( 14) as By collecting the weak forms for all tetrahedrons Ω i of Ω h and all triangular surfaces S l of S h and Γ j of Γ h , one can obtain the part of the matrix system that represents the HDG component of the HDG-BI solver [27,29].This matrix system and the hybrid method used to efficiently solve it are described in Section 2.4.

BI
The formulation of the governing equations for the BI component of the proposed solver starts with the well-known relationship between the scattered fields E sca and H sca and the equivalent currents J and M that are introduced on the computation domain boundary Γ [65]: Here, the integral operators L S {X}(r) and K S {X}(r) are given by where G 0 is theGreen's function of the unbounded medium with wavenumber k 0 .Note that Inserting ( 17) and ( 18) into the current-field relationships and J = n0 × (H inc + H sca ) on Γ yields the electric-field integral equation (EFIE) and the magnetic-field integral equation (MFIE), respectively [48,[65][66][67]: To obtain a better-conditioned matrix upon discretization, linear combinations αEFIE( 19)+ 20) are used to yield the electric current combined-field integral equation (JCFIE) and the magnetic current combined-field integral equation (MCFIE), respectively [48,66,67]: Here, the combined-field integral operator C α S {X}(r) is defined as [67] and α is the weight that should be selected as 0 ≤ α ≤ 1.In the rest of formulation and the code that implements this formulation, α = 0.5.Accordingly C α S {X} is simplified as C S {X}.
JCFIE (21) and MCFIE (22) are approximated on Γ h that discretizes Γ: The continuity of the fields on Γ h is enforced using the numerical flux The governing BI equations are obtained via two linear combinations: JCFIE(24) + 1 2 (26) and MCFIE(25) + 1 2 π × Γ h {(26)}.This yields Let j and m represent the testing functions corresponding to J and M, respectively.j and m are the well-known Rao-Wilton-Glisson (RWG) basis functions [68] that are defined on pairs of Γ i .Let each of these pairs be represented by T j .Then, one can express the weak forms of ( 27) and (28) as By collecting the weak forms for all pairs of triangular surfaces T j of Γ h , one can obtain the part of the matrix system that represents the BI component of the HDG-BI solver.This matrix system and the hybrid method used to efficiently solve it are described in Section 2.4.

Matrix System
On Ω h , the local unknowns E h and H h are expanded as where e and h are the 3D zeroth-order vector edge basis functions [15] and Ē and H are the vectors storing the corresponding expansion coefficients.Similarly, on S h and Γ h , the global unknown Λ h is expanded as where η is the 2D zeroth-order vector edge basis function [15] and ΛS and ΛΓ are the vectors storing the coefficients of the expansions on S h and Γ h , respectively.
On Γ h , J h and M h are expanded as where j and m are the well-known RWG basis functions [68] as mentioned earlier.
Inserting the expansions in ( 31)-( 35) into the weak forms ( 11), ( 12), (15), and ( 16), and collecting the resulting equations for all tetrahedrons Ω i of Ω h and all triangular surfaces S l of S h and Γ j of Γ h , one can obtain the part of the matrix system that represents the HDG component of the HDG-BI solver: Here, the entries of the matrices Ā, F , B, L, DΛJ , and DΛM are given by To decrease the computational cost, the HDG scheme relies on reducing the size of the matrix system it solves.This is done by inverting (36) for Ē H T and inserting the resulting expression into (37).This operation yields: Here, the dimension of the matrix L − B Ā−1 F is equal to the number of degrees of freedom in the expansion of the global unknown Λ h as done in (33).Let this number be represented by N HDG , and let the total number of degrees of freedom in the expansions of E h and H h as done in (31) and (32) be represented by N DG .Since N HDG < N DG , the computational cost of the HDG scheme is significantly smaller than that of the traditional DG schemes [27,29].
Inserting the expansions in ( 33)-( 35) into the weak forms ( 29) and ( 30) and collecting the resulting equations for all pairs of triangular surfaces T m of Γ h yield the part of the matrix system that represents the BI component of the HDG-BI solver.Combining this system with the one for the HDG component in (44) yields the final matrix system of the HDG-BI solver as Here, the entries of the matrices DJΛ , DMΛ , CJJ , CJM , CMJ , and CMM and the entries of the right-hand side vectors bJ and bM are given by CJM mn = j m , The dimension of the matrix system ( 45) is N HDG + N BI , where N HDG is already defined above and N BI is the total number of degrees of freedom in the expansions of J h and M h as done in (34) and ( 35).This matrix system is solved for the unknown vectors ΛS , ΛΓ , J, and M using the scheme described in Section 2.5.

Hybrid Solver
The matrix system (45) can be re-written in a more compact form as: where matrices Q, DΛX , DXΛ , and C represent the four blocks of the matrix in (45)  Ideally, the matrix system (54) could be solved using a Krylov subspace-based iterative method assuming that the matrix-vector product associated with C is accelerated using MLFMA [55][56][57][58][59][60][61][62][63].However, Q is not well-conditioned [25] and as a result this iterative solution converges very slowly.Indeed, this is the reason why the traditional HDG schemes (in frequency domain) almost always rely on a direct (but sparse) matrix solver [28,29].
However, for the HDG-BI scheme developed in this work, using only a direct solver on (54) would be computationally expensive since C, which represents the BI component, is a full matrix.
To this end, in this work, a "hybrid" scheme is developed to efficiently solve the matrix system (54).The first row of ( 54) is inverted for Λ and the resulting expression is inserted into the second row to yield: This "reduced" matrix system of dimension N BI is solved using the hybrid scheme as described next step by step: 1. Apply LU decomposition to the sparse matrix Q as Q = L Ū and store matrices L and Ū .
2. Start the iterations of a Krylov subspace-based iterative scheme to solve (55).The matrix-vector multiplication subroutine of this iterative scheme is implemented as described by steps (a)-(c) below.Let x0 be the guess vector of this matrix-vector multiplication.
3. Continue the iterations until the relative residual error reaches the desired level.
The iterative solver used above is the general minimal residual method (GMRES) [69].The

LU decomposition in
Step (1) and the backward and forward substitutions in Step (2c) are carried out using the sparse matrix direct solver PARDISO [70].The efficiency and the accuracy of this hybrid solver are demonstrated by the numerical experiments described in Section 3.

Comments
Several comments about the formulation of the proposed HDG-BI solver detailed in Sections 2.1, 2.2, 2.3, 2.4, and 2.5 are in order: 1.The proposed solver allows for the surface of the dielectric object (which is determined by ε r and µ r ) to fully overlap with the computation domain boundary Γ.In such cases, the formulation detailed above stays the same without requiring any modifications.
This type of flexibility is especially important for a concave scatterer since enforcing the BI equations on this type of scatterer's surface significantly reduces the size of the computation domain (compared to using ABCs or PML) [53,54].
2. The proposed solver can be easily modified to efficiently account for disconnected scatterers.In this case, the BI equations should be enforced separately on the surface of each scatterer.This approach eliminates the need to discretize the space around the scatterers resulting in a very efficient solver especially for scenarios where the scatterers are well separated (compared to using ABCs or PML which would require a computation domain that encloses all scatterers and call for its full discretization) [53,54].

Perfect electrically conducting (PEC) objects possibly present in the computation do-
main Ω can easily be accounted for with a small modification of the numerical flux described by (13).Assume that the triangular surface S l ∈ ∂Ω i has all its three corners on the PEC surface, then ( 13) should be updated as [28,29] π 4. The formulation described in this section assumes a conformal mesh, i.e., the triangular surfaces of any two neighboring tetrahedrons match (share the same three nodes) and similarly the triangular surfaces of the tetrahedrons match to those discretizing the computation domain boundary.The expansions of the local HDG unknowns E h and H h , the global HDG unknown Λ h , and the BI unknowns J h and M h are carried out independently and "connected" to each other using the numerical flux.Theoretically, this approach allows for non-conformal meshes to be used to discretize each of these sets of unknowns.Such an approach would require defining/computing the numerical flux on overlapping regions of non-matching triangular faces of different mesh sets.To the best of authors' knowledge, an HDG method, which can account for non-conformal meshes, has been developed for only 2D problems [71,72].The possibility of extending this method to 3D problems and to account for non-conformal surface meshes (for incorporation of BI) will be investigated in a future publication.
5. Electromagnetic scattering problems are often analyzed using integral equation solvers (for example see [73][74][75]).Surface integral (SIE) solvers [73,74] when accelerated using MLFMA [55][56][57][58][59][60][61][62][63] result in the most computationally efficient methods for scattering analysis.However, their applicability is limited to problems where the material properties are piecewise homogenous.In problems where the scatterer is inhomogeneous, one can switch to a volume integral equation (VIE) solver [75], but this type of solvers requires a volumetric discretization of the scatterer.The HDG-BI solver can set the computation domain boundary Γ on the surface of the scatterer, ensuring that only the scatterer is discretized using a volumetric mesh.Under this condition, one can expect that the HDG-BI solver would be more efficient than the VIE solver.This is fundamentally because the volumetric discretization by the HDG-BI solver results in a sparse matrix while the volumetric discretization by the VIE solver results in a dense matrix.Indeed, the numerical results provided in Section 3.5 show that the proposed HDG-BI solver is faster that than a volume-surface integral equation (VSIE) solver in a problem where the scatterer is a PEC object embedded in a layered dielectric cube.

Numerical Results
In this section, several numerical examples are presented to demonstrate the accuracy, efficiency, and applicability of the proposed HDG-BI solver.In all examples, the scatterers are non-magnetic (µ r = 1 in the whole computation domain) and the background medium is the free space with permittivity ε 0 and permeability µ 0 .In all simulations, the excitation is a plane wave with electric and magnetic fields where E 0 = 1 V/m is the electric field amplitude, p is the unit vector along the direction of the electric field, k is the unit vector along the direction of propagation, and k 0 = 2πf /c 0 , η 0 = µ 0 /ϵ 0 , and c 0 are the wave number, impedance, and speed in the background medium, respectively.Here, f is the frequency of excitation, and the wavelength in the background medium at this frequency is given by λ 0 = f /c 0 .

Dielectric Coated PEC Sphere
In the first example, electromagnetic scattering from a dielectric coated PEC sphere is analyzed.The radius of the sphere is 0.3 m and the thickness of the coating is 0.1 m.The boundary of the computation domain (as denoted by Γ) is the outer surface of the coating.
The excitation parameters are f = 0.3 GHz, p = x, and k = ẑ.A total of 12 simulations are carried out using the HDG-BI solver for three different values of the coating's relative permittivity as 2.0, 4.0, and 8.0 and four different discretizations of the computation domain with average edge lengths 0.1λ 0 , 0.075λ 0 , 0.05λ 0 , and 0.03λ 0 .Table 1 provides the values of N HDG and N BI (as used by the HDG-BI solver) and N DG (as a reference) for these four different levels of discretization.In all simulations, the iterations of the GMRES method used in solving the matrix system (55) are terminated when the relative residual error reaches 0.001.At the end of each simulation, the L 2 -norm of the relative error in radar cross section (RCS) is computed using Here, σ is the RCS computed using J and M obtained by the HDG-BI on Γ, σ ref is the reference RCS computed using the Mie series solution, N = 180, ∆θ = 1.0 • , and ϕ = 0. Fig. 3 plots error σ versus average element length for three different values of the coating's relative permittivity.As expected, the error decreases with the increasing mesh density regardless of the value of the coating's relative permittivity.

Dielectric Plate
In the second example, electromagnetic scattering from a dielectric plate is analyzed.The  55) is solved using the GMRES method without using a preconditioner.(ii) The matrix system (55) is again solved using the GMRES method but this time a sparse approximate inverse (SAI) preconditioner [76,77] is used.This preconditioner is constructed using only ¯C.In both cases, the iterations of the GMRES method are terminated when the relative residual error reaches 0.001.
Fig. 4(b) plots the number of GMRES iterations versus N BI for these two cases.For both cases, the slope of the iteration number curve flattens with increasing N BI , which means that the problem size does not have much effect on the efficacy of the iterative solver.Also, the figure shows that the SAI preconditioner can reduce the number of iterations.But it is worth mentioning here that for a more complicated scatterer (complex shape, inhomogeneous permittivity, etc.), this type of preconditioning may not be as effective [77].fact that the computation domain boundary where the ABC is enforced has to be located away from the randome surface to achieve the same accuracy level as the HDG-BI solver.This increases the computation domain size and, accordingly the computational requirements of the HDG-ABC solver.

Aircraft Head
In this example, electromagnetic scattering from a coated aircraft head model is analyzed.
The aircraft head's cross section on the xz-plane is shown in Fig. 6.The relative permittivity of the coating is 2.0 − 0.5j.The excitation parameters are f ∈ {1.8, 3.6} GHz, k = −ẑ, and p ∈ {x, y}.Three sets of simulations are carried out.(i) HDG-BI: The boundary of the computation domain (as denoted by Γ) is the outer surface of the coating, i.e., HGD-BI discretizes only the volume of the coating and its surface.Two levels of discretization with average edge lengths 0.07λ 0 (resulting in N HDG = 813 357 and N BI = 81 192) and 0.075λ 0 (resulting in N HDG = 8 696 010 and N BI = 316 086) are used for the simulations with f = 1.8 GHz and f = 3.6 GHz, respectively.The iterations of the GMRES method used in solving the matrix system (55) are terminated when the relative residual error reaches 0.001.(ii) HDG-ABC: The computation domain is a sphere of radius 0.8 m.This sphere fully encloses the coated aircraft head and the first-order ABC is enforced on its surface.
The computation domain is discretized using elements with an average edge length of 0.07λ 0 resulting in N HDG = 8 505 207 for the simulation with f = 1.8 GHz.The HDG matrix system is solved using the sparse LU solver PARDISO [70].Note that for the simulation with f = 3.6 GHz, the HDG-ABC solver is not used because of its prohibitive computational requirements.(iii) MoM: The multi-trace surface integral equation solver described in [78,79] is used.Furthermore, the computational requirements of the HDG-BI solver are significantly lower than those of the HDG-ABC solver: The HDG-BI solver requires 6.7 GB of memory and completes the simulation in 165 s while the HDB-ABC solver requires 210.1 GB of memory and completes the simulation in 3.58 h.Note that the number of GMRES iterations required by the HDG-BI solver is only 45.The large difference in the computational requirements of these two solvers can be explained by the fact that the computation domain of the HDG-ABC and the degrees of freedom required its discretization (as represented by N HDG ) are significantly larger than those of the HDG-BI solver.

PEC Cylinder Embedded in a Layered Dielectric Cube
In the last example, electromagnetic scattering from a PEC cylinder embedded in a layered dielectric cube is analyzed.The geometry of the scatterer is shown in Fig. 8  The iterations of the GMRES method used in solving the matrix system (55) are terminated when the relative residual error reaches 0.001.(ii) VSIE: The commercially available software package FEKO [80] that solves a coupled system of VIE (enforced inside the cube) and SIE (enforced on the surface of the cylinder) is used.The average edge length in the software is set to 0.05λ 0 which results in 74 162 degrees of freedom for the VSIE solver.Note that neither the HDG-BI solver nor the VSIE solver is accelerated using MLFMA.(18.12 m to compute the matrix and 59.6 m to solve the matrix system).This comparison shows the benefits of the proposed HDG-BI solver over the VSIE solver, which mainly stems from the fact that the volumetric discretization by HDG results in a sparse matrix while the volumetric discretization by VSIE results in a dense matrix.

Conclusions
A method, which couples the HDG and the BI equations is developed to efficiently analyze electromagnetic scattering from inhomogeneous/composite objects.The coupling between these two sets of equations is realized using the numerical flux operating on the equivalent current and the global unknown of the HDG.This approach yields sparse coupling matrices upon discretization.Inclusion of the BI equation ensures that the only error in enforcing the radiation conditions is the discretization.Furthermore, the computation domain boundary, where the BI equation is enforced, can be located very close, even conformal, to the surface of the scatterer without any loss of accuracy.This significantly reduces the number of unknowns to be solved for compared to the traditional HDG schemes that make use of ABCs or PML to truncate the computation domain.
However, the discretization of the BI equation yields a dense matrix, which prohibits the use of a direct matrix solver on the overall coupled system as often done with traditional HDG schemes.To overcome this bottleneck, a "hybrid" method is developed.This method uses an iterative scheme to solve the overall coupled system but within the matrix-vector multiplication subroutine of the iterations, the inverse of the HDG matrix is efficiently accounted for using a sparse direct matrix solver.The same subroutine also uses the multilevel fast multipole algorithm to accelerate the multiplication of the guess vector with the dense BI matrix.Numerical examples show that the proposed HDG-BI solver has clear advantages over the traditional HDG schemes with ABCs and a VSIE solver.
As future work, a domain decomposition method and high-order vector basis functions will be incorporated into the HDG-BI solver to further improve its efficiency and accuracy and its applicability to large-scale problems.Additionally, a discretization scheme that can account for non-conformal meshes will be formulated and implemented within the framework of the HDG-BI solver.

Figure 1 :
Figure 1: Description of the electromagnetic scattering problem.
yields the tangential components of u(r) on surface S, and (ii) π × S {u}(r) = n(r) × u(r)| S that yields the twisted tangential components of u(r)

Figure 2 :
Figure 2: Description of the mesh supporting the local unknowns E h and H h and the global unknown Λ h .

Figure 3 :
Figure 3: L 2 -norm of the relative error in RCS error σ [computed using (58)] for different values of ε r versus average edge length of discretization.
and the vectors Λ, X, and b represent the corresponding parts of the unknown and the right-hand side vectors.L, B, and F are sparse matrices and Ā and Ā−1 are block-diagonal matrices.Therefore, Q is a sparse matrix.DΛX and DXΛ are also sparse matrices since their blocks DJΛ , DMΛ , DΛJ , and DΛM are all sparse matrices.C is a full matrix since its blocks CJJ , CJM , CMJ , and CMJ are all full matrices.

Figure 4 :
Figure 4: Electromagnetic scattering from a dielectric plate.(a) Description of the geometry and the excitation.(b) Number of iterations required by the GMRES method (without a preconditioner and with SAI preconditioner) versus N BI .
dimensions of the plate are L × L × h as shown in Fig. 4(a), and its relative dielectric permittivity is 2.0.The boundary of the computation domain (as denoted by Γ) is the surface of the plate.Four simulations are carried out for four different values of L, L ∈ {3.0, 6.0, 12.0, 24.0} m.In all simulations, h = 0.1 m and the excitation parameters are f = 0.3 GHz, p = x, k = −ẑ.The average edge lengths in the discretizations on the surface and in the volume of the plate are 0.1λ 0 and 0.075λ 0 , respectively, resulting in N HDG ∈ {16 468, 65 168, 259 965, 1 037 761} and N BI ∈ {4 174, 15 550, 59 870, 234 990} for four values of L. Two cases are considered for each simulation: (i) The matrix system (

Figure 5 :
Figure 5: Electromagnetic scattering from a dielectric radome.(a) Description of the geometry (cross section) and the excitation.(b) Real and (c) imaginary part of J and M computed by HDG-BI and MoM on the surface of the radome.(d) RCS obtained using J and M that are computed by HDG-BI, HDG-ABC, and MoM.

Figure 6 :
Figure 6: Electromagnetic scattering from a coated aircraft head.(a) Description of the geometry (cross section) and the excitation.(b) Real and (c) imaginary part of J and M computed by HDG-BI and MoM on the outer surface of the dielectric coating at for p = x at f = 1.8 GHz.RCS obtained using J and M that are computed by HDG-BI and MoM for (d) p = x and (e) p = ŷ at f = 1.8 GHz.

Figure 7 :
Figure 7: Electromagnetic scattering from an aircraft head.(a) Real and (b) imaginary part of J and M computed by HDG-BI and MoM on the outer surface of the dielectric coating for p = x at f = 3.6 GHz.

Fig. 6 (
Fig. 6 (b) and (c) compare the real and the imaginary of parts of J and M obtained by HDG-BI and MoM solvers on the outer surface of the dielectric coating for the simulation with f = 1.8GHz and p = x.The results agree well.RCS is computed for θ ∈ [0 180 • ] and ϕ = 0 using J and M obtained by the HDG-BI, HDG-ABC, and MoM solvers in two simulations with p = x and p = ŷ.Fig. 6 (d) and (e) plots RCS computed by these three solvers versus θ in the simulations with p = x and p = ŷ, respectively.The figure shows that the results obtained by the HDG-BI solver agree very well with those obtained by the MoM solver while the results obtained by the HDG-ABC solver do not.The inaccuracy of the HDG-ABC solver can be explained by the fact that the first-order ABC is used to truncate the computation domain.As expected, the HDG-BI solver does not suffer from this bottleneck.
(a) and (b) compare the real and the imaginary of parts of J and M obtained by HDG-BI and MoM solvers on the outer surface of the dielectric coating for the simulation with f = 3.6 GHz and p = x.The results agree well.The HDG-BI solver completes this simulation in 1 993 s and requires 99.5 GB of memory.The number of GMRES iterations is only 35.
(a).The relative permittivities of the four layers (ordered from top to bottom) are 3.0, 2.0, 3.0, and 2.0.The excitation parameters are f = 0.3 GHz, p = x, and k = −ẑ.Two simulations are carried out: (i) HDG-BI: The boundary of the computation domain (as denoted by Γ) is the surface of the cylinder and the outer surface of the cube.The computation domain is discretized using elements with an average edge length of 0.05λ 0 resulting N HDG = 201 294 and N BI = 11 856.

Fig. 8 (Figure 8 :
Fig. 8 (b) plots RCS computed for θ ∈ [0 180 • ] and ϕ = 0 in these two simulations.Results agree well demonstrating the accuracy of the proposed HDG-BI solver.For this where Ω i is the ith tetrahedron.The boundary of tetrahedron Ω i , which consists of four triangular surfaces, is represented by ∂Ω i .E and H in Ω are approximated by E h and H h that are expanded on Ω i of Ω h .
h : Γ h = i Γ i , where Γ i are the triangular surfaces of

Table 1 :
N HDG , N BI , and N DG of the discretizations used in the electromagnetic simulation of scattering from the coated sphere.