Rigorous Coupled Wave Analysis of Parallel-Plate Waveguides Loaded With Glide-Symmetric Dielectric Structures

Recently, there has been a growing interest in dielectric antennas and electromagnetic (EM) components based on periodic structures, since such components can be realized using new additive manufacturing technologies. The design of these components starts from the geometry of the unit cell, which determines not only the value of the effective refractive index but also the parameters of the entire considered component, like for example, the bandwidth of operation. For this reason, periodic structures with higher symmetry are of particular interest, since they can exhibit a wider bandwidth and other improved parameters. To analyze and design such components, we present a semianalytical method for calculating the guiding properties of parallel-plate waveguides (PPWs) loaded with periodic glide-symmetric dielectric structures. It is shown that the presence of glide symmetry can be simply included in the dispersion characteristic equation based on the transverse resonance condition. The derived expressions reveal that due to glide symmetry only the Bloch-Floquet modes of the same index parity interact and contribute to the dispersion properties. The accuracy of the discussed method is verified through the analysis of several 1-D and 2-D structures suitable for the design of planar lens antennas with a large frequency band of operation.


I. INTRODUCTION
P ERIODIC structures are commonly used in microwave/ optical regimes for a wide range of applications because they allow the propagation of electromagnetic (EM) waves to be tailored by modifying the geometric parameters of the unit cells.If the unit cell is invariant to the composition of a translation and another geometric operator, it is said to have higher symmetry [1], [2], [3], [4].The introduction of higher symmetries into the unit cell, such as glide symmetry, allows some propagation properties to manifest independently of the specific unit cell geometry [1], thus encouraging researchers to find new functionalities of EM periodic structures.
In the last decade, recently developed technologies such as additive manufacturing have inspired the development of The authors are with the Faculty of Electrical Engineering and Computing, University of Zagreb, 10000 Zagreb, Croatia (e-mail: dubravko.tomic@fer.hr;zvonimir.sipus@fer.hr).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TMTT.2024.3352488.
Digital Object Identifier 10.1109/TMTT.2024.3352488new EM components [5], [6].For example, in pure dielectric structures, the appropriate refractive index profile is obtained by means of division of the structure into periodic unit cells where the appropriate value of the refractive index is achieved by a different proportion of typically two dielectrics (one of them is most often air).At the same time, the geometry of the unit cell affects parameters such as the achievable range of the effective refractive index, the bandwidth of operation, and the isotropic properties of the considered structure (especially important in lens antenna design).Glide-symmetric dielectric structures have been demonstrated to have advantages compared to standard periodic structures and have been used in the realization of, for example, various 2-D lens antennas [7], [8].
When analyzing and designing such structures, it is desirable to have a specialized computer program that can efficiently and accurately calculate the parameters of the considered structure and that will also provide insight into physical phenomena present.
More specifically, a glide-symmetric structure is a periodic structure invariant to the combination of translation operation for half a period along one or two directions and reflection with respect to the so-called "glide" plane.This can be written using the glide operator-mathematical description of a geometrical transformation.For example if the glide plane is plane z = 0 and structure is 2-D periodic with periods d x and d y in xand y-directions, the glide operator is Structure which is invariant to this operator is said to be glide-symmetric.An example of such structure made out of dielectric material is shown in Fig. 1.Parallel-plate waveguides (PPWs) loaded with glide-symmetric periodic structures, like the one in Fig. 1, became popular over the past couple of years for the realization of planar gradient index (GRIN) lens antennas, mostly because the dispersion properties of glide-symmetric structures enable the realization of frequency wideband lenses [9].So far, such planar GRIN lenses for microwave frequencies have been proposed with fully metallic [10], [11], [12], [13], [14], printed circuit board [15], [16], [17], and fully dielectric 3-D printed unit cells embedded in PPW [7], [8], [18].In this article we consider the latter solution, i.e., the 0018-9480 © 2024 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.loading GRIN media can be manufactured by commercially available additive manufacturing technology.Such technology offers multiple advantages, including quick and affordable prototyping, good performance for 5G frequency ranges, improved isotropy, and reduced minimal achievable effective refractive index due to glide symmetry.Apart from lens antennas, filtering structures and leaky-wave antennas can also be designed using glide-symmetric structures and structures with broken glide symmetry (which makes it possible to obtain additional stopbands) [19], [20], [21].These novel possibilities have also encouraged the development of novel analysis methods that will enable computationally efficient design of such structures.
Dispersion analysis of periodic structures is of great importance in wave propagation studies, since it provides insight into the behavior of the artificial material made of such periodic structures.From the dispersion diagram, it is possible to evaluate important parameters of the analyzed structure, such as the effective refractive index, the phase and group velocities, the level of isotropy/anisotropy, the presence of bandgaps, and the radiation and attenuation losses [22], [23].Analytic methods based on the mode matching technique have been developed for holey metallic glide-symmetric structures, optionally the presence of a homogeneous dielectric is included.Originally, the method was developed for the square hole profile [24], but it was generalized in [25] to any hole profile for which the mode distribution can be found.So far, the generalized formulation has been applied to circular [25], triangular [26], and elliptical holes [27].Effective refractive index of glide-symmetric configurations can also be obtained using the quasi-static homogenization techniques proposed in [28], [29], and [30].All of these methods are applicable for metallic periodic structures, but for dielectric structures it is possible to either use commercial eigenmode solvers or the recently proposed, semianalytic method, called the multimode transfer-matrix method (MMTMM) [31], which was later linearized and extended to 3-D structures in [32].The main drawback of this method is its dependence on commercial general EM simulators needed to obtain the scattering parameters of the unit cell.
In this article, we propose a semianalytic method for dispersion analysis of PPWs loaded with glide-symmetric periodic structures made out of dielectric materials.The formulation is based on rigorous coupled wave analysis (RCWA) and it does not require knowledge of any modal EM field distributions dependent on the dielectric spatial modulation.RCWA was originally proposed in [33] for the analysis of scattering properties of optical diffraction gratings.It was further developed in [34] to enable scattering and guiding analysis for arbitrary direction of the incident wave.It is a common tool for analyzing various optical structures such as distributed feedback and distributed Bragg reflector lasers [35], [36], [37], [38].
This article is structured as follows.Section II gives a detailed outline of the RCWA method for 1-D and 2-D periodic geometries, as well as the relationship with the transverse resonance condition, required for obtaining the characteristic equation.In Section III, the glide symmetry is included in the previously obtained transverse resonance condition.Finally, in Section IV, various 1-D and 2-D structures are analyzed using the proposed method.

II. RIGOROUS COUPLED WAVE ANALYSIS METHOD
The structure studied in this article is comprised of two dielectric layers with a periodic spatially modulated permittivity distribution ε(x, y) = ε(x + d x , y + d y ).The two layers have same permittivity distribution but with a spatial shift of half a period between them, thus making the complete structure glide-symmetric.These dielectric layers are inserted between two infinite perfect electric conductor (PEC) plates forming PPW.An illustration of a section of the structure is shown in Fig. 1, where in Fig. 1(a) the boundaries between the dielectric regions are carved out to provide a clear view of the dielectric distributions.
The rigorous coupled wave method, as outlined in [34], can be used to obtain the guiding and scattering characteristics of optical diffraction gratings [34], [39], [40], [41].The method is based on considering EM waves that may appear in every layer of a multilayered structure.These waves are assumed to be traveling in direction normal to the layering of the structure, z-direction in Fig. 1, while in the transversal plane the fields are varying according to the Bloch-Floquet theorem.The Fig. 2. Equivalent transverse networks for the coupled wave analysis of PPWs loaded with periodic, arbitrarily modulated, glide-symmetric, and dielectric structures shown in Fig. 1.The gray boxes correspond to different layers of dielectric structure and the ports are labeled based on the corresponding wavenumber k z,i .solution is then obtained by determining a linear combination of these waves that satisfies the boundary conditions.
Using the network representation for RCWA, as outlined in [39], the structure from Fig. 1 can be represented as a transverse multiport network shown in Fig. 2. The equivalent network is comprised of two infinite-port components, where each port corresponds to one space harmonic propagating in the ±z-direction with the wavenumber k z,i .Due to their periodic nature, the dielectric layers perform energy coupling of each space harmonic k z,i entering at one side into all other possible space harmonics leaving the structure (or reflecting back) at the other side.All one-side ports of the first component and all ports on the other side of the second component have short-circuit termination corresponding to the top and bottom PEC walls of PPW.Nonterminated ports of output/input sides of components are connected to each other on the basis of a common Vertically propagating space harmonic.Note that the wave numbers k z,i of the space harmonics are defined by the periodic structure itself, and procedure for obtaining them as well as coupling between them is given in Sections II-A and II-B.Once the wave numbers are obtained, it is possible to express the tangential electric and magnetic fields in terms of forward-and backward-traveling space harmonics of unknown amplitudes.
To determine the guiding characteristics of the structure under analysis, the transverse resonance condition is applied.This condition requires knowledge of the input impedances in the +z and −z directions at arbitrarily chosen plane z = z 0 .For the considered case, plane z = 0 is selected and thus the input impedance of a thin single dielectric layer backed by a PEC has to be evaluated.The input impedance for the considered case is determined in Section II-C.

A. Field Representations in 1-D Periodic Dielectric Layer
In this section, we will analyze 1-D periodic dielectric structures, i.e., a layer with periodic modulation only in the x-direction, embedded inside the PPW whose height is significantly smaller than the wavelength.In case of y-independent EM field distribution (k y = 0), the lowest order propagating modes in such structures are the TM x modes, i.e., modes with H x = 0.The axis in which the structure has finite thickness (the z-axis as shown in Fig. 1) will be called the normal direction, thus the plane orthogonal to it is the tangential plane.In order to use the transverse resonance method, tangential (x y plane) components of electric and magnetic fields need to be obtained.
Field components in this single dielectric layer with periodically modulated permittivity ε(x) = ε(x +d x ) can be expressed as a series of Bloch-Floquet harmonics where i ∈ {x, y, z} specifies the Cartesian component of the field and k x,m = k x,0 + 2π m/d x , m ∈ Z is the Bloch wavenumber.Since the structure has finite thickness in the zdirection, the EM field variation in the z-direction is described with E i,m (z) and H i,m (z), i.e., with terms that are different for each mth spatial harmonic and include dependence on the normal component k z of the wave vector.
The TM x nature of the modes propagating in the x-direction (i.e., k y = 0) implies that the only nonzero field components are E x , E z , and H y .Connections between these components are obtained from the source-free Maxwell's curl equations.Inserting the Bloch-Floquet series ( 2) and ( 3) into the equation Taking the inner product over the unit cell with the Bloch-Floquet mode exp + jk x, p x , p ∈ Z, leads to a system of equations where the orthogonality property of Bloch-Floquet modes is used Here, δ mp is the Kronecker delta function.
Similarly, inserting series ( 2) and ( 3) alongside with and In the first equation, special care is taken to avoid products of complementarily discontinuous functions which have Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
been shown in [42] to have an impact on the convergence of modal expansions.When taking the inner product with exp + jk x, p x , p ∈ Z, it is important to note that the permittivity distribution ε(x) is also a function of x.This results with two sets of equations where η mp and ϵ mp are defined as Before rewriting these equations in a vector-matrix form, an appropriate notation is introduced.All spatial harmonics of the EM field components are grouped into a vector with the form φ i (z) = φ i,−P , . . ., φ i,P T , φ ∈ {E, H}, i ∈ {x, y, z}.
Additionally, three square matrices are introduced: η = η pm , ϵ = ϵ pm and a diagonal square matrix These vectors and matrices should be infinite since p = −∞, . . ., ∞ and m = −∞, . . ., ∞, but for implementation purposes they are truncated to P Fourier modes, thus the vectors are of dimension [(2P + 1) × 1] and the matrices are of dimension [(2P + 1) × (2P + 1)].Using this notation, the sets of ( 5), (9), and (10) are rewritten as Inserting ( 16) into ( 14), then taking the derivative with respect to z and replacing dH y dz with (15), results in the wave equation This equation describes the couplings of all space harmonics inside the periodic layer.The solution of this vector-matrix equation can be represented using the set of eigenvalues and eigenvectors corresponding to the matrix A as follows.The matrix A can be factorized as A = QDQ −1 where Q is a square matrix whose ith column is the eigenvector q i of A, and D is the diagonal matrix with the corresponding eigenvalues d ii .With this result, the tangential components of the electric field can be written as where exp(− jK z z) is a diagonal matrix whose ith element equals exp − j √ d ii z , i.e., it is related to ith diagonal element d ii of the eigenvalue matrix D. Vectors f and g are amplitudes of forward-and backward-propagating waves along the z-axis represented using the Bloch-Floquet expansion.
Inserting this result into ( 14), together with ( 16), yields the expression for tangential magnetic field component Now tangential components of both electric and magnetic fields inside the periodic dielectric layer are expressed in terms of waves propagating along the z-axis, and input impedance expression can be obtained.

B. Field Representations in 2-D Periodic Dielectric Layer
In this section, the analysis procedure for 1-D structures is extended to dielectric layers with periodicity in two directions (2-D).In the same way as before, a single dielectric layer with finite thickness in the z-direction and infinite in the x y plane is analyzed with the goal of evaluating tangential components of electric and magnetic fields.Unlike for 1-D case, 2-D case requires the presence of both TE and TM modes in the traveling direction, thus all six components of electric and magnetic fields have to be taken into account.The considered dielectric layer has spatially modulated periodic permittivity distribution ε(x, y) = ε(x + d x /2, y + d y /2).Field components can again be expressed as a series of Bloch-Floquet harmonics where i ∈ {x, y, z} specifies the Cartesian component of the field, ρ = x x + y ŷ is the vector in the tangential x y plane, k t,mn = k x,m x + k y,n ŷ is the tangential part of the wave vector with k x,m = k x,0 + 2π m/d x and k y,n = k y,0 + 2πn/d y , m, n ∈ Z, being the Bloch wave numbers in xand y-directions, respectively.Inserting these expressions into the source-free Maxwell's curl equations now yields six equations.Taking an inner product over the unit cell with the Bloch-Floquet mode exp + jk x, p x + jk y,q y , p, q ∈ Z, leads to a system of equations.The first triplet of equations obtained from ∇×E = − jωµ 0 H is and the second triplet of equations obtained from m n η pq,mn jk x,m H z,mn (z) where η pq,my and ϵ pq,mn are given by Rewriting these equations in a vector-matrix form allow us to formulate the problem in a form of the wave equation.Notation from Section II-A has to be varied in order to account for both TE and TM modes.All vertical space harmonics in the z-direction are placed in a supervector Matrices η and ϵ are square matrices of order (2P + 1)(2Q + 1) and are filled according to indexation where ψ ∈ {η, ϵ}.Using these notations, ( 23)-( 28) can be rewritten in the vector-matrix form as where K x and K y are block diagonal matrices with (2P + 1) blocks k i, p along the diagonal Blocks k i, p are diagonal square matrices of order (2Q + 1), and depending on i = x or i = y they are populated differently Manipulating the set of equations in analogous way as in Section II-A, the following form of the wave equation is obtained: The solution of this equation is again obtained by decomposition A = QDQ −1 , where A is the product of the square matrices on the right side of (41) scaled with k 2 0 .Expressions for the tangential field components can be written in the same form as in (18) and (19)  Additionally, it is necessary to modify the expression for P that in the 2-D case is given with the equivalent equation to (20) in which the first matrix is replaced by the first matrix from the right side of (41).

C. Input Impedance
From the expressions for the tangential components of the electric and magnetic fields inside the periodic dielectric layer (18), (19), the input impedance for PEC backed periodic dielectric layer can be obtained.
For the periodic dielectric layer placed between z = 0 and z = t g , as in Fig. 1(a), the tangential components of the electric field vanish at the plane z = t g due to PEC termination.Therefore, tangential E-field components, given by (18), at the plane z = t g satisfy where 0 is the null-vector.The reflection matrix that relates forward-and backward-propagating waves along the z-axis can be expressed as One should note that (43) ensures that the PEC boundary condition is satisfied.
The transverse resonance condition is applied in the symmetry plane z = 0, and the input impedance matrix for the periodic dielectric layer toward the PEC termination needs to be obtained.The input impedance matrix is obtained from Using the reflection matrix and expressions for electric and magnetic fields in E t (z = 0) = Z in H t (z = 0), the expression for input impedance matrix of a single periodic dielectric layer with thickness t g backed by PEC wall is obtained and it is given by The input impedance matrix Z in does not depend on the direction along z-axis, i.e., it is irrelevant if the PEC termination is placed on a plane z > 0 or on a plane z < 0 where superscripts (up) and (dn) denote the direction along z-axis in which the impedance matrix is perceived.It is important to note that if the same structure is analyzed in both directions, the input impedance matrices are identical in .Applying the transverse resonance method at plane z = 0, the characteristic equation is obtained det

III. GLIDE-SYMMETRIC STRUCTURES
In this section, the derivation of the characteristic equation of PPW loaded with a periodic structure that includes glide symmetry is presented.The obtained expression gives an interesting insight into the coupling of modes at the glide interface.
As stated before, the EM problem can be described by the equivalent transverse network shown in Fig. 2 comprised of two infinite-port components representing the dielectric periodic layers.If the top layer was mirror-symmetric instead of glide-symmetric, the image theorem could be applied at z = 0, thus reducing the two component network to a single component network with short-circuit connection on the top and bottom, i.e., at z = t g and z = −t g .In the rest of this section, we construct the single network component model that will account for glide symmetry as well.

A. Effect of Glide Symmetry
Presence of symmetries in any problem can often be used to simplify it, and the same applies here.Before proceeding with the RCWA method, it is analyzed how glide symmetry reduces the complexity of the problem.The input impedance matrix, given by (44), is a function of Q, P, and K z , all of which are obtained from the matrix A which, in the end, is defined by permittivity matrices η and ϵ.Consequently, only these two matrices have an effect on the input impedance matrix and it is sufficient to analyze how glide symmetry affects them.
The spatial variation of the relative permittivity of the glide shifted layer is given by ε (glide) (x, y) = ε(x − d x /2, y − d y /2), where in the 1-D periodic case the dependence on y is omitted.Matrix η (glide) is obtained from (11) or (29).Evaluating the integrals in the considered equations while replacing ε(x, y) with ε(x − d x /2, y − d y /2) yields η (glide) mn, pq = (−1) (m− p) (−1) (n−q) η mn, pq for 1-D and 2-D cases, respectively, showing how glide symmetry affects matrices.In an analogous way, the same results are obtained for ϵ (glide) mp and ϵ (glide) mn, pq .Therefore, connection between nominal layer and glide-shifted layer permittivity matrices can be written as where ⊙ is the Hadamard product operator [43], and G is Toeplitz symmetric matrix whose elements on the main diagonal and on all even diagonals are 1, and on all odd diagonals are −1.For the 2-D case, the matrix G is 2 × 2 block matrix with each block being equal to the 1-D G matrix.
Using this representation, it is possible to find a relation between all matrices that describe the nominal and glideshifted layers.The characteristic equation of a single layer is described with eigenvector and eigenvalue matrices Q and K z , respectively.The eigenvalue matrix K z will remain identical for both layers since it describes the same problem and the waves that are supported along the z-axis have to be identical, i.e., they have the same wavenumber k z,i .The reflection matrix depends only on eigenvalues, thus it also stays the same for both layers = . (51) The relation between eigenvector matrices can be obtained using the following two properties of Hadamard products with matrix G: where X and Y are arbitrary square matrices, and additionally X is invertible.The proof is given in the Appendix.Now, it is possible to prove that A (glide) = G ⊙ A by using these two properties and the definition equations ( 17) and (41).To obtain the matrix Q, the eigenvalue decomposition of matrix A is performed, and from the uniqueness of this decomposition, the following connections between the eigenvector matrices for the nominal and the glide-shifted layers are obtained: The complete derivation of these equations is also given in the Appendix.
Using this property alongside with (51), (53), and the expression for Z in [see (44)], the connection between input impedance matrices is derived (55)

B. Glide-Symmetric Transverse Resonance Condition
The transverse resonance condition ( 47) is rewritten at the interface between nominal and glide-shifted layers as and is used to obtain the guiding characteristics of such PPWs.All values in the final matrix are defined by geometry, except frequency k 0 and Bloch-wavenumber k x,0 , i.e., the mode propagation constant along the x-direction.Solving (56) by fixing the frequency k 0 and finding the complex root k x,0 = β x − jα x of (56) gives the dispersion properties of the structure.In the 2-D case, it is necessary to fix the relation between k y,0 and k x,0 , in addition to fixing k 0 , when solving (56).
Applying the transverse resonance condition at the glide plane allows simplification of the dispersion equation.In more detail, the sum of two impedance matrices at the glide plane is equal to where J is all-ones matrix having the same dimension as G.By introducing G ′ 1-D/2-D = (1/2)(J + G), the general dispersion equation that takes into account the presence of glide symmetry can be written as The matrix G ′ 1-D is Toeplitz symmetric matrix, which on the main diagonal and all even diagonals has ones and on all odd diagonals has zeroes.For the 2-D case, the matrix G ′ 2-D is the 2 × 2 block matrix, with each block equal to G ′ 1-D .The matrix G ′ 1-D/2-D shows how the presence of glide symmetry modifies the nonglide case, i.e., it gives physical insight into glide symmetry.It shows that only Bloch-Floquet modes of the same index parity interact and contribute to the dispersion characteristics, i.e., odd/even modes interact only with odd/even modes.
Additionally, the 2-D transverse resonance condition can be modified to be used in the analysis of 1-D periodic structures illuminated with oblique incident plane waves.This can be done by replacing the matrix K y with a constant value of k y , which is possible since the 1-D media is uniform in the y-direction.The value of k y is determined by the relation k y = k 0 sin φ, where φ is the angle of propagation of the incident wave relative to the x-axis.

IV. NUMERICAL RESULTS
The proposed formulation can be applied to glide-symmetric structures with any spatial modulation of the dielectric material by following simple steps.
1) Evaluate integrals to fill matrices η and ϵ.This can be done either analytically or numerically.2) For given k t,0 and k 0 , fill in matrix A and then perform eigendecomposition on it to obtain Q and K z .It is important to note that depending on the subroutine used, the resulting eigenvector and eigenvalue matrices are not necessarily ordered properly.3) Using ( 44), for given k t,0 and k 0 , evaluate the input impedance of a single dielectric layer backed by PEC. 4) For the frequency of interest k 0 , solve (58) for the complex propagation constant k t,0 .For the 2-D case, either k x,0 or k y,0 or their ratio has to be additionally fixed to a constant value, and then the complex root k t,0 is found.To validate the proposed analysis method, we apply it with the aim of computing the dispersion diagrams of different 1-D and 2-D periodic structures.The study cases are chosen to demonstrate the versatility of the method, since it can be used to analyze structures with continuous variation of relative permittivity, as well as discrete ones.The results are compared with those obtained with the MMTMM, where unit cells with continuously varying permittivity distributions are discretized in a sufficient number of constant-permittivity dielectric layers until the obtained propagation constants relatively converge.

A. 1-D Glide-Symmetric Dielectric Grating With Rectangular Permittivity Modulation
The first study case is PPW loaded with rectangularly modulated dielectric grating, i.e., glide-symmetric alternating dielectric slabs with relative permittivities ε 1 and ε 2 , with period d x and with a being the width of the first slab.The permittivity distribution ε(x) of the lower layer is defined as The elements inside the matrices η and ϵ for this selection of ε(x) are and they are obtained from (11) and (12).Without loss of generality, air is taken as one of the dielectric materials, that is, ε 2 = 1, and for easier reading Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the permittivity ε 1 will be referred to as ε slab .A unit cell of the structure with air as second dielectric material is illustrated in Fig. 3(a); for easier visualization of glide symmetry, it is illustrated for a/2 ≤ x ≤ d x + a/2.
In Fig. 4, dispersion diagrams for different values of the width of dielectric slabs are presented.Two values of the slab width a are tested: a = 5 mm and a = 10 mm with ε slab = 3 and d x = 15 mm.The effective permittivity is closer to ε slab for structure with a = 10 mm, as expected, since in that case the periodic structure has significantly higher amount of dielectric filling.Specifically, the effective permittivity at lower frequencies is equal to 1.35 ε 0 and 2.02 ε 0 for these two considered cases, respectively.These values correspond to the effective permittivity values of the filling of the parallel plate capacitor-the considered unit cell of the periodic structure.Furthermore, the attenuation constant in the bandgap frequency range is also higher for larger slab widths, although we have considered lossless dielectric material.It seems that a larger deviation from the light line (i.e., a larger effective permittivity) also causes a larger attenuation in the stopband.Even though the dependence of attenuation constant on the slab width appears to be linear, this dependence is much more complex, as will be demonstrated.Fig. 5 shows the dispersion diagrams for different selections of dielectric permittivity ε slab of the dielectric slabs.Three permittivity values are compared: ε slab = 2, ε slab = 3, and ε slab = 4. Slight discrepancy is observed in higher order modes for the highest value of ε slab .This is due to the need of larger number of modes for the result to converge, and this problem is well documented in the RCWA literature [42], [44], [45], [46].Like in the previous case, higher selected permittivities also cause higher attenuation constant in the stopband.
The dependence of the maximum value of attenuation constant and stopband bandwidth on the ratio of dielectric slab width and periodicity is shown in Fig. 6.Special care has to be taken into account when designing such structures for EM bandgap operation (i.e., to prevent propagation of EM waves) since for certain ratios a/d x the attenuation constant is significantly lower and the bandgap is narrow.One possible  explanation for this occurrence is that for a/d x ≈ 0.5 the structure can simply satisfy the required phase shift necessary for propagation and the structure behaves as a homogeneous medium.This property is also present in the second stopband of the equivalent nonglide structure (note that glide symmetry is closing the first stopband [4]).
Finally, it should be noted that in all considered cases there is an excellent agreement between calculated results obtained by RCWA method and MMTMM.All RCWA results are obtained by taking into account nine Bloch-Floquet modes (P = 4 ), which is sufficient to achieve accurate results for all modes below and above the first stopband.To achieve a good convergence up to the beginning of the stopband, it is enough to consider the lowest five Bloch-Floquet modes (P = 2 ).

B. 1-D Glide-Symmetric Dielectric Grating With Sinusoidal Permittivity Modulation
The second 1-D periodic structure analyzed is the sinusoidally modulated dielectric grating shown in Fig. 3(b).The dielectric distribution of the lower layer is given by where ε g and δ are parameters which define the permittivity modulation, and d x is the period.The condition on parameter δ ensures that ε(x) is always greater than 1.It is known that the RCWA method numerically performs better for continuously modulated structures.This permittivity modulation is chosen to demonstrate this feature of the method.As a first step in the implementation of the RCWA method, it is necessary to evaluate the elements of the matrices η and ϵ.By solving the integrals ( 11) and ( 12), the expressions for matrix elements are obtained as where ξ = (1 − 4δ 2 ) 1/2 .The derivation of these expressions is given in the Appendix.Note that the terms in (63) decay exponentially, in contrast to the linear decay in (60) and ( 61).The values of ε g and δ are selected such that the highest values of dielectric permittivity are equal to ε max = 2, ε max = 3, and ε max = 4, respectively.Simultaneously, the lowest value of dielectric permittivity in all cases is equal to ε min = 1.The dispersion diagram obtained with the proposed method is shown in Fig. 7.The results are compared with the results obtained by the MMTMM, where the sinusoidally modulated layer is discretized in 15 layers.As expected, for higher values of ε g , the structure has a higher effective dielectric permittivity, although the modulation factor δ is also increased.Furthermore, the maximum value of the attenuation constant in the stopband is increased for higher values of the modulation factor δ and the maximum permittivity ε max .
For sinusoidally modulated gratings, a smaller number of Bloch-Floquet modes is sufficient to achieve convergence when compared to step gratings with the same ratio ε max /ε min .The results shown in Fig. 7 are obtained using five Bloch-Floquet modes (P = 2 ) for the case where δ = 1/6.In case of δ = 0.25 and δ = 0.3, seven Bloch-Floquet modes (P = 3 ) are used to achieve convergence.This reduction in the number of modes necessary for convergence is due to the fact that ε(x) does not change abruptly.

C. 2-D Glide-Symmetric Dielectric Grid
To validate the 2-D case, the method is applied to rectangular glide-symmetric dielectric grid proposed in [7] and [8].The unit cell with all characteristic dimensions marked is shown in Fig. 8.Such a unit cell is typically used for its guiding characteristics below the first stopband, e.g., for realization of planar lenses.For this reason, in the subsequent analysis, only the real part of k x,0 or k y,0 is discussed.Poyanco et al. [7] used the unit cell with parameters d x = d y = 4.2 mm, t g = 1 mm, and ε slab = 2.25 and to achieve different values of effective refractive index, the parameter a x = a y = a is varied from 0.4 to 2.2 mm.Fig. 9 shows the dispersion diagrams over the irreducible Brillouin zone (shown in the inset inside the figure) for the minimum and maximum values of a x = a y obtained by the proposed method (full line) and by the MMTMM (markers).Again, the agreement between calculated results is excellent.Effective refractive index of PPW loaded with rectangular glide-symmetric dielectric grid shown in Fig. 8. Parameters are identical to parameters from Fig. 10.Solid line corresponds to the → X direction and dashed line corresponds to the → Y direction in Brillouin zone, also illustrated on insets.varied from 0.5 to 2.5 mm, while keeping the other parameters equal to a x = 0.5 mm, t g = 1 mm, and ε slab = 2.25.The irreducible Brillouin zone has four important directions to be analyzed, and they are shown in Fig. 10 (inset).The rectangular unit cells are usually selected with an intention to obtain anisotropic properties, as for example, in the case of the holey fully metallic realization of the equivalent PPW structure [24].However, in our case, there are no anisotropic effects present, as shown in Fig. 11 in which the effective refractive index n eff in two principal directions is given.This can be explained by considering the unit cell in the lowfrequency limit, in which we can treat the structure as a parallel plate capacitor with effective dielectric filling, i.e., as a structure without anisotropic properties.
The results shown in Figs. 9 and 10 are obtained with 25 Bloch-Floquet modes (P = 2, Q = 2).This is consistent with the 1-D case where for good agreement at frequencies below the first stopband it is sufficient to use five Bloch-Floquet modes, and therefore for 2-D structures five Bloch-Floquet modes are considered both in the xand y-directions.
Since the materials used in additive manufacturing have losses, the proposed method is validated against the MMTMM results for the lossy unit cell.The assumed losses of dielectric material are tan δ = 0.004 [7].The structure from Fig. 8 with a y = 2.2 mm is used as a reference case.Only the direction → X is analyzed and the obtained real and imaginary parts of k x,0 are shown in Fig. 12. Above the stopband, two modes close to each other appear, for this reason they are shown in different colors.The value of α x /k 0 in the first passband is approximately 0.0025, which is consistent with the rough estimate α x /k 0 ≈ n eff tan δ/2 = 0.0028.Again, the agreement between the results obtained by RCWA and MMTMM is excellent.
Finally, if we compare the needed computation time, the proposed RCWA computation for 25 Bloch-Floquet modes (P = 2, Q = 2) takes on average 20 s per direction in the dispersion diagram.The MMTMM computation time is primarily dependent on the CST frequency solver and takes approximately 350 s for the complete dispersion diagram.This is for the structure with losses; in the case without losses, the CST simulation in the frequency domain takes up to 900 s.All tests were performed on an AMD Ryzen 9 5950X@3.4GHz with 64 GB of RAM.

V. CONCLUSION
This article presents the extension of the RCWA method to periodic structures with glide symmetry.In addition to accurately calculating the dispersion diagram, the considered analysis method provides physical insight in the coupling between shifted periodic layers, showing that the coupling exists only between the modes of the same index parity.The accuracy of the method is verified by comparing the calculated dispersion diagrams with the diagrams obtained by an independent method (MMTMM).Through several examples, such as corrugated structures for band-gap applications and anisotropic 2-D periodic structures, it was shown that the properties of the dispersion diagram are not intuitive and that it is necessary to rigorously analyze such structures in order to avoid fundamental errors in the design process.The RCWA method itself is very general and additional extensions can be Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
easily made in order to expand its applicability, for example, to analyze open glide-symmetric waveguide structures like leaky wave antennas.

C. Derivation of Matrices for Sinusoidally Modulated Dielectric Grating
Sinusoidally modulated dielectric grating has dielectric permittivity distribution given by ε(x) = ε g 1 + 2δ cos 2π x d x (81) where d x is the period and ε g and δ are modulation constants.
To have a physically realizable grating, δ must be chosen so that ε(x) is not less than 1.The RCWA method requires the evaluation of the matrices ϵ and η.Elements of matrix ϵ are simple to evaluate by direct integration of ( 12) with the Euler exponential form of cosine function (83) This results with symmetric tridiagonal matrix with ε g on the main diagonal and δε g on the superdiagonal and subdiagonal.Analytically, obtaining the elements of η is not so straightforward.First, the expression can be simplified by substituting θ = 2π x/d x into η mp = 1 2π ε g 2π 0 e − j ( p−m)θ 1 + 2δ cos(θ ) dθ. (84) Integration of this expression can be performed in the complex domain using the residue theorem [47].Transitioning into complex domain with substitution z = e jθ transforms the integration path into unit circle |z| = 1 and cos(θ ) → (z + 1/z)/2.Now the integral of interest becomes The function under integral has different number of singularities depending on m and p.In case m ≥ p, there is only one singular point, otherwise there are two.From the residue theorem follows: where ξ = (1 − 4δ 2 ) 1/2 , and Res(z 0 ) is the residue of the considered function at the singular point z 0 .Evaluating the residues results with (87)

Fig. 1 .
Fig. 1.Section of a PPW loaded with periodic, arbitrarily modulated, glide-symmetric, and dielectric structure.Blue and orange sections represent dielectric material with different relative permittivities, but these may be continuous and not only discrete.Certain sections are carved out to provide clear view of dielectric distribution in one unit cell for top and bottom grating.A unit cell of the structure is shown in (b), where the dashed red outline corresponds to the dielectric change in the top layer and the solid blue outline corresponds to the dielectric change in the bottom layer.(a) Isometric view.(b) Top view.
where the vectors E x and H y are replaced with supervectors −E T y E . Here the subscript t generalizes the equation for both 1-D and 2-D cases, where in the 1-D case the tangential components are E x and H y and in 2-D they are in supervector form −E T y E

Fig. 4 .
Fig. 4. Dispersion diagram of PPW loaded with glide-symmetric dielectric grating with rectangular permittivity modulation shown in Fig. 3(a) for different widths of dielectric slabs.Parameters: d x = 15 mm, t g = 1 mm, and ε slab = 3. RCWA (solid lines) and MMTMM (markers) results of the complex propagation constant k x,0 are presented.(a) Normalized real part and (b) normalized imaginary part.

Fig. 5 .
Fig. 5. Dispersion diagram of PPW loaded with glide-symmetric dielectric grating with rectangular permittivity modulation shown in Fig. 3(a) for different dielectric permittivities ε slab .Parameters: a = 10 mm, d x = 15 mm, and t g = 1 mm.RCWA (solid lines) and MMTMM (markers) results of the complex propagation constant k x,0 are presented.(a) Normalized real part and (b) normalized imaginary part.

Fig. 6 .
Fig. 6.Parametric sweep of (a) maximum value of attenuation constant α x /k 0 and (b) stopband bandwidth over the dielectric slab thickness a x for rectangularly modulated glide-symmetric dielectric grating.The shaded area in (b) corresponds to stopband frequency range for a particular a x .The period of the unit cell is d x = 15 mm.

Fig. 7 .
Fig. 7. Dispersion diagram of PPW loaded with glide-symmetric sinusoidally modulated dielectric grating shown in Fig. 3(b).Parameters: d x = 15 mm, and t g = 1 mm.RCWA (solid lines) and MMTMM (markers) results of the complex propagation constant k x,0 are presented.(a) Normalized real part and (b) normalized imaginary part.

Fig. 12 .
Fig. 12. Dispersion diagram of the glide-symmetric unit cell of dielectric square grid enclosed in PPW shown in Fig. 8. Parameters: a x = a y = 2.2 mm, d x = d y = 4.2 mm, t g = 1 mm, ε slab = 2.25 and tan δ = 0.004.RCWA (solid lines) and MMTMM (markers) results of the complex propagation constant k x,0 are presented.(a) Normalized real part and (b) normalized imaginary part.