A Fast-Converging Scheme for the Electromagnetic Scattering from a Thin Dielectric Disk

In this paper, the analysis of the electromagnetic scattering from a thin dielectric disk is formulated as two sets of one-dimensional integral equations in the vector Hankel transform domain by taking advantage of the revolution symmetry of the problem and by imposing the generalized boundary conditions on the disk surface. The problem is further simplified by means of Helmholtz decomposition, which allows to introduce new scalar unknows in the spectral domain. Galerkin method with complete sets of orthogonal eigenfunctions of the static parts of the integral operators, reconstructing the physical behavior of the fields, as expansion bases, is applied to discretize the integral equations. The obtained matrix equations are Fredholm second-kind equations whose coefficients are efficiently numerically evaluated by means of a suitable analytical technique. Numerical results and comparisons with the commercial software CST Microwave Studio are provided showing the accuracy and efficiency of the proposed technique.


Introduction
Many applications in the framework of frequency-selective surfaces, graphene structures, radome design, microstrip antennas and antenna arrays, just as examples, involve thin layers of dielectric/conducting materials [1][2][3][4][5][6][7][8][9][10]. For these reasons, such structures have been extensively studied in the past decades and even more recently. The literature devoted to this subject is principally focused on finding approximate solutions based on Rayleigh-Gans or physical optics techniques, or on proposing suitable numerical methods to achieve more accurate solutions. In the last case, when the thickness of the involved structure is much less than the free-space wavelength, the classical 3D approaches reveal to be not particularly suitable in terms of computation time, storage requirements and they can even suffer from lack of convergence.
On the other hand, the problem can be significantly simplified if the scatterer is approximated as a surface with the fields satisfying suitable boundary conditions, generally called generalized boundary conditions [11,12]. In such a case, surface integral equation formulations for the effective electric and/or magnetic current densities, taking into account the radiation conditions of the fields, are particularly attractive because the integral equations and the unknowns are defined on a finite support [12,13].
The discretization of the obtained equation is, however, a key point because the convergence of the approximate solution to the exact one cannot be generally established [14]. Moreover, it can lead to a huge ill-conditioned matrix to be numerically inverted. Sometimes, compression techniques [15,16] and/or suitable discretization schemes [17][18][19] can be usefully applied. However, a general way to overcome these problems is to recast the integral equation in hand as a Fredholm second-kind integral equation by analytically inverting a part of the operator containing the most singular one [20]. Then, due to the Fredholm's theory, every discretization scheme preserving the Fredholm's nature of the obtained second-kind integral equation is a convergent scheme [21].
A different approach consists in properly select the discretization scheme in order to obtain a guaranteed convergence matrix equation [20]. This means that regularization and discretization are condensed in a single step. The wide literature devoted to this subject, in the context of the analysis of propagation, radiation and scattering problems, have demonstrated that one way to reach this goal is the application of the Galerkin method with a suitable selection of the expansion functions such that the most singular part of the discretized operator is invertible with a continuous two-side inverse, and the remaining part is a compact operator [22][23][24][25][26][27][28][29][30][31][32]. Moreover, if the physical behavior of the unknowns is properly reconstructed by the expansion functions, the convergence is even fast, i.e., highly accurate results can be achieved by means of a small size coefficients' matrix, with an obvious advantage in terms of computation time and storage requirement.
In this paper, an analytically regularizing approach for the analysis of the electromagnetic scattering from a thin dielectric disk in a homogeneous medium is presented. Two independent surface integral equations for the effective electric and magnetic currents, respectively, are obtained by imposing the generalized boundary conditions on the disk surface. Due to the revolution symmetry of the problem, all the involved functions are conveniently expanded in Fourier series and the surface integral equations are recast as sets of one-dimensional integral equations in the spectral domain for the components of each harmonic of the vector Hankel transforms (VHTs) of the effective currents. Helmholtz decomposition allows to further simplify the problem by introducing new scalar unknowns to be discretized. Galerkin method with complete sets of orthogonal eigenfunctions of the static parts of the involved integral operators, reconstructing the behavior of the unknowns at the edge and around the center of the disk, as expansion bases is used to discretize the integral equations. In this way, the Galerkin projection acts as a perfect preconditioner and the obtained matrix equations are the Fredholm second-kind matrix operator equations. Moreover, the convolution integrals are reduced to algebraic products in the spectral domain. Hence, the resulting matrix coefficients are combinations of one-dimensional improper integrals efficiently evaluated by means of a suitable analytical procedure. The presented numerical results and the comparisons with the commercial software CST Microwave Studio (CST-MWS) clearly show the accuracy and efficiency of the proposed technique. This paper is organized as follows: Section 2 is devoted to the formulation and solution of the problem, Section 3 shows the numerical results, Section 4 is dedicated to the conclusions and an Appendix A concludes the paper.

Integral Equations in the Spectral Domain
In Figure 1, a thin dielectric disk of radius a, thickness τ, dielectric permittivity ε and magnetic permeability µ is immersed in free space of dielectric permittivity ε 0 and magnetic permeability µ 0 . A cylindrical coordinate system (ρ, φ, z) with the origin at the center of the disk and the z axis orthogonal to it is introduced.
are the effective electric and magnetic currents, respectively, and: are the electric and magnetic resistivities of the disk, respectively. It is interesting to observe that the integral Equations (1a) and (1b) are decoupled for the problem in hand [11]. The uniqueness of the solution is guaranteed by requiring that the fields components satisfy, additionally, the edge conditions and the radiation conditions [20].
The revolution symmetry of the problem allows to expand all the involved functions in Fourier series. Hence, the surface integral equations for the effective electric and magnetic currents can be equivalently reduced to two infinite sets of independent one-dimensional integral equations: Figure 1. Geometry of the problem: (a) a thin dielectric disk, of radius a, thickness τ, dielectric permittivity ε and magnetic permeability µ, immersed in free space of dielectric permittivity ε 0 and magnetic permeability µ 0 ; (b) a plane wave impinging onto the disk surface with propagation direction (θ 0 , φ 0 ).
A plane wave, i.e., E inc (r) = E 0 e −jk·r and H inc (r) = ε 0 /µ 0k × E 0 e − jk·r , where k = −k 0 (sin θ 0 cos φ 0x + sin θ 0 sin φ 0ŷ + cos θ 0ẑ ) is the propagation vector, k 0 = 2π/λ = ω √ ε 0 µ 0 is the free-space wavenumber, λ is the free-space wavelength, ω is the angular frequency, (θ 0 , φ 0 ) identifies the propagation direction and r = xx + yŷ + zẑ, impinges onto the disk surface generating a scattered electromagnetic field (E sc (r), H sc (r)). If εµ ε 0 µ 0 , τ a and τ λ, the disk can be approximated as a flat zero-thickness surface, on which the following generalized boundary conditions have to be satisfied [11,12,18] for ρ ≤ a, where E tot (r), H tot (r) denotes the total electromagnetic field: are the effective electric and magnetic currents, respectively, and: are the electric and magnetic resistivities of the disk, respectively. It is interesting to observe that the integral Equations (1a) and (1b) are decoupled for the problem in hand [11]. The uniqueness of the solution is guaranteed by requiring that the fields components satisfy, additionally, the edge conditions and the radiation conditions [20].
The revolution symmetry of the problem allows to expand all the involved functions in Fourier series. Hence, the surface integral equations for the effective electric and magnetic currents can be equivalently reduced to two infinite sets of independent one-dimensional integral equations: for ρ ≤ a, where the apex n denotes the n-th term of the Fourier series, the symbol: has been introduced, and the n-th harmonics of the scattered electric and magnetic fields, respectively, can be expressed in the spectral domain as follows [32]: where: J n (·) and J n (·) are the Bessel function of the first kind and order n and its first derivative with respect to the argument, respectively [33]: with r ∈ {e, m}, χ e = ε 0 , χ m = µ 0 , k 2 0 − w 2 = −j −k 2 0 + w 2 , and the following definition of the VHT of order n (VHT n ) has been introduced [34]: Observing that [34]: the Equations (4) can be immediately rewritten as follows: where I is the identity matrix operator.

Helmholtz Decomposition
Therefore, the problem in hand has been formulated in terms of two infinite sets of one-dimensional integral equations in the spectral domain for the VHT of the harmonics of the effective currents, which are vector functions. On the other hand, the problem can be further simplified by means of the Helmholtz decomposition which leads to scalar unknowns instead of vector ones.
Let P(ρ, φ) be the effective electric/magnetic current density. Due to the Helmholtz decomposition, it can be written as the superposition of a surface curl-free contribution D} are suitable potential functions [35]. Moreover, the n-th harmonic of P T (ρ, φ) with T = C or T = D can be immediately written as a function of the n-th harmonic of the corresponding potential function, i.e., It is possible to demonstrate (see Appendix A) that the VHT n of the function P where: Henceforth, the functions (13) are assumed as new unknowns in the spectral domain.

Discretization of the Integral Equations
The obtained integral equations are discretized by means of the Galerkin method. Following the line of reasoning in [32], it is possible to show that the behavior of the n-th harmonic of the effective currents at the edge and around the center of the disk is correctly reconstructed by expanding the unknown functions ∼ P (n) T (w) in complete series of weighted Bessel functions [32,36], i.e., where δ n,m is the Kronecker delta, γ (n) T,h denote the expansion coefficients, η (n) C,−1 for n 0 can be established by imposing the currents to be vanishing outside the disk surface.
Due to the spectral domain formulation, the convolution integrals are reduced to algebraic products. Hence, the coefficients of the obtained matrix equations can be readily expressed as combinations of one-dimensional integrals of the kind: where: By extracting the asymptotic behavior of the kernel in (16), i.e., w , and using the well-known Weber-Schafheitlin discontinuous integral [37], the following alternative expression for the integral in (16) can be obtained: Since it is possible to demonstrate that the matrix operators associated to the first term at the right hand side of (18) are compact operators in l 2 [32], the obtained equations are the Fredholm second-kind matrix operator equations in l 2 [21], and the convergence of the discretization scheme is guaranteed.
The numerical convergence of improper integrals of oscillating functions with an algebraic decay like the one in (18) can be time-consuming when highly accurate results are searched for. This problem can be completely overcome by adopting the analytical procedure in the complex plane developed in [38], which leads to the numerical evaluation of fast converging proper integrals.

Results and Discussion
An approximate solution can be obtained by truncating the obtained infinite matrix equations. In order to show the fast convergence of the presented method, the following normalized truncation error is introduced: where 2N − 1 is the number of considered harmonics estimated as in [39], · is the usual Euclidean norm and x (n) M is the vector of all the expansion coefficients evaluated by using M expansion functions for each unknown. All the simulations are performed on a laptop equipped with an Intel Core i7-8550u 1.8 GHz CPU and 16 GB RAM, running Windows 10 and the integrals of the coefficients' matrix evaluated by means of an adaptive Gauss-Legendre quadrature routine implemented in the Matlab environment.
In Figure 2, the normalized truncation error, when a plane wave with incidence direction θ 0 = 45 • , φ 0 = 0 • , electric field amplitude E 0 = 1 V/m, TE incidence and TM incidence impinges onto a disk with a normalized radius k 0 a = 1, 3, 5, thickness τ = 0.05a, dielectric permittivity ε = (10.5 − j0.3)ε 0 and magnetic permeability µ = µ 0 , is plotted as a function of M for N = 5, 9, 11, respectively. θ =°, 0 0 φ =°, electric field amplitude 0 1V /m E = , TE incidence and TM incidence impinges onto a disk with a normalized radius As can be seen in Figure 2, the convergence is very fast in all the examined cases. Indeed, in the worst case, an error less than 2 10 − is achieved form M = 10, while M = 18 allows to obtain an error less than 3 10 − . Moreover, an increased normalized radius leads only to a greater number of harmonics to be used. This is because the values of the electric and magnetic resistivities are functions of the normalized thickness ( 0 k τ ) and, therefore, a change in the normalized radius can increase or As can be seen in Figure 2, the convergence is very fast in all the examined cases. Indeed, in the worst case, an error less than 10 −2 is achieved form M = 10, while M = 18 allows to obtain an error less than 10 −3 . Moreover, an increased normalized radius leads only to a greater number of harmonics to be used. This is because the values of the electric and magnetic resistivities are functions of the normalized thickness (k 0 τ) and, therefore, a change in the normalized radius can increase or decrease the rate of convergence. In Figure 3, the normalized truncation error is plotted as a function of M for three incidence directions, i.e., for θ 0 = 0 • , 45 • , 90 • and φ 0 = 0 • , and for E 0 = 1 V/m, TE incidence and TM incidence, N = 2, 9, 11, k 0 a = 3, ε = (10.5 − j0.3)ε 0 , µ = µ 0 . As expected, a change in the incident angle θ 0 affects the number of harmonics to be used in order to achieve a given accuracy. On the other hand, for TE incidence, the change of θ 0 does not significantly affect the number of expansion functions to be used while, for TM incidence, the convergence is even faster if θ 0 is larger. To conclude, it is worth noting that a normalized truncation error less than 10 −3 can be achieved in the worst of the cases examined in this paper by filling a coefficient matrix of about 50,000 elements, with a computation time less than 2 min.  As clearly stated before, the convergence of the proposed method is guaranteed. Hence, comparisons with CST-MWS are just provided in order to validate the numerical implementation and to show the accuracy and the efficiency of the proposed approach.
In Figure 4   As clearly stated before, the convergence of the proposed method is guaranteed. Hence, comparisons with CST-MWS are just provided in order to validate the numerical implementation and to show the accuracy and the efficiency of the proposed approach.
In Figure 4 the Bistatic Radar Cross Section (BRCS) for a disk with k 0 a = 3, τ = 0.05a, ε = (10.5 − j0.3)ε 0 and µ = µ 0 , on which a plane wave orthogonally impinges with E 0 = E 0ŷ , obtained for M = 19 and N = 2 in order to achieve a normalized truncation error less than 10 −3 , is plotted and compared with the results provided by CST-MWS. The agreement is quite good. It is interesting to note, however, that CST-MWS requires a computation time of 80 min and 9.4 M mesh-cells to reconstruct the plotted solution. After all, to the best of the authors knowledge, CST-MWS does not provide a 2D model for dielectric objects and an accurate simulation of the considered 3D object has turned out to be time-consuming and particularly burdensome in terms of memory requirement. As a result, the proposed approach drastically outperforms CST-MWS in terms of both computation time and storage requirements.   The proposed method allows the simple reconstruction of near field and far field parameters. For the sake of completeness, in Figures 5-7, the effective electric and magnetic currents and the BRCSs for a disk with k 0 a = 3, τ = 0.05a, ε = (10.5 − j0.3)ε 0 and µ = µ 0 , when a plane wave with E 0 = 1 V/m, TE incidence and TM incidence, and three incidence directions, i.e., θ 0 = 0 • , 45 • , 90 • and φ 0 = 0 • , impinges onto the disk surface, are shown.    In Figure 5, where the TE incidence is considered, the absolute values of the non-vanishing components of the effective electric and magnetic currents along the x axis (  In Figure 5, where the TE incidence is considered, the absolute values of the non-vanishing components of the effective electric and magnetic currents along the x axis (  Figure 6, the absolute values of the non-vanishing components of the effective currents along the x axis obtained for TM incidence, i.e., the radial component of the effective electric current and the azimuthal component of the effective magnetic current, are plotted as a function of a ρ . As expected, the radial component of the effective magnetic current for TE incidence and the radial component of the effective electric current for TM incidence vanish at grazing incidence ( 0 90 θ =°) due to the symmetry of the problem with respect to the z axis. To conclude, in Figure 7, the BRCSs for both TM incidence and TE incidence are plotted in the plane 0 , 180 φ =°° as a function of the observation angle θ and compared with quite good agreement with the results provided by CST-MWS by means of the same discretization parameters used to reconstruct the behavior in Figure 3. It is worth noting that the low number of expansion functions and harmonics used to reconstruct the currents and the BRCSs with the proposed approach, i.e., M = 19 and N = 11, is, however, enough to guarantee a normalized truncation error less than 10 −3 . In Figure 5, where the TE incidence is considered, the absolute values of the non-vanishing components of the effective electric and magnetic currents along the x axis (φ = 0 • , 180 • ), i.e., the azimuthal component of the effective electric current and the radial component of the effective magnetic current, are plotted as functions of ρ/a. Similarly, in Figure 6, the absolute values of the non-vanishing components of the effective currents along the x axis obtained for TM incidence, i.e., the radial component of the effective electric current and the azimuthal component of the effective magnetic current, are plotted as a function of ρ/a. As expected, the radial component of the effective magnetic current for TE incidence and the radial component of the effective electric current for TM incidence vanish at grazing incidence (θ 0 = 90 • ) due to the symmetry of the problem with respect to the z axis. To conclude, in Figure 7, the BRCSs for both TM incidence and TE incidence are plotted in the plane φ = 0 • , 180 • as a function of the observation angle θ and compared with quite good agreement with the results provided by CST-MWS by means of the same discretization parameters used to reconstruct the behavior in Figure 3. It is worth noting that the low number of expansion functions and harmonics used to reconstruct the currents and the BRCSs with the proposed approach, i.e., M = 19 and N = 11, is, however, enough to guarantee a normalized truncation error less than 10 −3 .

Conclusions
In this paper, a new analytically regularizing scheme for the analysis of the electromagnetic scattering from a thin dielectric disk has been presented. The method is accurate and efficient, leads to a quick evaluation of near-field and far-field parameters and drastically outperforms the general-purpose commercial software CST-MWS. Future perspectives are the generalization of the presented method to the analysis of arrays of disks and disks in layered media.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The aim of this appendix is to demonstrate Equations (13) and (14). Starting from the definition of VHT n , the following relation can be immediately stated: On the other hand, the following expression for the first component of the vector function (A2) can be obtained by integrating by parts and using the Bessel's equation [33]: Hence, Equations (13a) and (14) for T = C have been demonstrated. Following the same line of reasoning, it is possible to justify even Equations (13b) and (14) for T = D.