Method of Moments for the Dispersion Modeling of Glide-Symmetric Periodic Structures

A modeling methodology to obtain the dispersion characteristics of mirror- and glide-symmetric structures is presented. A novel Green’s function is proposed as the integration kernel of the electric-field integral equation solved by the method of moments (MoM). Key aspects of implementation, such as adapting the Ewald acceleration, accurate computation of singular integrals, and a zero-search algorithm to obtain solutions, are presented. The proposed methodology is applied to fully metallic 2-D periodic unit cells with arbitrary geometries. The results of the method are found to be in very good agreement with reference results from the literature. Compared to the conventional MoM analysis, the proposed approach obtains results in half the time and gives additional information about the modal properties.

results in increased attenuation of electromagnetic waves in materials and free space.One way to mitigate increased losses in materials is to utilize fully metallic structures, thus completely removing wave attenuation in dielectric materials.On the other hand, free-space losses can be mitigated with high-gain antennas.
A particular topic of interest, in recent years, to design highly efficient structures, has been the use of glide-symmetric metasurfaces [4].Glide-symmetric structures are periodic and invariant with respect to a translation by half a period and a mirroring.The effects of introducing such higher symmetry on propagation of waves in periodic structures have first been studied in the 1960s and 1970s [5], [6].Recently, their properties have been exploited to design highly efficient fully metallic graded index lens antennas [7], [8], leaky wave antennas [9], [10], [11], and filters [12].Furthermore, glide symmetries have been used for suppression of leakage between two connecting parts of a waveguide [13] and flanges [14].
As a result of increased interest, methods to accurately compute the dispersion diagram (i.e., the propagation constants of Bloch modes) of glide-symmetric structures have also garnered attention from the antenna and microwave communities.A commonly applied technique is mode matching [15], [16], which is particularly suited to holey structures.At discontinuities, the coupling of modes in different parts of the structure is done by enforcing boundary conditions.If the structure has known analytical field expansions, mode matching is a computationally efficient and a straightforward method to implement.Furthermore, at lower frequencies, quasi-static approximation can be used with mode matching to faster obtain the effective refractive index [17].Recently, a hybrid method named multimodal transfer matrix method (MMTMM) has been proposed for the analysis of glide-symmetric structures [18].In this method, waveguide ports are placed at the edge of the unit cell and a full-wave solver is used to compute the scattering parameters.In this step, the definition of multiple modes on each port is necessary for each waveguide in order to obtain correct results.A postprocessing procedure is used to obtain the propagation characteristics of Bloch modes in the structure.The method can find both real and complex solutions and has recently been applied to study the stopband attenuation of glide-symmetric unit cells with circular holes [19].
In this work, we propose a method of moments (MoM) modeling approach tailored to fully metallic glide-symmetric structures.MoM has been previously used for modal analysis of printed structures [20], [21], [22], [23], electromagnetic bandgap structures [24], and 1-D periodic glide-symmetric parallel-plate waveguides [25].In this article, we present a novel integration kernel for 2-D glide-symmetric structures, which reduces the computational domain to the bottom half of the unit cell only.The new kernel is applied to analyze a rectangular holey structure [15] and a circular holey structure [19].The obtained results are compared to a commercial eigenvalue solver and to the MMTMM analysis [19], finding, in both cases, a very good agreement.Several benefits are offered by the method presented here over other approaches.Most importantly, the MoM is able to compute the attenuation constant in the stopbands and for complex modes, both of which are currently unavailable from commercial eigenvalue solvers.Furthermore, the method provides additional information on mode parity and can be easily applied to arbitrary geometries, an advantage over mode matching.As opposed to MMTMM, it does not require truncation of waveguide modes at the edges of the unit cell [19].Hence, less physical insight, necessary to understand which modes to include in the computation, is required to obtain accurate results.Furthermore, since commercial software includes waveguide-port modes by following a strict order of cutoff frequencies and cannot arbitrarily select the modes to retain, using MMTMM can result in an increased computational effort since many higher order modes must be included in the simulation in order to keep a few modes of very high order, as in [26].
This article is organized as follows.In Section II, the integral equation formulation and the novel Green's function are presented.Then, in Section III, the equation discretization, the computation of singular integrals, and the zero-search algorithm are described.In Section IV, numerical results are compared to a commercial eigenvalue solver and MMTMM analysis.The conclusion and perspectives are given in Section V. Preliminary results of the proposed modeling of 2-D periodic glide-symmetric structures were presented in [27].

A. Integral Equation and Periodic Green's Function
Let us consider the task of finding modes in a fully metallic periodic structure whose surface is described with the geometry S. If losses due to finite conductivity of metal are negligible, the surface can be described as a perfect electric conductor (PEC) in free space without significantly impacting the accuracy of the solution.On S, the tangential component of the electric field is zero where ω is the angular frequency, ω = 2π f , and the field is decomposed into the vector and the scalar potentials, A and , respectively.The vector potential represents the contribution to E tan due to the electrical surface current density J and the scalar potential represents the contribution of the surface charge density Both integrals are computed over the entire surface S. In ( 2) and ( 3), ε and µ are the permittivity and the permeability of the surrounding medium, respectively, J(r ′ ) is the unknown electrical surface current density at the source point r ′ , and G(r, r ′ ) is the scalar Green's function at the observation point r.The dyadic Green's function G is defined as a product of the identity dyad I = x x + ŷ ŷ + ẑẑ and the scalar Green's function as Since the structure is periodic and thus spans to infinity in both dimensions, the computation domain can be restricted to the finite size of one single unit cell (the surface S in (2) and (3) becomes the metallic surface of a single cell) if the periodicity is captured in the Green's function.This is achieved by introducing the free-space periodic Green's function (FSPGF) that can be written as a double infinite sum of free-space Green's functions [28] G r, r where Here, k t00 is the wave vector controlling the (possibly complex) propagation between different points in the lattice.The distance R mn is given by where Here, ρ mn describes the translation from the unit cell indexed by m = n = 0 to cell indexed by m and n of the lattice, s 1 and s 2 being the lattice periodicity vectors.

B. Glide-Symmetric Periodic Green's Function
We will now extend the concept of FSPGF to introduce a novel Green's function, which exploits the structure's internal symmetry.For mirror-and glide-symmetric structures, we introduce the scalar higher-symmetric periodic Green's function (HSPGF) where G B represents the bottom array of sources as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and G T represents the top array of sources as The Green's function is represented in Fig. 1 as a double array of point sources.To complement the scalar higher symmetric periodic Green's function (HSPGF), its dyadic counterpart is defined as where the mirroring dyad R z is defined as In the computation of G T , the additional factor, e − j k t00 •ρ g , is added due to the translation of the source points of the top part by where b 1 and b 2 are translation factors in terms of vectors s 1 and s 2 , respectively.For mirror-symmetric structures, b 1 = b 2 = 0, and for glide-symmetric structures, b 1 = b 2 = 0.5.Only a combination of the values of 0.0 and 0.5 is possible, as in other cases, the higher symmetry is broken and the current of the top part of the structure cannot be expressed with a phase shift and attenuation of the current in the bottom part [6].The terms G mn,t are defined as similar to ( 5), but with the distance R mn changed to The choice of summing or subtracting in (9) represents two different types of symmetries.For mirror-symmetric structures, these correspond to perfect magnetic conductor (PMC) (summation) and PEC (subtraction) mirrored modes, but for glide-symmetric structures, we will refer to them as the plus (+) and minus (−) symmetries.

C. Ewald Acceleration
Expressions ( 9) and ( 12) could be used in simulations as they are, but, in practice, their use is prohibited by a very slow convergence [29], [30], [31], [32], [33].In this work, we extend the Ewald method [34], [35] to the concept of HSPGF to allow for efficient computation of the Green's function.
First, we consider the evaluation of ( 10) and (11) as two separate sums.Then, the individual sums are separately split into spatial and spectral sums [35], [36] as where all sums are double sums with summation indices m and n for the spatial sums and p and q for the spectral sums.When numerically computing the Green's function, the series in (17) are truncated.In this work, the truncation limits are determined by the procedure described in [37].For producing the plots in this article, we have not observed values larger than ±3.The Ewald method introduces the splitting parameter E, which controls the convergence properties.For frequencies where the wavelength is larger than the periodicity of the unit cell, the value of E opt = √ π/A, which provides fastest convergence, can be used [38].Otherwise, when the wavelength is smaller than the periodicity of the unit cell, the values of the spectral and spatial parts become similar in magnitude but different in sign.Due to the subtraction of two large numbers and the finite precision arithmetic, this leads to a loss of significant digits.The issue can be alleviated by changing the value of E at the cost of slower convergence.The procedure from [37] is used in this work to choose E which ensures automatically, without the need of setting a threshold value of frequency, a set number of significant digits even at higher frequencies.The values of the summands in (17) are and for (10).For the contribution of the top part of the structure ( 11) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.and Here, is the unit cell area, and k t pq is the transverse phasing wave vector of the p and q indexed mode The cost of computing the Green's function can be then further reduced by exploiting the properties of the complementary error function erfc(a− jb) = erfc(a + jb) (23) and where a and b are real numbers and bar represents conjugation.Using ( 23) and ( 24) reduces the evaluations of the complex complementary error function, which is the most computationally expensive task in the evaluation of the Ewald sums.

III. NUMERICAL SOLUTION A. Discretization of the Equation
Now, we can proceed with numerically solving (1).First, the surface inside the enclosed structure is discretized into triangles.During meshing, it is necessary to ensure that the edge triangle vertices coincide perfectly after a translation by a period.Then, additional triangles are added, as shown in Fig. 2.These triangles must be exact copies of their periodic counterparts as they ensure continuity of current between two adjacent unit cells.Then, we discretize (1) by expanding the surface current density J(r ′ ) into where n are basis functions and I n are their scaling coefficients.Here, Rao-Wilton-Glisson (RWG) basis functions are used, which are defined as linear interpolation functions over triangle pairs with a shared edge [39].Then, a system of equations is obtained with a testing procedure.In this work, the testing functions used on (2) and ( 3) are also RWG basis functions, performing the so-called Galerkin testing.Ultimately, we obtain expressions for and which are combined to form the elements of the where each element is equal to To evaluate the dispersion diagram of the structure under analysis, we assume no incident field (i.e., no impressed field is necessary for modes to exist in the structure), obtaining the following system of equations: where [I ] is an N × 1 vector collecting the coefficients I n of ( 25) and [0] is the N × 1 vector of zeros.The system in ( 29) has a nontrivial solution when the value of the wave vector k t00 at a frequency f is such that the matrix [Z ] is singular.However, due to finite precision of evaluation of integrals in ( 26) and ( 27), the matrix is never truly singular.When the matrix is evaluated at frequency and wave-vector pairs close to the solutions, its determinant goes to zero.The solutions of the system are therefore the set of values of k t00 at f where the determinant is small and at a local minimum.

B. Evaluation of Singular Integrals
An important aspect of evaluating the impedance matrix in ( 29) is a careful treatment of integrals (26) and (27).For basis functions that are overlapping or are adjacent, the integrals become singular and integration by standard Gaussian quadrature rules cannot produce an accurate result.To efficiently and accurately compute singular integrals, a singularity cancellation procedure named radial-angular method has been proposed in [40].This method was developed to work with free-space Green's function, which only contains one singularity.Thus, in the HSPGF, the application of the method cancels the singularity corresponding to 1/R 00 in (6).Since we have a function with multiple singularities, this method must be modified such that the most important singularities are canceled.Note that the method computes interaction between individual cells (here triangles) and not the entire basis functions.The matrix elements Z mn are thus obtained by summing all four individual interactions between the elements of the two basis functions.In the following, we describe how to apply the radial-angular method in the case of the proposed HSPGF to cancel singularities arising from periodic source replicas and quasi-singularities in ( 6) and (15) (i.e., sources from adjacent unit cells) and the quasi-singularities in the case of a small gap between the surfaces arising from the top sources, contained in (11).Here, we describe the two possibilities in the following.
First, the periodicity of the function requires a modification, which ensures that the singularity closest to the observation element is canceled.For example, if the source element, shown in blue in Fig. 3, is on the opposite side of the unit cell than the test element (orange), the periodic image of the source (dashed blue) is the dominant singularity.Thus, the singularity of the periodic image element 1/R 11 is not appropriately canceled as the radial-angular method cancels only the 1/R 00 singularity.To rectify this problem, these element interactions are treated by first shifting the test elements with substitution of where the factors of shifting M 0 and N 0 can be obtained through and In Fig. 3, M 0 = N 0 = −1.After the translation, the dominant singularity becomes 1/R 00 and the radial-angular method can be applied to properly cancel the singularity, as shown in Fig. 3(b).Since shifting one of the elements causes it to be outside of the original unit cell, the result of the integration is different from a multiplicative constant.Then, the integral is backpropagated using the Floquet theorem where Z kl mn here represents the contribution of interaction between elements k and l to the total interaction between the two basis functions.The procedure described here could be adapted to shift the source and not the test element to ensure that 1/R 00 is the dominant singularity.In our case, we choose to shift the test element as it corresponds to the outer loop of the impedance matrix evaluation and is thus easier to implement in our code.
Second, there are cases when a quasi-singular behavior of the top Green's function must be canceled for accurate integration.This is particularly true for structures with a small gap between the top and bottom plates, where singularities from both contributions to Green's function of top and bottom sources are appreciable, such as the example presented in Fig. 4. Here, the blue and dashed blue triangles are the bottom and top source elements, and the test element (orange) is sufficiently close to both such that 1/R 00 and 1/R 00,t need to be canceled.In these cases, the radial-angular method is applied twice, once to the contribution of the bottom part of the structure, using only the Green's function in (10), and once to the contribution of the top part, with the Green's function in (11).For the bottom part, the radial-angular method is applied directly.For the top part, the test element is first mirrored by mirroring the test element in z, as shown in Fig. 4. The integral is then evaluated with (10), which is the Green's function of the bottom part.In this way, the roles of R 00 and R 00,t are exchanged and the radial-angular method cancels the singularity.

C. Zero-Search Algorithm
To find the solutions of (29), we search for all numerical singularities in the lowest eigenvalue, λ min , of the impedance matrix [Z ] within a given frequency region.Here, both plus and minus branches of ( 9) and ( 12) must be considered in the search.To accelerate obtaining the results, we use interpolation of the individual matrix elements to obtain a finer sampling of the impedance matrices.We start to fill the dispersion diagram with the case of real modes (i.e., k t00 is real).The first three steps are shown in Fig 5 .The procedure is given as follows.
1) Compute N init impedance matrices in linearly equispaced frequency points (at cross markers in Fig. 5) in the frequency range of interest.2) Use spline interpolation to obtain a finer sampling of impedance matrices with N interpol equispaced frequency points.3) Compute the smallest eigenvalue λ min of these matrices.
The magnitude of the eigenvalues is presented in full line and the phase is presented in the dashed line in Fig. 5. 4) Find solutions where the phase of λ min crosses zero and the magnitude is a local minimum.5) Use bisection on the phase of λ min to refine the solution up to the desired precision.In the second step, the higher frequency solution converges, so only the lower frequency solution is added.
6) Compute impedance matrices at the found unconverged solutions (plus markers in Fig. 5).7) Reobtain the solution frequencies with the procedure listed in steps from 2) to 5) and with new impedance matrices added to the set of previous initial points N init .8) Repeat steps 6) and 7) until the scheme converges.
The frequency solution is considered to have converged when (34) where we select the pairs of previous solution frequency f previous and next solution frequency f next as the closest two frequencies.In this work, the precision value is set to δ search = 10 −3 .The entire procedure is then repeated for the next value of k t00 .
Once real modes are tracked, we can move to the analysis of complex modes: the procedure is analogous, but a 2-D search is needed.We start by identifying the maximum and minimum frequencies in the computed real part of the dispersion diagram.In many cases, they lie at the edge of the irreducible Brillouin zone (usually , X, or M).In these cases, it is possible to do only a 1-D search as described previously by fixing the real part of components of k t00 according to the considered stopband and finding frequencies of all singular matrices for the desired level of attenuation α.However, when the beginning of a new complex mode is not at an edge of the Brillouin zone, both real and imaginary parts are unknown.In this case, we fix a frequency value close enough to the last frequency of the mode already tracked.Then, we follow the procedure.
1) Compute N init (initial points) impedance matrices in a grid in the search region in the real and imaginary space of k t (cross markers in Figs. 6 and 7).2) Use natural interpolation [41] to obtain finer sampling of impedance matrices in N interpol points.6) Evaluate the impedance matrix in the solution point and repeat steps from 1) to 5) until the scheme converges.First, three steps can be seen in Figs. 6 and 7.Each time, the finer sampling of matrices is obtained with one more initial point, located at the previous solution.7) Move to the next frequency.Estimate the next solution with linear extrapolation and center the search region around it.The convergence criterion to terminate the search is applied to both real real k t,previous − real k t,current real k t,previous < δ search (35) and imaginary parts separately Thus, (35) and (36) need to be satisfied for achieving convergence.

IV. NUMERICAL RESULTS
In this section, we present the results of the modeling procedure described in Sections II and III.Two different geometries, for both mirror-and glide-symmetric unit cell, are considered.Solutions are computed for different levels of mesh refinement to ensure the convergence and stability of the scheme.The obtained results are compared to CST Eigenvalue Solver (CST ES) [42] for propagative modes and to MMTMM [19] for evanescent modes (including the real part of the complex modes).
Note that, differently from the CST ES and the MMTMM methods, singularities in MoM arise when modes are supported either inside or outside of the structure.Using the methodology presented in this article, modes corresponding to solutions in the unbounded pin structure in the external region can appear.However, they can be easily identified as having zero fields inside the bounded region and nonzero fields in the external region.As they are not of interest to this work, these solutions are not shown here.Furthermore, the evanescent modes of the interior region not presented in [19] are also not included for clarity of the plots.

A. Rectangular Holey Structure
First, we consider a rectangular holey structure, whose meshed bottom part is shown in Fig. 8.This geometry is chosen as it consists of flat surfaces only, and hence, a direct subdivision of one triangle to four results in a refined mesh with four times more elements.Such division is applied to the original mesh in Fig. 8(a) once to obtain the mesh of Fig. 8(b) and twice for Fig. 8(c).From the coarsest to the most refined mesh, the number of unknowns N is 177, 708, and 2832.These values correspond to ratios of wavelength to the largest triangle edge of approximately 3.5, 7, and 14 at 60 GHz (the highest considered frequency).The solutions of the first two meshes were computed using N init = 12 starting points in the frequency from 0.5 to 60 GHz and N interpol = 500 interpolation sampling points.The finest mesh refinement was only computed in a 1-GHz region surrounding the coarser solution with N init = 4 and N interpol = 40.
Full dispersion diagrams for the irreducible Brillouin zone obtained with the proposed HSPGF are shown in Fig. 9(a) for the mirror-symmetric unit cell and in Fig. 9(b) for the glide-symmetric unit cell.In both cases, a very good agreement is obtained with respect to the solutions obtained with CST ES, already with the coarsest mesh (cross markers), although it can be observed that the accuracy of the solution slightly deteriorates at the highest frequencies.Refining the mesh once (plus markers) and twice (square markers) improves the result.In the mirror-symmetric case, all modes are PEC-symmetric, while, in the glide-symmetric case, the first mode is minus-symmetric, while the second and third modes are plus-symmetric.Here, the second mode is in fact a Floquet harmonic of the first mode as reported in [19] and could also be obtained by extending the k t00 search space of the first mode.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

B. Circular Holey Structure
For the second structure, we aim to verify the applicability of the code to find complex and evanescent modes.For this purpose, the geometry was chosen to compare the results in [19,Fig. 8].Two meshes were generated to have a ratio of wavelength to largest triangle edge of approximately 4 and 8 at 100 GHz, corresponding to 255 and 831 basis functions, as shown in Fig. 10(a) and (b).Here, due to the curved surface, the entire structure is remeshed instead of splitting each triangle into four.For the generation of the dispersion diagrams, the frequency range considered is 0.5-100 GHz with N init = 20 initial points in the interpolation scheme.For evanescent modes, the frequency ranges are set to coincide with the stopbands at , X, and M. The number of initial points N init varies depending on the width of the stopband such that the impedance matrices are computed every 5 GHz and the number of interpolation points is N interpol = 500.
The full dispersion diagram for the mirror-symmetric unit cell can be seen in Fig. 11.The attenuation in the stopband is shown in Fig. 12(a) and (b) for the cases of X and M regions, respectively.Overall, the results agree very well with [19], with a reduction of accuracy as the frequency increases for the coarsest mesh.Here, the difference is not only due to the larger relative size of the unit cells when compared to the wavelength but also due to approximation of the curvilinear surface with planar triangular elements.Only PEC-mirrored modes are found in both original (cross markers) and refined (plus markers) meshes, which is expected as the gap between top and bottom plates is small enough to suppress any PMC-mirrored modes.
The dispersion diagram for the glide-symmetric unit cell can be seen in Fig. 13.The attenuation in the stopband is presented in Fig. 14(a) and (b) for the cases of X and M , respectively.Four modes are found, two minus-symmetric (yellow and green markers) and two plus-symmetric (orange and purple markers).As before, the original mesh is presented with cross markers and the refined mesh is presented with plus markers.It is worth noting that, due to the extra information of the mode symmetry, utilizing the HSPGF helps with interpretation of the mode evolution.For example, the plus-mirrored evanescent mode in Fig. 14(b) is found at M for the lower frequencies and transitions to a propagative mode at about 30 GHz.Then, with increasing frequency, the mode gradually moves from M to .At , it becomes again evanescent and does not transition into a propagative mode in the frequency range considered.

C. Computational Time Analysis
The computational time to evaluate the impedance matrices is presented in Table I.The structure studied is the rectangular holey unit cell from Fig. 8 in both mirror-and glide-symmetric configurations.Times for both FSPGF and the HSPGF proposed in this article are presented.All computations were done on an HP ProBook 430 G8 Notebook with Intel i7-1165G7 and 16 GB of RAM.Note that the times presented are obtained with a single-threaded program.As the computational time  has slight frequency dependence, the numbers are the average time of ten evaluations of linearly spaced frequencies from 0.5 to 60 GHz and rounded to 1 s.For a fair comparison, the mesh of the full structure (corresponding to the use of FSPGF) is obtained by a mirror and a translation of ( 14) and thus has twice the number of basis functions.It can be observed from Table I that using the HSPGF results in twice faster evaluation of the matrix.While the matrix is reduced to one-fourth number of entries, introducing the top part of the structure in HSPGF results in about double the computational effort in evaluating the Green's function.The average times necessary to find all the solutions (frequencies) for a given value of k t00 in Fig. 9 is presented in Table II.The difference in solution times when comparing the mirrorand glide-symmetric structures can be attributed to a larger computational load of the zero-search algorithm since the first and the second mode are not separated in the XM region.The simulation time for obtaining the entire dispersion diagram in Figs. 9, 11, and 13 with CST ES is from 7 to 9 min for 30 values of k t00 .

V. CONCLUSION AND PERSPECTIVES
In this work, a novel Green's function with a modeling approach to obtain properties of modes of mirror-and glide-symmetric structures has been presented.Using the proposed Green's function with the MoM, we obtained accurate results for both propagative and evanescent modes when compared to results from CST ES and MMTMM [19].Moreover, the Green's function evaluation accelerated via Ewald and the adapted use of the radial-angular method [40] to correctly cancel singular integrals were presented.The solutions were found via a technique for the search of singularities of a matrix function.Using the proposed Green's function was found to reduce the computational time of one matrix evaluation to one half.Furthermore, the Green's function provides additional information on the mode parity, which was found helpful in understanding and tracking the evolution of modes.The presented method can be easily integrated into an existing MoM code and is able to obtain information not directly available from commercial solvers, such as attenuation in the stopband.In the future work, we aim to extend the code to work with dielectric materials and utilize it to study hexagonal unit cells.

Fig. 2 .
Fig. 2. Meshing process.The inside of the structure is meshed and additional triangles (orange) are added by translating their periodic counterparts (blue) in the original unit cell.

Fig. 3 .
Fig. 3. Translation of the test element (orange) for proper cancellation of periodic singular elements.(a) Before translation, the periodic image element indexed 11 (blue dashed) is closest to the test element.(b) After translation, the source element indexed 00 (blue) becomes the closest.

Fig. 4 .
Fig. 4. Mirroring of the test element (orange) for proper cancellation of singular elements from G T .The top mesh is faded as it is present through the HSPGF.(a) Before mirroring, applying radial-angular method cancels the bottom singularity (full blue).(b) After mirroring, the top singularity can be canceled.

Fig. 5 .
Fig.5.algorithm for N init = 12 and N interpol = 600.From top to bottom, found singularities are added to the interpolation scheme.In the second step, the higher frequency solution converges, so only the lower frequency solution is added.

Fig. 6 .
Fig. 6.Magnitude in first three steps of the complex zero-search algorithm with N init = 36 and N interpol = 1600.

Fig. 7 .
Fig. 7. Phase in first three steps of the complex zero-search algorithm with N init = 36 and N interpol = 1600.

Fig. 8 .
Fig. 8. Meshes of the rectangular holey structure.From left to right, they are (a) original, (b) refined 1, and (c) refined 2 meshes.Here, z min = −1.25 mm and z max = −0.25 mm.The parameters periodicity is p = 4 mm, the depth of hole is h = 1.5 mm, the width of hole is w = 3 mm, and the gap is 0.5 mm.

Fig. 9 .
Fig. 9. Dispersion diagrams for (a) mirror-and (b) glide-symmetric rectangular holey structures with the geometry simulated and the irreducible Brillouin zone depicted in the insets.

Fig. 10 .
Fig. 10.Rectangular holey structure meshes for (a) original mesh and (b) refined mesh.Here, z max = −0.025mm and z min = −1.025mm.The dimensions marked are periodicity p = 3.2 mm, hole diameter d = 2.56 mm, hole depth h = 1 mm, and the gap between the top and bottom plates is 0.05 mm.The meshing software used is Gmsh [43].

Fig. 11 .
Fig. 11.Dispersion diagram for the irreducible Brillouin zone and the mirror-symmetric unit cell, both depicted in the insets.

Fig. 12 .
Fig. 12.(a) Attenuation in the X region for the mirror-symmetric unit cell depicted in the inset.The solutions from 30 to 60 GHz were found at X, from 76 to 85 GHz at and the mode from 85 to 89 GHz is a complex mode with the real part shown in Fig. 11.(b) Attenuation in the M region for the mirror-symmetric unit cell depicted in the inset.The solutions from 49 to 60 GHz were found at M and from 75 to 100 GHz at .

Fig. 13 .
Fig. 13.Dispersion diagram for the irreducible Brillouin zone and the glide-symmetric unit cell, both depicted in the insets.

Fig. 14 .
Fig. 14.(a) Attenuation in the X region for the glide-symmetric unit cell depicted in the inset.The solutions from 46 to 59 GHz are located at and the solutions from 59 to 87 GHz are a complex mode with the real part presented in Fig. 13.(b) Attenuation in the M region for the mirror-symmetric unit cell depicted in the inset.The solutions from 49 to 60 GHz were found at M and the solutions from 75 to 100 GHz at .

TABLE I AVERAGE
COMPUTATIONAL TIME (IN SECONDS) FOR ONE IMPEDANCE MATRIX EVALUATIONTABLE II AVERAGE COMPUTATIONAL TIME (IN MINUTES) FOR OBTAINING THE SOLUTIONS (FREQUENCIES) IN FIG. 9 FOR FIXED k t00