Analysis of the Scattering from a Two Stacked Thin Resistive Disks Resonator by Means of the Helmholtz–Galerkin Regularizing Technique

: In this paper, the scattering of a plane wave from a lossy Fabry–Per ó t resonator, realized with two equiaxial thin resistive disks with the same radius, is analyzed by means of the generalization of the Helmholtz–Galerkin regularizing technique recently developed by the author. The disks are modelled as 2-D planar surfaces described in terms of generalized boundary conditions. Taking advantage of the revolution symmetry, the problem is equivalently formulated as a set of independent systems of 1-D equations in the vector Hankel transform domain for the cylindrical harmonics of the effective surface current densities. The Helmholtz decomposition of the unknowns, combined with a suitable choice of the expansion functions in a Galerkin scheme, lead to a fast-converging Fredholm second-kind matrix operator equation. Moreover, an analytical technique speciﬁcally devised to efﬁciently evaluate the integrals of the coefﬁcient matrix is adopted. As shown in the numerical results section, near-ﬁeld and far-ﬁeld parameters are accurately and efﬁciently reconstructed even at the resonance frequencies of the natural modes, which are searched for the peaks of the total scattering cross-section and the absorption cross-section. Moreover, the proposed method drastically outperforms the general-purpose commercial software CST Microwave Studio in terms of both CPU time and memory occupation.


Introduction
The great interest in studying electromagnetic scattering from thin objects originates from the huge number of applications in which they are involved: frequency selective surfaces, antennas, radomes, tree leaves models, etc. More recently, the discovery of graphene [1], which is a planar monolayer of carbon atoms with amazing electronic, optical, thermal, and mechanical properties, and the consequent countless number of devices imagined to be made with such a material, have been leading to an impressive push towards the accurate study of lossy thin structures.
Boundary and radiation conditions together with the local power boundedness condition allow to formulate a uniquely solvable boundary value problem for the Maxwell equations, for which, in general, no closed form solutions are available, except for problems involving suitable canonical structures [2]. In other cases, only approximate solutions are provided based on Rayleigh-Gans or physical optics [3][4][5]. Consequently, the analysis of general structures or the search for more accurate solutions require the use of numerical methods. Amongst them, the classical Moments' method approach is generally preferred because the unknowns and the integral equations are defined on finite supports, the radiation condition is automatically satisfied by the proper selection of the Green's function of the problem and the local power boundedness condition is guaranteed by the correct definition of the functional space at which the unknowns belong.
Unfortunately, the discretization of a 3-D problem is generally onerous in terms of storage requirements. On the other hand, when the thickness is much smaller than the wavelength, the scatterer can be modelled as an infinitesimally thin structure described in terms of suitable generalized boundary conditions [6]. Just for an example, this is what happens when dealing with a graphene object, for which the surface impedance can be obtained from the graphene surface conductivity provided by the Kubo formalism [7]. According to this approach, the boundary value problem for the Maxwell equations at hand can be equivalently formulated in terms of a uniquely solvable system of 2-D singular integral equations for the effective surface current densities [8,9]. It is worth noting that, due to the singular nature of the obtained integral equations, nothing can be established a priori about the existence of the solution. Moreover, the convergence of the approximate solution, obtained by adopting a discretization scheme, to the exact one, if it exists, cannot be, in general, predicted [10]. Anyway, ill-conditioned and/or dense matrices can result from a direct discretization.
One way to overcome all these problems is represented by the methods of analytical regularization [11]. Such methods are aimed at individuating a suitable singular part of the integral operator, containing the most singular part of the operator itself, to be analytically inverted. The selection of such an operator depends on the problem at hand. Usually, it can be a suitable asymptotic part (e.g., the static part or the high-frequency part, etc.) or a canonical-shape part, which can be inverted by means of functional techniques such as Wiener-Hopf, Cauchy, Titchmarsh, Abel, Riemann-Hilbert Problem, and the separation of variables techniques [12][13][14][15][16][17][18][19][20]. Following this line of reasoning, the resulting integral equation is of the Fredholm second-kind; hence, the Fredholm theory [21], generalized by Steinberg for operators [22], can be applied. As a result, the existence originates from the uniqueness and the convergence of a discretization scheme preserving the Fredholm second-kind nature of the obtained integral equation is guaranteed. On the other hand, a proper choice of the discretization scheme can lead directly to a Fredholm second-kind matrix equation. This is what happens when: (1) the Galerkin scheme is adopted, and (2) the selected expansion functions are orthonormal eigenfunctions of a suitable operator containing the most singular part of the integral operator at hand. Such an approach, appropriately called method of analytical preconditioning, is very effective, as clearly shown in the literature devoted to the study of the scattering, propagation, and radiation problems [23][24][25][26][27][28][29][30][31][32][33]. Another way to obtain a guaranteed-convergence consists in solving numerically the singular integral equation by means of a Nyström-type discretization scheme taking into account the singularity of the integral equation and the behavior of the unknowns at the edges [34][35][36][37][38].
When dealing with thin disks, the complexity of the problem can be significantly reduced [39][40][41][42][43]. Indeed, the Fourier series expansion of the unknowns and the orthogonality property of the cylindrical harmonics allow recasting the surface integral equations as infinite sets of independent 1-D integral equations for the harmonics of the effective surface current densities. The strategy of combining the vector Hankel transform (VHT) [44] and the Helmholtz decomposition [45] is particularly suitable because the spectral domain counterpart of the surface curl-free and divergence-free contributions of the currents are scalar functions. After assuming such functions as new unknowns in the spectral domain, the method of analytical preconditioning can be readily applied. The selected expansion functions in the spectral domain, diagonalizing the most singular part of the integral operator, have a closed-form spatial domain counterpart reconstructing the behavior of the unknowns at the disk rim and around the centre of the disk. As a result, the projection integrals reduce to algebraic products and the convergence is very fast. Moreover, the 1-D improper integrals of the coefficient matrix can be quickly and accurately evaluated by adopting the analytical technique developed in [46] and generalized in [40].
In this paper, the method described above is successfully generalized to the analysis of a lossy Fabry-Perót resonator realized with two equiaxial thin resistive disks with the same radius. As will be clearly shown: (1) The proposed method is fast convergent and allows Appl. Sci. 2021, 11, 8173 3 of 16 one to accurately and efficiently reconstruct both near-field and far-field parameters; (2) The resonance frequencies of the structure are simply individuated by the peaks of the total scattering cross-section (TSCS) and the absorption cross-section (ACS), which are expressed in closed form; (3) The proposed method drastically outperforms the general-purpose commercial software CST Microwave Studio (CST-MWS) in terms of both CPU time and memory occupation.
This paper is organized as follows. Section 2 is devoted to the description of the proposed method. In Section 3, the obtained numerical results are presented, and the conclusions are summarized in Section 4.

Formulation of the Problem and Proposed Solution
Two equiaxial resistive disks with the same radius a, thickness τ i and electric conductivity σ i , where i = 1, 2 identifies the i-th disk, in free space are considered (see Figure 1). A Cartesian coordinate system, (x, y, z), and a cylindrical coordinate system, (ρ, φ, z), are introduced so that the z axis coincides with the axis of the disks and the median surfaces of the disks are located at the abscissas z = z i . Henceforth, the symbols ε 0 , µ 0 , f , ω = 2π f , λ, k 0 = 2π/λ = ω √ ε 0 µ 0 and Z 0 = µ 0 /ε 0 will be used to denote the dielectric permittivity and magnetic permeability of free space, the frequency and the angular frequency, the wavelength, the free space wavenumber and intrinsic impedance, respectively. For high conductivity and thin disks, i.e., for σ i ωε 0 , τ i a and τ i λ, the disks can be approximated with infinitesimally thin flat surfaces located at the abscissas z = z i [6]. Effective surface current densities, J i (ρ, φ), are excited on the disks by an impinging plane wave, ,ρ,φ,ẑ denote the unit vectors in the ρ, φ, z directions, respectively, and the incidence angles θ 0 and φ 0 are defined in Figure 1, which, in turn, generate scattered fields, (E sc i (r), H sc i (r)), such that the overall scattered field (E sc (r), H sc (r)) = (E sc 1 (r) + E sc 2 (r), H sc 1 (r) + H sc 2 (r)) satisfies the conditions of a lossy Fabry-Perót resonator realized with two equiaxial thin resistive disks with the same radius. As will be clearly shown: (1) The proposed method is fast convergent and allows one to accurately and efficiently reconstruct both near-field and far-field parameters; (2) The resonance frequencies of the structure are simply individuated by the peaks of the total scattering cross-section (TSCS) and the absorption cross-section (ACS), which are expressed in closed form; (3) The proposed method drastically outperforms the general-purpose commercial software CST Microwave Studio (CST-MWS) in terms of both CPU time and memory occupation. This paper is organized as follows. Section 2 is devoted to the description of the proposed method. In Section 3, the obtained numerical results are presented, and the conclusions are summarized in Section 4.

Formulation of the Problem and Proposed Solution
Two equiaxial resistive disks with the same radius a, thickness i  and electric con- will be used to denote the dielectric permittivity and magnetic permeability of free space, the frequency and the angular frequency, the wavelength, the free space wavenumber and intrinsic impedance, respectively. For high conductivity and thin disks, i.e., for note the unit vectors in the ,, z  directions, respectively, and the incidence angles 0  and 0  are defined in Figure 1, which, in turn, generate scattered fields,  A uniquely solvable boundary value problem for the Maxwell equations [8] is obtained by imposing the radiation condition, the local power boundedness condition and the following generalized boundary conditions on the disks [6]: for ρ ≤ a, 0 ≤ φ < 2π and i = 1, 2, where (E(r), H(r)) = E inc (r) + E sc (r), H inc (r) + H sc (r) denotes the total field and R i = 1/(σ i τ i ) is the electric resistivity of the i-th disk.
It is well-known that the problem at hand can be equivalently formulated in terms of a system of 2-D integral equations for the effective surface current densities [9]. Moreover, the radiation condition can be automatically satisfied by the proper selection of the Green's function of the problem whereas the local power boundedness condition is guaranteed by the edge behavior prescribed for the fields.
Taking advantage of the revolution symmetry of the problem at hand, the fields can be expanded in Fourier series, i.e., In this way, the orthogonality of the Fourier series terms allows to reduce the system of surface integral equations to an infinite set of independent systems of 1-D singular integral equations for the harmonics of the effective surface current densities. It is possible to show that the system of equations associated to the n-th harmonics of the currents can be written in the VHT domain as in the following: is the VHT of order n (VHT n ) of the n-th harmonic of the current on the i-th disk, H (n) (wρ) is the kernel of the VHT n [44], with k 2 0 − w 2 = −j −k 2 0 + w 2 is the spectral domain Green's function of the problem at hand, δ i,j is the Kronecker delta, I is the identity matrix, and E inc(n) (ρ, z i ) is the n-th harmonic of the incident electric field in the i-th disk plane.
By means of the Helmholtz decomposition, the n-th harmonic of the current on the i-th disk is replaced by the superposition of a surface curl-free contribution, J i,C (ρ), and a surface divergence-free contribution, J (n) i,D (ρ) [45]. This is an advantageous choice because it allows us to obtain scalar unknowns in the spectral domain. Indeed, it is simple to show that the VHT n of J (n) i,T (ρ) with T = C, D are scalar functions [41], i.e., Now, the system of Equation (3), which cannot admit a closed form solution, has to be solved by resorting to a discretization scheme. It is important to note that fast convergence and analytical regularization can be simultaneously achieved by means of the Galerkin method providing suitable expansions of the unknowns.
In order to guarantee fast convergence, orthogonal bases of the functional spaces at which the unknowns belong are the best choice. Such spaces can be completely characterized by the edge behavior [47] and the behavior around the origin of the harmonics of the currents, and observing that the sources are off the disks: i,t (ρ) are well behaved functions. By selecting the following expansion functions: where J n (·) is the Bessel function of the first kind [48], which are orthonormal functions on the interval (0, +∞) with the weight function w 2p T −1 [49], it is possible to demonstrate that the complete and non-redundant Neumann series [50] J where γ , p C = 3/2 and p D = 1, reconstruct the physical behavior in (8) [40].
On the other hand, by extracting the asymptotic behavior of the kernel for i = j in (3), the system of Equation (3) can be rewritten as in the following: for ρ ≤ a and i = 1, 2. The discretization of (12) by means of the Galerkin method with the expansions (10) leads to a matrix operator equation. It is simple to show that the matrix operator associated with the first term at the left-hand side of (12) is diagonal due to the orthonormality property of the expansion functions. The general element of the matrix operator A (n) associated with the second term at the left-hand side of (12) can be expressed as a linear combination of integrals of the kind Looking at Formula (13), it is clear as the integrands of the integrals associated to the mutual contributions, i.e., for i = j, have an asymptotic exponential decay related to the distance between the disks, z i − z j . On the other hand, the elements for i = j coincide with the ones obtained for a single resistive disk analyzed in [40]. Hence, following the same line of reasoning presented in [40], it can be shown that A (n) is a compact operator in l 2 . To conclude, the obtained matrix equation is of the Fredholm second-kind in l 2 because the free term has a bounded l 2 -norm [39,40].

Numerical Results
The aim of this section is the validation of the method proposed in Section 2. The approximate solution obtained by truncating a Fredholm second-kind matrix equation converges to the exact solution of the problem as the truncation order increases. In principle, the accuracy of the solution is limited only by the machine precision; in practice, it depends on the accuracy in the numerical evaluation of the coefficient matrix. It will be shown in the following that the solution proposed in this paper, together with the analytical technique shown in [40] devised to efficiently evaluate the 1-D improper integrals of oscillating and algebraic decaying functions filling the coefficient matrix, makes it possible to accurately reconstruct the solution with very low CPU time and memory occupation. It is worth noting that all the simulations are performed by means of an inhouse software code running on a laptop equipped with an Intel Core i7-10510U 1.8 GHz, 16 GB RAM.
The rate of convergence will be analyzed by introducing the following relative computation error as a function of the truncation order of the matrix equation: where 2N − 1 is the number of cylindrical harmonics considered according to the estimation formula in [51], · denotes the Euclidean norm, and x (n) M is the vector related to the first M expansion coefficients of the general harmonics of the unknown functions.
In Figure 2, err N (M) is plotted for different values of the radius (a = λ/2, λ, 2λ) and the distance between the disks (d/a = 0.1, 1, 10, ∞). As can be clearly seen, the error is substantially independent of the distance between the disks and it depends slightly on the radius of the disks in all the examined cases. It is apparent as the convergence is really very fast. Indeed, in the worst case examined (a = 2λ and d/a = 0.1), a relative error less than 1% is obtained for M = 6 and N = 15 with a computation time (t) of only 9 s. Moreover, more accurate solutions can be quickly obtained. Indeed, a relative error less than 0.1% and 0.01% is obtained for at most M = 11, t = 25 s and M = 22, t = 65 s, respectively. In the examples examined in the following M will be chosen in order to guarantee a relative error less than 1%. In Figure 3, for the sake of completeness, the reconstructed behavior of the components of the effective surface current densities for d/a = 1 and for different values of the radius (a = λ/2, λ, 2λ) is shown. Moreover, in Figure 4, the behavior of the components of the current on the disk 1, plotted for a = λ and for different values of the distance between the disks (d/a = 0.1, 1, 10, ∞), clearly shows that the proposed method tends to reconstruct the solution obtained in the case of a single disk (d/a = ∞) as the distance between the disks increases. shows that the proposed method tends to reconstruct the solution obtained in the case of a single disk ( da) as the distance between the disks increases.
It is interesting to observe that the far scattered electric field can be expressed in closed form as a function of the unknowns in the spectral domain by means of the stationary phase method [ can be immediately evaluated once the solution of the problem is known. In Figure 5, BRCS is plotted for the cases considered in Figure 3, i.e., for 1 da and for 2, ,2 a     . Comparisons with CST-MWS are provided for all the examined cases. As can be clearly seen, the agreement is really very good. However, the proposed method drastically outperforms the Integral Equation solver of CST-MWS, which requires tens of thousands of cells and a few minutes to find a reasonable approximation of the solution.
It is interesting to observe that the far scattered electric field can be expressed in closed form as a function of the unknowns in the spectral domain by means of the stationary phase method [9]: i,D (k 0 sin θ)e jz i k 0 cos θ .
As a result, the bistatic radar cross section (BRCS), can be immediately evaluated once the solution of the problem is known. In Figure 5, BRCS is plotted for the cases considered in Figure 3, i.e., for d/a = 1 and for a = λ/2, λ, 2λ. Comparisons with CST-MWS are provided for all the examined cases. As can be clearly seen, the agreement is really very good. However, the proposed method drastically outperforms the Integral Equation solver of CST-MWS, which requires tens of thousands of cells and a few minutes to find a reasonable approximation of the solution.  A resonant behavior can be observed due to the multiple reflection between the disks, i.e., the system of two stacked resistive disks works like a lossy Fabry-Pérot resonator. One way to individuate the resonance frequencies on the natural modes consists in searching for the peaks of TSCS and ACS, which are defined as in the following: where abs P is the power absorbed by the disks and   denotes the real part of a complex number. As clearly shown in [43], by means of Parseval's formula [52], Weber-Schafheitlin discontinuous integral [49] and forward scattering theorem [53], even TSCS and ACS can be expressed in closed form and then quickly evaluated.
In order to clearly appreciate the resonant behavior of the proposed structure, it is assumed 12 1 RR    , i.e., the losses are reduced with respect to the cases examined A resonant behavior can be observed due to the multiple reflection between the disks, i.e., the system of two stacked resistive disks works like a lossy Fabry-Pérot resonator. One way to individuate the resonance frequencies on the natural modes consists in searching for the peaks of TSCS and ACS, which are defined as in the following: where P abs is the power absorbed by the disks and {·} denotes the real part of a complex number. As clearly shown in [43], by means of Parseval's formula [52], Weber-Schafheitlin discontinuous integral [49] and forward scattering theorem [53], even TSCS and ACS can be expressed in closed form and then quickly evaluated.
In order to clearly appreciate the resonant behavior of the proposed structure, it is assumed R 1 = R 2 = 1 Ω, i.e., the losses are reduced with respect to the cases examined above, providing the conditions for sufficiently high Q-factors of the resonant modes. Moreover, orthogonal incidence, i.e., θ 0 = 0 • , is considered so that only the harmonics for n = ±1 contribute to the representation of the solution. In Figure 6, TSCS and ACS for d/a = 1 are plotted for varying values of the normalized frequency, k 0 a. An almost periodic sequence of peaks can be clearly identified. In Figures 7 and 8, the near E-field in the xz plane and z/d = −0.5 plane is shown at the resonant frequencies in the two periods marked by the rectangular windows in Figure 6. It is interesting to observe that the field behavior is reconstructed by using at most 10 expansion functions for each unknown. Moreover, a given number of oscillations along the z axis is associated with each period, whereas the resonances in a period differ from each other by the number of oscillations along the radial direction. As expected, the resonances associated with the most pronounced peaks show a better confinement of the field between the disks.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 16 above, providing the conditions for sufficiently high Q-factors of the resonant modes.
Moreover, orthogonal incidence, i.e., 0 0   , is considered so that only the harmonics for 1 n  contribute to the representation of the solution. In Figure 6, TSCS and ACS for 1 da are plotted for varying values of the normalized frequency, 0 ka. An almost periodic sequence of peaks can be clearly identified. In Figures 7 and 8, the near E-field in the xz plane and 0.5 zd plane is shown at the resonant frequencies in the two periods marked by the rectangular windows in Figure 6. It is interesting to observe that the field behavior is reconstructed by using at most 10 expansion functions for each unknown. Moreover, a given number of oscillations along the z axis is associated with each period, whereas the resonances in a period differ from each other by the number of oscillations along the radial direction. As expected, the resonances associated with the most pronounced peaks show a better confinement of the field between the disks.  above, providing the conditions for sufficiently high Q-factors of the resonant modes.
Moreover, orthogonal incidence, i.e., 0 0   , is considered so that only the harmonics for 1 n  contribute to the representation of the solution. In Figure 6, TSCS and ACS for 1 da are plotted for varying values of the normalized frequency, 0 ka. An almost periodic sequence of peaks can be clearly identified. In Figures 7 and 8, the near E-field in the xz plane and 0.5 zd plane is shown at the resonant frequencies in the two periods marked by the rectangular windows in Figure 6. It is interesting to observe that the field behavior is reconstructed by using at most 10 expansion functions for each unknown. Moreover, a given number of oscillations along the z axis is associated with each period, whereas the resonances in a period differ from each other by the number of oscillations along the radial direction. As expected, the resonances associated with the most pronounced peaks show a better confinement of the field between the disks.

Conclusions
In this paper, the analysis of the scattering from a two stacked thin resistive disks resonator has been successfully carried out by means of a generalized version of the Helmholtz-Galerkin regularizing technique recently developed by the author for the analysis of PEC, resistive/graphene, and dielectric disks. It has been clearly shown that the proposed method is accurate and efficient in reconstructing the near-field and the far-field parameters even at the resonance frequencies of the structure, and drastically outperforms CST-MWS. Future perspectives are the focus on the physical issues of the considered problem and the generalization of the proposed approach to the analysis of arrays of non-equiaxial resistive/dielectric/composite disks with different radii in both homogeneous and layered media.

Data Availability Statement:
The data presented in this study have been obtained by means of an in-house software code implementing the proposed method.