Spectral Dyadic Green’s Function for Multilayer Periodic Bi-Anisotropic Media

A computational technique for evaluating the spectral dyadic Green’s function utilized in the analysis of multilayer structures composed of periodic bi-anisotropic layers is presented. The Green’s function is achieved with the help of a fully vectorial rigorous matrix formulation based on a generalized transmission line (TL) formulation of Maxwell’s equations. This matrix formulation simplifies the mathematics involved in the modeling of the coupled diffraction orders within a periodic bi-anisotropic layer. To examine the soundness of this modeling, the extracted Green’s function is used in an integral equation for the analysis of metallic gratings on homogeneous/inhomogeneous anisotropic and chiral materials. The resulting integral equation is then solved by the Method of Moments (MoM) formulated in terms of sub-domain basis and Galerkin’s test functions. As verification, several examples are analyzed and the results obtained by this method are compared with those available in the literature or obtained by EM-solvers. The efficiency assessment of this method is carried out by considering its convergence rate and cost-time in terms of truncation orders. It is revealed not only the metallic grating with inhomogeneous periodic bi-anisotropic layers can be analyzed by this method but also analysis time, memory, and CPU occupancies are decreased in comparison with commercial EM-solvers.


I. INTRODUCTION
Computation of the dyadic Green's function is an essential step in solving a scattering problem in electromagnetics [1], [2]. This function is required in solving the scattering problems of surfaces in layered media involving homogeneous isotropic [3], planar stratified structures with homogeneous complex media such as bi-anisotropic [4], gyroelectric [5], anisotropic-uniaxial [6], and chiral [7] materials. In [8], interaction of electromagnetic waves with bi-anisotropic homogeneous slabs is examined and the computation of the dyadic Green's function is carried out with the formulation explained in [4]. The technique explained in [4] is based on the seminal work of [9], and on the mathematical and computational techniques discussed in [10]. Analysis of multilayer structures and metasurfaces by using dyadic Green's function The associate editor coordinating the review of this manuscript and approving it for publication was Ladislau Matekovits . is paramount importance and is encountered in numerous applications [11], [12]. The dyadic Green's function technique has also been applied to the analysis of periodic planar structures such as Frequency Selective Surfaces (FSSs) on anisotropic [13], [14], [15], [16], chiral [17], [18], and ferrite [19] substrates. The aforementioned works and other reported Green's functions in [20], [21], [22], [23], [24], [25], [26] can be used for analyzing the stratified-multilayered structures with homogeneous layers and are not efficient for inhomogeneous periodic anisotropic substrates. Recently, the advantage of using periodic inhomogeneity in the FSS substrate has been studied [27], [28], [29], [30]. It has been demonstrated that mounting the metallic grating on a periodic substrate or so-called artificial anisotropic dielectric [30] can improve the characteristics of the resultant FSS. Drilling periodic holes in the host substrate and/or filling them with other materials changes the effective permittivity and permeability of the substrate. Thus, this measure allows controlling the dielectric VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ and magnetic properties of the FSS. Furthermore, a periodic substrate allows coupling of energy among different diffraction orders, which is not occurred once a homogeneous substrate is used.
Herein, by periodic we mean the boundary between two alternating materials of the inhomogeneous substrate is perpendicular to the surface of metallic grating (like a photonic crystal slab). To exploit the full potential of FSS, accurate and efficient techniques are needed. The conventional Method of Moments (MoM) has been developed for scattering analysis of metallic grating on homogeneous anisotropic substrates [13], [14], [15], [16], [19] and is not applicable for ones on multilayer substrates containing periodic bi-anisotropic inhomogeneities. Therefore, developing a fast solution and full-vector scheme is required to analyze the mentioned complex structures.
To the best of our knowledge, the dyadic Green's function of such an inhomogeneous periodic bi-anisotropic substrate has not been reported in the literature. Unlike homogeneous anisotropic substrates, the diffraction orders are mutually coupled within periodic bi-anisotropic layers and thereby the dyadic Green's function is no longer diagonal. Thus, when substrates with periodic inhomogeneities are involved, evaluation procedure of the dyadic Green's function is different from [3], [14] and complicated such that the available commercial EM-solvers do not offer any analysis for them. Obviously, solving the scattering problem of metallic grating with inhomogeneous periodic bi-anisotropic substrates, an efficient technique has to be developed for evaluating the required dyadic Green's function in the first step.
The basic idea of combining periodic MoM with coupled wave analysis of an artificial anisotropic substrate has been introduced by authors for the analysis of nano-antenna grating with periodic, isotropic, and non-magnetic thin films [30]. The present work concerns with the generalizing this idea for the scattering analysis of FSS which contains metallic grating printed on an electromagnetic crystal [31] with both electric and magnetic anisotropy as well as chiral medium. To this end, a rigorous semi-analytical technique based on an Equivalent Transmission Line Model (ETLM) in the Fourier domain [30], [32] is exploited to obtain the spectral dyadic Green's function. With the help of this technique, the mutual coupling among diffraction orders inside each periodic bi-anisotropic layer is taken into account. The obtained dyadic Green's function is thereafter utilized in an integral equation formed for solving scattering problems of conducting gratings on homogeneous anisotropic, and inhomogeneous periodic bi-anisotropic substrates. The resultant integral equation is consequently solved by the MoM with sub-domain basis and Galerkin's test functions. The novelty of this study is in the application of the equivalent TL modeling [32], [33] for an electromagnetic crystal to evaluate the dyadic Green's function. Note that our method is different from the Effective Medium (EMT) for analysis [34]; the presented formulation serves as intuitive physical explanation. Finally, the efficiency, and accuracy of the proposed FIGURE 1. Typical scattering problem involving a multilayer structure with periodic bi-anisotropic layers. The diffraction fields are calculated as a function of an incident plane-wave. The unit cell of the periodic substrate is shown as an inset. R 1 and R 2 are the regions containing dielectric rods and host medium, respectively. method is examined and it is shown this technique is a fast and efficient solution for scattering analyzing of FSS with homogeneous or periodic bi-anisotropic layers.

II. FORMULATION AND METHOD OF ANALYSIS
A typical geometry of a periodic bi-anisotropic multilayer structure is shown in Fig. 1. As depicted in this figure, a periodic substrate is stacked on a homogeneous substrate. The thickness of the periodic substrate is h, and is composed of dielectric rods (R 1 ) embedded in a background medium (R 2 ). Such a periodic substrate is like that a photonic crystal slab [32].
In order to extract the spectral dyadic Green's function of this structure, we first study the EM field behavior in the inhomogeneous periodic layer, which can be composed of alternating two different bi-anisotropic materials. This layer is characterized by a generalized anisotropic constitutive relation with a Tellegan representation [35]: The 3 × 3 permittivity tensor ε, the permeability tensor µ, and the magneto-electric tensors ξ and ζ are complex valued. Therefore, the following analysis is applicable to artificial anisotropic or homogeneous complex materials such as arbitrarily biased ferrite, chiral, orthorhombic, or tetragonal crystals (bi-anisotropic materials). Fig. 2 illustrates the subdivided multilayer structure with homogeneous or periodic substrates into parallel layers. Although TL model is derived by using an inhomogeneous periodic bi-anisotropic layer (Layer i in Fig. 2) in this section, but, the resultant model should be applied to all of layers (both of homogeneous and inhomogeneous substrates) for diffraction analysis of the multilayer structure. As shown in Fig. 2, inhomogeneous i-th layer is made of two periodically alternating bi-anisotropic materials. Obviously, the periodicity of the structure in the x-and y-directions must be taken into account. For an incident plane wave represented by

A. DERIVATION OF AN EQUIVALENT TL MODEL
where r = xx + yŷ + zẑ is the position vector and L = L xx + L yŷ in which L x and L y are the dimensions of a unit cell in the x-and y-directions, respectively (see Fig. 1). The vector k = k xx + k yŷ denotes the x-and y-components of the wave vector of the incident plane wave. In other words, for a plane wave arriving at (θ, ϕ), this vector is given by: where k 0 is the free-space wave number. The electric E(r) and magnetic H(r) fields in (4) and (5) can be represented by the following pseudo-Fourier series: where (α n , β m ) are defined as By assuming we define a K ×1 matrix [E 1 (z)] whose elements are E (n,m) 1 (z). Here, K = (2N + 1) × (2M + 1), and the coefficient E (n,m) 1 (z) occupies a specific location in the column matrix [E 1 (z)] according to [32]. Similarly, one can define the column matrices ]. To express the constitutive relations, we observe that the structure can be subdivided into z-invariant layers as shown in Fig. 2. In each of these layers, the relative permittivity, permeability, and chirality can be expressed as a bi-periodic function in x and y, i.e., where u, v ∈ {1, 2, 3} and τ ∈ {ε, µ, ξ, ζ } which are constitutive dyadics (i.e., Cartesian tensors) that can be interpreted as 3 × 3 matrices as follows In a homogeneous isotropic or bi-anisotropic layer, the coefficients τ (n,m) uv (z) are reduced to τ uv δ nm in which δ nm denotes the Kronecker delta. For an inhomogeneous periodic one, the coefficients τ (n,m) uv (z) are given by which can be evaluated either analytically or numerically on the unit cell (UC). Once the expansion coefficients τ (n,m) uv (z) are determined, the constitutive relations with Tellegan representation ( (1) and (2)) in the spatial domain can be reformulated as VOLUME 10, 2022 In these expressions u ∈ {1, 2, 3}, whereas * signifies twodimensional convolution, i.e.
The convolution provides the desired expansion coefficients of D(r) and B(r), and can be converted into a matrix expression [32]. To this end, three single-column matrices [D u (z)] and [B u (z)] with u ∈ {1, 2, 3} are introduced. Using these matrices, one can obtain where the matrix [N τ uv ](u, v ∈ {1, 2, 3} and τ ∈ {ε, µ, ξ, ζ }) in the above relations is a K × K Toeplitz matrix with the following entries Now, we can derive the exact TL formulation using the introduced matrix representation. Here, the mentioned solutions given by equations (4) and (5), along with the relations (19) and (20) are substituted in Maxwell's equations. Therefore, three following relations equivalent to the equation ∇ × E = −jωB are obtained: A similar procedure for the equation ∇ × H = jωD will lead to the following three relations: and [C] are 2K × 2K dimensional matrices and are given by the following relations: in which , , where with and with In relations (30)- (33), the remaining matrices are given by:  Fig. 2. Moreover, all the necessary boundary conditions will be satisfied if the continuity of the line voltages and currents is fulfilled. This is because the line voltages and currents in (29) are the tangential field components for every plane parallel to the XY-plane. To apply the derived TL formulation to the modal analysis of a periodic bi-anisotropic layer, one has to decouple the voltages and currents on various lines of the equivalent TL model. After decoupling of the algebraic system represented by (29), the voltages [V (z)] and currents [I (z)] are to satisfy the following relation: where the single-column matrices [p], [q], and the propagation constant in the z-direction (k z ) are to be determined. By substituting (34) in (29), one arrives at: The relation (35) is an eigenvalue problem. By solving this equation, one can find 4K eigenvalues (k zn ) and the corresponding eigenvectors ([p i ], [q i ]) T . Various eigenvalues can be put in two groups: one set of eigenvalues represents waves VOLUME 10, 2022 traveling in the +z direction and the other one includes eigenvalues for the traveling waves in the −z direction. Therefore, the voltage and current in each layer can be expressed in terms of eigen-solutions as where a n and b n are constant complex values. The parameters z i−1 and z i denote the lower and the upper boundaries of the layer i, respectively. Note that the primed parameters (k zn ) and ([p n ], −[q n ]) T denote the traveling waves in the −z-direction. (36) can be rewritten in a normalized compact matrix form By using (38) at the lower interface (z = z 1 ), [ (2) ] can be obtained as follows: (39) where [P 1 ] and [Q 1 ] correspond to the medium existing in the first layer (Layer 1 in Fig. 2). At the uppermost or at an internal interface, e.g. interface z i−1 , one obtains [ (i) ] in layer i as follows: The other matrices in relations (39) and (40) are given by , where the exponential matrices in the above relations are 2K × 2K diagonal matrices. It is obvious that by using this equivalent TL model and the proper boundary conditions, one can compute the electromagnetic field components in each layer for a given incident wave.

B. CALCULATION OF THE SCATTERED WAVE FROM THE METALLIC GRATING WITH THE MULTILAYER PERIODIC BI-ANISOTROPIC SUBSTRATE
In this section, we evaluate the spectral dyadic Green's function of the multilayer periodic bi-anisotropic structure. To do so, it is assumed that a metallic grating is mounted on such an inhomogeneous substrate at the uppermost interface in Fig. 1 (z = z i ). In the Fourier domain, through the matrix equation Unlike for multilayer metallic grating with homogeneous layers, the diffraction orders are coupled within inhomogeneous periodic ones. In other words, the p-th Fourier order of the electric current on the metallic grating can affect the q-th Fourier order of the diffracted electric field. Therefore, the four sub-matrices in [ G] are no longer diagonal. In this case, two admittance matrices looking upward and downward at the metallization screen, which are shown in Fig. 3, are needed to evaluate the spectral dyadic Green's function as follows [30]: In (44), for calculating [ (i+1 * ) ] by means of the developed formulation, it is assumed that the plane wave is incident from the lowermost layer. Therefore, the dyadic Green's function matrix [ G] is evaluated with the help of (42). According  to (41), the calculation of the scattered field also requires to compute the induced surface electric current densities. To this end, the following boundary condition must be satisfied on the conducting patch surface: where E e t = [E e x , E e y ] T and E s t = [E s x , E s y ] T are tangential components of the excited and the scattered electric fields, respectively. Z s is the surface impedance of the conducting patches. Substituting

III. RESULTS
To verify the rigorous full-wave method developed in the previous section, the reflection characteristics of various periodic structures are evaluated.
The computed numerical results are compared with those obtained by a full-wave EM-solver or published in the literature. Note that by using this method, one can calculate the amplitude and phase of the diffraction orders. However, if the higher-order diffraction orders do not carry active power, the transmitted power is just due to the zero diffraction order (L x , L y < λ inc /(1 + sin θ inc ) [14]).
Fourth layer: cylindrical cavities with (ε 11 , ε 22 , ε 33 ) = (3.4, 3.4, 5.12) and (µ 11 , µ 22 , µ 33 ) = (1, 1, 2) which are surrounded by air. Fifth layer: air gap. Sixth layer: ε r = 2.58, µ r = 1. Last layer: Free space). The parameters of a unit cell as demonstrated in Fig. 4(a). The metallic grating is in the interface of fifth and sixth layers. Note that for cavities depicted in Fig. 1, one may obtain a closed form of the coefficients τ (n,m) uv (z) in Eq. (13). For instance, the cavities in the fourth layer of Fig. 4(a) are of cylindrical shapes and their locations are in the center of the unit cell. Thus, for |x| ≤ L x /2 and |y| ≤ L y /2: in which τ uv ∈ {ε uv , µ uv }, (u, v = 1, 2, 3) are elements of permittivity and permeability tensors of cavities surrounded by air in Layer 4 of depicted structure in Fig. 4(a). By substituting (47) in (13), one can obtain τ (n,m) wherein Jinc(x) = x −1 J 1 (x) which J 1 is the Bessel function of the first kind and first order. As can be seen in Fig. 4(b), the calculated reflected power is in good agreement with the one obtained by a full-wave EM-solver based on Finite-Element Method (FEM). But, the EM-solver is only used for a medium with diagonal constitutive tensors. Since this studied example is a periodic structure, the periodic conditions are used for its simulation by FEM-based EM-solver. Therefore, it is not needed to determine the number of elements unlike finite FSS [36]. Table 1 compares the computational efficiency of the proposed method with FEM in terms of computation speed, CPU usage and memory occupancy. In this table, the speed of computation is sum of the needed times for evaluating the dyadic Green's function and MoM solution in each frequency sample. In other words, the calculation approach of them is mixed together. It should be added that our computer code was running on an Intel (R) 7 core CPU computer with a processing capacity of 3.06GHz and 24.0GHz RAM to generate the results reported in this paper. In the second example, an array of PEC Jerusalem-cross on an electrically and magnetically anisotropic substrate is investigated. The definition of angles included in this example is drawn in Fig. 5(a). As seen in the same figure, θ, and φ determine the propagation direction, (ξ ε , η ε , ζ ) are the principal axes for the permittivity tensor (ε) and (ξ µ , η µ , ζ ) are the principal axes for the permeability tensor (µ). ϑ and ϑ are the misalignment angles. The geometrical parameters of the structure shown in Fig. 5(b) are L x = L y = 24mm, l 1 = 15mm, l 2 = 3mm, l 3 = 9mm, and h = 3mm. In this example, the substrate is made out of an electrically and magnetically homogeneous substrate with (ε ξ ξ , ε ηη , ε ζ ζ ) = (3.4, 5.12, 5.12) and (µ ξ ξ , µ ηη , µ ζ ζ ) = (2, 1, 1). The principal axes lie in the xy-plane. Rotation angles ϑ and ϑ for principal axes are assumed 45 • , and 90 o , respectively. In this case, such an air backed substrate is characterized by:   [15]. Now, the algorithm is tested for the case of a TE-polarized plane wave normally incident to an ungrounded periodic structure consists of a 2D lattice of dipole-shaped patches printed on the homogenous ferrite substrate with a permittivity ε r = 12.8, and a saturation magnetization of 4πM s = 1780G. The unit cell configuration and its design parameters are depicted in the inset of Fig. 6. The design dimensions are L x = 7.6mm, L y = 15.2mm, l = 13.5mm, w = 2.54mm, and the substrate thickness of h = 0.2mm. With y-directed biased magnetic field, the permeability dyad has the form [37]: where γ = 2.8MHz/Oe (is the gyromagnetic ratio divided by 2π), f 0 = γ H 0 (is Larmor frequency), H 0 = 5000Oe is the magnetic bias field and f is the frequency of the incident plane wave. In Fig. 5, the computed reflected power is compared with one reported in [19]. In this extraordinary mode, the obtained Larmor frequency is 14GHz, in which µ, and κ have infinite values based on (50). In this singular point and after it, some methods may be unstable and have a highly oscillatory behavior for reflection coefficient responses [38]. As shown in this figure, our proposed method is stable at about this singular point and has a stable behavior after it. For demonstrating stability of the proposed method, this example is deliberately selected. In order to validate the developed method for a periodic structure on a bi-isotropic substrate, the Chiro-FSS investigated in [18], is analyzed as the next example. This structure consists of a two-dimensional array of the cross-shaped patch shown in the inset of Fig. 7 with the parameters L = 10mm, W = 1mm, and L x = L y = 10mm backed by a chiral slab (ε r = 1.06, and |ξ r | = 0.0026 < µ −1 ε [39]) with the thickness of h = 4mm. As shown in this figure, our computed reflected power of Co-polarized wave and the one reported in [18] are compared when the structure is illuminated by a normally incident plane wave. It should be emphasized the mentioned methods in the above can be used for analyzing FSSs on homogeneous anisotropic substrates and are not efficient for metallic grating analysis with inhomogeneous periodic ones as depicted in Fig.1 and Fig. 2.

B. ANALYSIS OF A METALLIC GRATING EMBEDDED IN AN INHOMOGENEOUS PERIODIC ANISOTROPIC UNIAXIAL CHIRALS
For assessment the evaluated dyadic Green's function in the analysis of the multilayer periodic bi-anisotropic structure, a depicted unit cell in Fig. 8(a) is selected as the last example. In this case, L x = L y , and the metallic array consists of square-shaped patches (L = 0.5L x ) sandwiched by two periodic layers in which air-cavities drilled within an anisotropic uniaxial chiral (r = 0.25L x in Fig. 8(c)). It is assumed the optic axis of the uniaxial anisotropy inside the chiral region forms an angle of α with the z-direction in the XZ-plane as shown in Fig. 9(a). Hence, the constitutive tensors for the uniaxial chiral are expressed by: where a ∈ {ε, ξ }, a E denotes the values of the host medium constants in the direction of the optic axis (Extra-ordinary) and a O implies the values in the direction perpendicular to this axis (Ordinary). In this case, it is supposed ε O = 1.2, ε E = 1.35, ξ O = 0.0001 and ξ E = 0.0007 and all layers are nonmagnetic, and lossless. Such a structure is analyzed by means of MoM/ETLM for normal incidence of TM-polarized plane-wave. Fig. 9(b) shows the reflected powers of the (0, 0)-th spatial harmonics of TE-and TM-polarized waves in terms of periodicity variation with α = 20 • . As observed in the same figure, not only the co-polarization (TM waves), but also the cross-polarized (TE waves) appear in the reflected waves due to chirality of the periodic region. The result of Achiral structure analysis is also reported in this figure.
As shown, no cross-polarization appears for null chirality case because no coupling between TE and TM (0, 0)-th modes occurs.

C. EFFECIENCY OF THE PRESENTED METHOD
The developed spectral dyadic Green's function used for the analysis of metallic grating on multilayer, inhomogeneous periodic bi-anisotropic media contains the expansion of the electromagnetic quantities based on coupled plane waves (fields within the medium) and basis functions (induced currents on the metallic grating). This technique is referred to as a semi-analytical method which, in general, is computationally more efficient than conventional methods used for homogenous structures. This fact has been confirmed throughout the numerical investigations performed in this paper. However, in this section, it focuses on the convergence rate of the expansion of electromagnetic fields within the homogeneous/periodic anisotropic substrates and rooftop basis functions. In this case, the mentioned second and last examples with homogeneous and inhomogeneous anisotropic substrates are respectively considered. The computed convergence rate of these examples are visualized in Fig. 10, which shows the truncation errors in terms of the truncation order M = N. Such an error is calculated for the first resonance of Fig. 5(b) and the first peak of blue-line in Fig. 9(b), as shown in Fig. 10(a), and Fig. 10(b), respectively. This relative error is defined as follows: Truncation Error = RP 00 (N) − RP 00 (N − 1) RP 00 (N) (52) in which, RP 00 (N) is the amplitude of the reflected power of (0, 0)-th mode. For increasing the accuracy of the expansion, the radiated field from the smallest element in the patch of unit cell should be modeled. The narrowest dimensions in the patch configuration of the second and last examples are the one-sixteenth and one-twelfth of the lattice constants, respectively. This leads to the truncation order of M = N = 16 for the second example and M = N = 12 for the last example to achieve relative errors less than 0.5%. Fig. 10 also shows the computation times as a function of the truncation order. By examining various examples, it is yielded that the following relations will give an acceptable accuracy: in which, nd Px and nd Py are the narrowest dimensions of the patch in the x-and y-directions, respectively. As shown in Fig. 10, the computation cost may be very high once the high truncation orders should be assumed. Based on (53), when the FSS contains patches with very fine dimensions, using high truncation orders is inevitable. In these cases, using MoM with entire-domain basis functions and Galekin testing functions is expected to be more efficient. Using methods based on Finite differences and Finite Elements are not efficient when the multilayer bi-anisotropic substrates are used due to their high computation costs.

IV. CONCLUSION
A semi-analytical method for a full-wave analysis of metallic grating loaded by multilayer, inhomogeneous periodic bi-anisotropic media has been introduced. The technique utilizes the MoM with sub-domain basis and Galerkin's test functions along with a generalized equivalent TL modeling described in [30] for the evaluation of the dyadic Green's function in the spectral domain. By using this method, various periodic structures on inhomogeneous periodic isotropic, anisotropic, chiral, and bi-anisotropic substrates have been analyzed to obtain their reflection characteristics. In order to verify the results produced by this method, they have been compared with previously published ones and those obtained by a full-wave EM-solver. The presented analysis incorporates arbitrary incident angles and optic axis for the periodic bi-anisotropic substrate. It was shown that the developed technique is valid for the analysis of FSSs on various inhomogeneous anisotropic substrates. It is shown that the reflection characteristics of a unit cell can be controlled by the direction of the optic axis among other parameters. In other words, FSSs design can be accomplished with a larger degree of freedom if the substrate of them is periodic bi-anisotropic.