A symmetry-based multimodal transfer-matrix method for the analysis of 2D-periodic structures

IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES 1 A Symmetry-Based Multimodal Transfer-Matrix Method for the Analysis of 2D-Periodic Structures Jesus M. Jimenez-Suarez , Graduate Student Member, IEEE, Francisco Mesa , Fellow, IEEE, and Oscar Quevedo-Teruel , Fellow, IEEE Abstract—We propose a systematic and efficient extension of the multimodal transfer-matrix method to obtain the dispersion diagram of structures with 2-D periodicity specifically targeted to primitive unit cells that possess internal symmetries. When symmetry planes can be applied, the study of the unit cell can be simplified to a number of 1D-periodic scenarios that depend on the boundary conditions imposed by the symmetry planes. The study of these 1D-periodic scenarios is simpler, more accurate, and requires less computational cost. The proposed methodology has been validated with different examples of periodic structures with different lattices (squared, rectangular, and hexagonal), symmetries, and motifs. Furthermore, this approach brings about a deeper understanding of the study of the Brillouin zone (BZ) and the relationship between phase shift and paths on its irreducible Brillouin zone (IBZ). Index Terms—Dispersion analysis, hexagonal lattice, multimodal analysis, periodic structure, scattering matrix, symmetry planes. I. INTRODUCTION PERIODIC structures are commonly used in microwave engineering and optical devices to control the propagation and scattering of electromagnetic waves [1],[2],[3],[4]. Modifying the elements and/or geometric characteristics of the unit cell of periodic structures enables these structures to have various functionalities that provide great versatility and applicability on different devices, such as filters, antennas, frequency-selective surfaces, and lenses [5],[6], [7],[8],[9]. Among periodic structures, those with hexagonal lattices have gained the attention of researchers due to the improvement of their electromagnetic properties relative to the square/rectangular lattice case [10],[11],[12]. In addition, the hexagonal lattice allows for a specific higher symmetry known as the mirror half-turn [13]. This higher symmetry has been shown to lead to higher refractive indices and wider stopbands [14]. The characteristics of periodic Received 26 December 2024; revised 25 February 2025; accepted 17 March 2025. The work of Jesus M. Jimenez-Suarez was supported by the Horizon Europe Research and Innovation Program and UKRI through the GENIUS Project, Marie Sklodowska-Curie Grant, under Agreement 101072560. The work of Francisco Mesa was supported in part by MICIU/AEI/10.13039/501100011033 under Grant PID2023-148281NB-I00 and in part by ERDF/EU. (Corresponding author: Oscar Quevedo-Teruel.) Jesus M. Jimenez-Suarez and Oscar Quevedo-Teruel are with the Division of Electromagnetic Engineering and Fusion Science, KTH Royal Institute of Technology, 100 44 Stockholm, Sweden (e-mail: [email protected]; [email protected]). Francisco Mesa is with the Department of Applied Physics 1, ETS Ingenier´ ıa Inform´ atica, Universidad de Sevilla, 41012 Seville, Spain (e-mail: [email protected]). Digital Object Identifier 10.1109/TMTT.2025.3554934. structures are typically represented by dispersion diagrams, which illustrate the relationship between wavevector components and frequency, providing insight into the pass/stopband structure, phase/group velocities, and attenuation levels [2], [3],[15]. The calculation of these dispersion diagrams is performed using analytical/numerical methods, such as equivalent circuits [16],[17], mode matching [18],[19], method of moments [14],[20],[21],[22], finite-difference timedomain modeling [23], or by means of commercial full-wave solvers. In recent years, some of the authors have successfully employed a hybrid approach, the multimodal transfer-matrix method (MMTMM), for this task [24],[25],[26], as well as for obtaining the effective constitutive material parameters [27]. The MMTMM combines and takes advantage of the strengths of full-wave commercial software to deal with the scattering parameters of general intricate structures and an ad hoc postprocessing to study the structure in a periodic environment [26]. Symmetries are crucial for solving electromagnetic problems because they enable us to reduce the complexity of the entire problem to a fraction. Furthermore, symmetry analysis provides deeper insight into mode classification and modal electromagnetic properties [28],[29],[30],[31]. This has been used, for example, to decompose the incident field (even and odd excitations), to reduce problem computation, and to simplify the formulation of numerical methods [31],[32], [33],[34],[35],[36],[37]. In the study of periodic structures, the first obvious simplification comes from the periodicity of the problem, and therefore, only the unit cell needs to be studied [2],[3]. However, if the unit cell is symmetric, its study can be further reduced by taking advantage of the use of symmetry operations within it. In this work, an enhanced version of the MMTMM is proposed for calculating the dispersion diagram of symmetric 2D-periodic structures. In particular, the MMTMM can be considerably simplified by applying symmetry planes within the unit cell so that only 1D-periodic situations have to be considered. These 1D-periodic scenarios depend on the boundary condition imposed by the type of symmetry plane (electric/magnetic) applied. This article is organized as follows. In Section II, symmetric unit cells are characterized, highlighting the symmetry planes that are used in this work. Section III presents an outline of the MMTMM to calculate dispersion diagrams. A brief description of the analysis of 2D-periodic structures is reported in Section III-A, while the proposed approach using symmetry ©2025 The Authors. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 2 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES Fig. 1. Examples of symmetric unit cells. (a) Square unit cell with a rectangular motif and twofold axial symmetry. (b) Square unit cell with a circular hole and fourfold axial symmetry. (c) Rectangular unit cell with an arrow-point motif and onefold axial symmetry. (d) Rectangular unit cell with a cross motif and twofold axial symmetry. (e) Hexagonal unit cell with a rectangular motif and twofold axial symmetry. (f) Hexagonal unit cell with a circular motif and sixfold axial symmetry. planes is given in Section III-B. The numerical results of the analysis of square/rectangular unit cells are exposed in Section IV and hexagonal unit cells are treated in Section V. Finally, the main conclusion is discussed in Section VI. II. CHARACTERIZATION OF UNIT CELLS Earlier research on waveguides has often taken advantage of symmetries to reduce the computational cost of modal analysis in these structures using numerical techniques such as the finite element method and mode matching [31],[38],[39]. A similar procedure can be used to analyze symmetric unit cells of 2D-periodic structures and leverage symmetries to reduce the computational effort to calculate dispersion diagrams [14], [22],[37]. The initial task is to find the symmetries displayed by the unit cell. A symmetry operation within an object is a spatial transformation that preserves the shape of the object. In the case of a 2-D unit cell located, say, in the xy plane, the symmetry operations can be categorized into two types: rotational symmetries about an axis aligned with the z-axis and reflective symmetries in planes perpendicular to the xy plane [40]. A unit cell exhibiting splanes of symmetry will allow us to reduce the study region of this unit cell to 1/2sof its original size. A series of examples with some common shapes of a planar unit cell with its motif are illustrated in Fig. 1, highlighting their symmetry planes with dashed lines. Considering that MMTMM is a hybrid method combining commercial software and an ad hoc algorithm, in this work, we will focus only on the symmetry planes perpendicular to the Cartesian axes that are allowed for most solvers and mesh types in commercial software, specifically those highlighted by blue dashed lines in Fig. 1. The region enclosed by these planes, colored light blue in the figure, defines the minimal study region to be considered for analysis. These symmetry planes correlate with the corresponding even and odd modes of the structure [28]. They also play a significant role in the analysis of wave propagation within periodic structures, as they are essential to define the irreducible Brillouin Fig. 2. Equivalent network of the unit cell of a 2D-periodic structure with periodicity along the xand y-directions. Nmodes are retained on each of the four geometrical faces. zone (IBZ) [3],[13],[15]. Although Fig. 1only shows the symmetry planes within the xy plane of the unit cells, an additional possible symmetry plane can exist perpendicular to the z-axis if the unit cell exhibits mirror symmetry in the vertical direction. III. ANALYSIS OF PERIODIC STRUCTURES A. MMTMM for 2D-Periodic Structures Assuming a time-harmonic regime characterized by an angular frequency ω, the study of structures with 2-D periodicity by means of the MMTMM can be carried out by treating their unit cells, for example in a square/rectangular lattice case, as a 4N-port network with four terminal planes located on the side faces of the unit cell [26]. As shown in Fig. 2, ports in sides 1 and 2 are considered input ports, and ports in sides 3 and 4 as output ports. In each side, we consider Nports associated with their corresponding first Nsignificant modes. A general-purpose full-wave commercial software is first used to obtain the frequency-dependent multimodal scattering matrix of the unit cell, from which the 4N×4Nmultimode transfer (ABCD) matrix [T(ω)] that characterizes the unit cell in the periodic environment can be easily obtained. Using this matrix and Bloch’s theorem, the following generalized eigenvalue problem can be found to describe the wave propagation within the unit cell: [T]ψ=λψ (1) where [T] is the transfer matrix, λis an eigenvalue, and ψis the associated eigenvector, representing the voltages and currents at the port of the unit cell. For a 2D-periodic structure with periodicity along xand y-directions, this eigenvalue problem can be written as follows [26]: T(ω)2 6 6 4 V1 V2 I1 I2 3 7 7 5 =2 6 6 4 V3 V4 I3 I4 3 7 7 5 =2 6 6 4 e−jkxpxV1 e−jkypyV2 e−jkxpxI1 e−jkypyI2 3 7 7 5 (2) This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. JIMENEZ-SUAREZ et al.: SYMMETRY-BASED MMTMM FOR ANALYSIS OF 2D-PERIODIC STRUCTURES 3 where Vnand Inare the voltage and current arrays related to the Nport modes in side n,pνis the period, and kν=βν−jαν is the wavenumber along the ν-direction (ν≡x,y), with βνand ανbeing their corresponding phase and attenuation constants. Since the previous eigenvalue problem does not conform to a typical linear form, solving it is not straightforward, requiring the application of a zero search algorithm [41], a troublesome numerical task, especially if searching for complex-value solutions. This problem was overcome by the linearization presented in [42], which using a permutation matrix makes it possible to rearrange the original eigenvalue problem into an equivalent linear eigenvalue problem when the phase shift condition is fixed in one direction. It allows us to solve the problem using standard packages of computational algebra available for any programming language. Despite this convenient linearization procedure, the numerical solution of the resulting eigenproblems can still be cumbersome showing spurious solutions and numerical noise. Some of these problems can be caused by the required numerical postprocessing required to setup the corresponding eigenproblem from the multimodal scattering matrix of the unit cell. Small errors in obtaining this scattering matrix by means of the generalpurpose electromagnetic simulator are likely to be magnified during its transformation into a transfer matrix and by the subsequent numerical operations needed to solve the eigenproblem. In our experience, another common numerical issue occurs due to the poor condition number of the 4N×4N transfer matrix in the initial 2-D eigenproblem when N&6, particularly if one of the selected input/output port modes is hardly relevant for characterizing the Floquet mode within the intended frequency range of the unit cell. In this case, it might instead introduce numerical noise in the computation of the solution. This issue will be discussed further when dealing with numerical results. B. Application of Symmetry Planes in the MMTMM As well reported in the literature [29],[30],[32],[33],[34], [35],[36],[37],[38], the presence of symmetry planes in the structure under analysis helps to simplify the corresponding electromagnetic problem. In our case, the symmetry planes within the unit cell introduce important simplifications in the eigenvalue problem discussed in Section III-A. If a unit cell contains a plane of symmetry, then any mode must exhibit symmetry or antisymmetry (i.e., even or odd characteristics) relative to this plane [2],[4],[34],[35],[37]. This requirement does not depend on the unit cell’s shape, design, or materials and implies the existence of a perfect magnetic conductor (PMC) or a perfect electric conductor (PEC) in this plane of symmetry. Enforcing these PEC/PMC planes means that waves must propagate parallel to these planes or be reflected by them. The main implication of this is that the corresponding electromagnetic propagation eigenproblem within the unit cell simplifies to a 1D-periodic problem in a specific direction, with its eigenvalues only depending on the modal wavenumbers in that direction. In this situation, the periodic structure is analyzed as the 2N-port network shown in Fig. 3with its Fig. 3. Equivalent network of the 1D-periodic scenario resulting from the application of 1SP, with Nmodes are retained on each of the two geometrical faces. Fig. 4. Examples of two symmetry planes in a square unit cell with a mirrored circular hole. These planes are identified by following the nomenclature proposed in this article. corresponding eigenvalue problem given by T(ω)V1 I1=V2 I2 =e−jkνpνV1 I1.(3) Dealing with 1D-periodic cases provides distinct advantages compared with 2D-periodic scenarios, as the involved matrices are smaller in size and exhibit considerably lower condition numbers. Consequently, it makes the resulting problem less computationally intensive, exhibits reduced numerical noise, and is more precise. Significantly, breaking down the original 2D-periodic problem into several 1D-periodic scenarios (see [24] for further information on 1-D cases) eliminates the need for a zero search in the complex plane and implicitly accomplishes the MMTMM linearization reported in [42]. In the full-wave simulator used to obtain the generalized scattering matrix of the unit cell of the original 2D-periodic problem, the symmetry planes are regarded as PMC/PEC boundary conditions. The nomenclature proposed to identify these planes is Aν d, where A identifies the type of boundary condition in the plane (M for PMC and E for PEC), νidentifies the axis perpendicular to the plane, and dis the position on the axis νwhere the plane is placed. For example, the planes colored in Fig. 4are denoted by following the above nomenclature as Ex 0, My 0, and Ez 0. The characteristic (PMC/PEC) chosen for the symmetry plane is directly related to which path is considered in the IBZ. Even/odd modes are excited when PMC/PEC are imposed, since even/odd modes correspond to 2nπ/(2n+1)πphase shifts [n=0,1,2, . . .] between the original boundaries of the primitive unit cell. To clarify this fact, the square unit This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 4 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES Fig. 5. (a) Square unit cell with a mirrored circular hole as motif and (b) first Brillouin and IBZ of the structure. (c) 2-D geometric configuration of the periodic structure whose unit cell is illustrated in (a). A rotated supercell is highlighted in pink solid line. (d) First BZ for the supercell is colored in pink. TABLE I PHASE-SHIFT CONDITIONS FOR THE PATHS IN A TRIANGULAR IBZ AS IN SHOWN IN FIG.5(B) cell shown in Fig. 5(a) and its IBZ will be simplified to different 1D-periodic scenarios. The IBZ for this structure is the region within the triangle edge ΓXMΓshown in Fig. 5(b), which can be studied under the phase shift condition given in Table I. If we impose My 0and My p/2planes to the unit cell, wave propagation is enforced along the x-direction with a zero phase shift between the faces y=−p/2 and y=p/2, as demonstrated in Appendix A. This scenario corresponds to study the path ΓX in the reciprocal space. If an Ey 0was used instead of the previous My 0, a phase delay of πin the y-direction would be imposed, which means solving the path YM in Fig. 5(b). Due to the symmetries in the structure under study, this path is equivalent to the path XM belonging to the IBZ. Thus, to study propagation along the xor y-direction, one can focus on a single direction and utilize symmetry planes orthogonal to the other. Finally, propagation along the path ΓM can be addressed as an equivalent 1D-periodic scenario by using symmetry planes within the rotated square supercell of size √2p, as illustrated in Fig. 5(c). Using symmetry planes and analyzing the propagation along the x0-axis within the supercell, it effectively parallels the solution of the path ΓM with respect to the primitive cell (PC) [see Fig. 5(d)]. In the previous example, we have reduced the original simulation region (the complete unit cell) to half its size. We can further reduce this half-size simulation region to only a quarter of the unit cell if there is another symmetry plane perpendicular to the direction of wave propagation under study. The application of this second symmetry plane is equivalent to having used even/odd excitations; that is, to having imposed PMC/PEC on this plane. This situation will be identified with the additional label “e/o” in the figure legends. As well reported in the literature [4],[28], in this case, the generalized scattering matrix can be obtained from the reflection matrices associated with even/odd excitations as follows: [S11]=1 2Se 11+So 11 (4) [S21]=1 2Se 11−So 11.(5) In the unit cell shown in Fig. 5(a), there is also another plane of symmetry perpendicular to the z-axis in the middle of the unit cell. This additional symmetry would allow us to reduce the simulation region to one-eighth of the original unit cell. In this study, we limit this symmetry plane to being electric, as our attention is on wave transmission in the unit cell similar to a perturbed parallel plate waveguide (PPW). Furthermore, if the height of the PPW (air gap) was less than λ/4 and the symmetry plane was magnetic, the waves could not propagate in the structure. As demonstrated, an apparent benefit of employing symmetry planes is to reduce the computational cost, since the simulation domain is minimized to 1/23of its initial volume. In addition, implementing symmetry planes ensures that only modal solutions compatible with the boundary conditions imposed in the simulation region are permitted. Consequently, modal solutions that exhibit a particular symmetry can only be computed from input/output modes that possess the same symmetry. An important consequence of this fact when applying the MMTMM is the significant decrease in the number of input/output modes required to solve the eigenproblem, which strongly improves the numerical efficiency of the method. In Sections IV and V, several examples of periodic structures with different shapes, motifs, and symmetries are studied to demonstrate the potential of using symmetries and the advantages of its use in the study of periodic structures. IV. SQUARE/RECTANGULAR UNIT CELL A. Square Unit Cell With a Circular Mirrored Hole The first periodic structure under study is the square unit cell with a circular mirrored hole shown in Fig. 5(a). In this and subsequent examples presented in this work, an air gap of g=0.05 mm has been employed. Numerous studies highlight the significance of this parameter for the system’s overall performance This unit cell has three symmetry planes as explained previously, and its IBZ is the colored-blue triangle region shown in Fig. 5(b). Following the guidelines discussed in Section III-B, Fig. 6(a) and (c) shows the phase shifts obtained for the path ΓX and XM, respectively. In these This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. JIMENEZ-SUAREZ et al.: SYMMETRY-BASED MMTMM FOR ANALYSIS OF 2D-PERIODIC STRUCTURES 5 Fig. 6. Dispersion diagrams of the 2D-periodic structure shown in Fig. 5(a) with dimensions (in milimeters) are p=3, h=1, d=2.4, and g=0.05. (a) Phase shift and (b) normalized attenuation constant (k0is the free-space wavenumber) along the path ΓX, and (c)–(d) along the path XM. Note that the label “e/o” means that even/odd excitations have been applied due to the existence of another symmetry plane along the propagation direction. figures, we can observe a good agreement between the phaseshift results provided by the CST Eigenmode Solver (ES) and those given by the MMTMM without using symmetry planes (w/o SP). In these and the following figures, MMTMM results using symmetry planes will be identified by the specific combination or the number of symmetry planes used. In Fig. 6(a) and (c), the results given by the MMTMM using symmetry planes also show an excellent agreement with those calculated with the MMTMM w/o SP, thus validating the proposed strategy of using symmetry planes. As discussed in [26] and [42], the CST ES does not directly provide solutions for the attenuation constant, and thus, only the results calculated with the MMTMM for the corresponding normalized attenuation constants are shown in Fig. 6(b) and (d). As reported in previous works [14],[26], the values of the attenuation constant provided by MMTMM were validated with in-house MoM codes and were found to converge with the same number of modes required for the convergence of the phase constant. The color of the dashed lines used for the attenuation constant in Fig. 6(b) and (d) corresponds to the color of the symbols used for the phase shift in Fig. 6(a) and (c). The same correspondence will be used in the following figures. In Fig. 6(a) and (b), we can see some spurious solutions that appear when using the MMTMM w/o SP (those square symbols scattered randomly in the figures). These spurious solutions are attributed to a numerical coupling between the modes that excite the structure in the postprocessing of the transfer matrix. This unexpected coupling is prevented by the use of symmetry planes, since only the modes compatible with the boundary conditions imposed by symmetry planes are excited. Furthermore, the symmetry planes reduce the computational cost not only because of the smaller size of the region to be simulated but also because of the fewer Fig. 7. Condition number of (a) scattering and (b) ABCD matrices used in the computation of MMTMM for the unit cell is illustrated in Fig. 5(a). Symmetry planes (SPs) refers to the combination Mx 0Mx p/2. modes needed for achieving satisfactory solution convergence. This reduction in the number of modes results from smaller waveguide port size and the further limitation on the number of modes compatible with the imposed boundary conditions of the symmetry plane. For the path ΓX, using MMTMM w/o SP requires six modes at each of the four ports in the initial 2D-periodic setup, resulting in a 24-sized transfer matrix. In contrast, incorporating symmetry planes in the MMTMM approach reduces the requirement to only two modes for input and output ports in the simplified 1D-periodic case, shrinking the transfer-matrix size to 4. Although size reduction is the same for the scattering and transfer matrices, their influence on their condition numbers varies between these two matrices, as illustrated in Fig. 7for the combination Mx 0Mx p/2. In this case, Fig. 7(a) shows that the scattering matrix condition number is reduced by a factor of 6 when one or two symmetry planes are imposed, and this condition number is further reduced to 1 when three symmetry planes are applied (the problem becomes monomodal in this case). When the scattering matrix is transformed into the transfer ABCD matrix (which is the fundamental matrix required for solving the eigenvalue problem that characterizes the periodic problem), Fig. 7(b) shows that the use of symmetry planes makes the condition number decrease by several orders of magnitude and smoother (with no peaks). A high condition number leads to an eigenvalue problem more sensitive to numerical noise; hence, small variations in the input matrix result in the appearance of spurious modes in the final solution. Some of these spurious samples are shown in Fig. 6(a) and (b) and always appear when applying the MMTMM w/o SP, although they will not be explicitly shown in the following examples. This situation concerning condition numbers also appears in other situations treated in this study and will not be further elaborated on. In general, we can find situations in which more than one combination of symmetry planes is necessary to compute all the modes present in the dispersion diagrams. For example, the first and third modes along the path XM in Fig. 6(c) and (d) are obtained when the combination Ey 0My p/2is applied while the second mode requires the combination My 0Ey p/2. Occasionally, flat modes appear in the dispersion diagram, such as the one shown in Fig. 6(c). In our experience, these modes are very difficult to compute by the 2-D implementation of the MMTMM w/o SP. However, the use of symmetry planes This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 6 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES Fig. 8. (a) Phase shift and (b) normalized attenuation constant along the ΓM path for the structure shown in Fig. 5(a). A supercell has been used to compute the dispersion diagram using symmetry planes. imposes convenient boundary conditions in the problem which allow us to find this type of mode. The convenience of using a square supercell to calculate the dispersion diagram along the path ΓM within a 1D-periodic framework was previously discussed in Section III-B. This strategy significantly affects how to obtain the dispersion diagram as shown in Fig. 8. Note that the phase and attenuation constant for propagation along this path (the 45◦direction) need to be scaled by a factor of √2 to make the phase shift vary in the range [0, π]. The use of a supercell for this purpose was previously analyzed in [13] by some of the authors, where the emergence and identification of additional modes were examined. In the present case, each modal solution calculated using the PC appears as two different solutions when the supercell is simulated. This occurs as a consequence of the compression of the Brillouin zone (BZ) of the supercell. Specifically, the size of the supercell BZ, colored pink in Fig. 5(d), is half that of the BZ of the PC, colored blue in the same figure. Following the discussion in [13], the extra mode is readily identified and the solution of the supercell is easily transformed into the one expected from the PC. The simplest scenario with the smallest simulation region appears when three symmetry planes are applied to the unit cell shown in Fig. 5(a). A very significant feature often found in this case is that the dispersion behavior of the fundamental modes of the structure can be obtained without solving an eigenvalue problem. In this situation, the symmetries of the modal solution in the reduced simulation region are compatible only with a single mode on the input port. This makes the problem become monomodal and, therefore, can be solved by the following equation [2],[4]: cos (kp)=A(ω)+D(ω) 2(6) where kis the complex wavenumber and A(ω) and D(ω) are elements of the 2 ×2 frequency-dependent monomodal transfer matrix. A comparison between the results obtained from CST and MMTMM with and without the use of (three) symmetry planes is depicted in Fig. 9, showing excellent agreement between all of them. The comparison of the computational time required by CST and MMTMM with and without the use of symmetry planes is shown in Table II. In this comparison, 200 samples have been set up for CST Eigensolver and 400 for MMTMM. The computational cost of the MMTMM includes the time required to obtain the initial Fig. 9. Dispersion diagram along the path ΓXMΓusing three symmetry planes of the 2D-periodic structure shown in Fig. 5(a). Dimensions (in millimeters) are p=3, h=1, d=2.4, and g=0.05. TABLE II COMPUTATION TIME TO OBTAIN THE RESULTS DEPICTED IN FIG.9 scattering matrices using the CST Frequency-Domain Solver (CST-FDS) and the resolution time of the eigenvalue problem. The resolution of the eigenproblem is implemented as a postprocessing step using standard algebraic packages available in MATLAB, and its computation time is negligible compared to that required by CST. The computational effort required for a given number of symmetry planes (# of SPs) does account for all necessary 1D-periodic scenarios (automating the interaction between MATLAB and CST is here recommended). This results in a clear reduction in computational cost when using an increasing number of symmetry planes. All calculations were performed using a desktop computer with an Intel i78700 and 64 GB of RAM. B. Square Unit Cell With an Elliptical Mirrored Hole The proposed methodology is now applied to the structure illustrated in Fig. 10(a), where an elliptical hole substituted the circular hole of the structure analyzed in Section IV-A. The geometry of the structure makes the IBZ no longer a triangle but a square region, as shown in Fig. 10(b). Due to the symmetries of the structure, the original 2D-periodic structure can be reduced to four distinct 1D-periodic cases, comprising two for the x-direction and two for the y-direction. The phaseshift conditions for the IBZ edge are given in Table III. In the present case, several combinations of symmetry planes are again necessary to calculate all the modes present in the dispersion diagram. For example, on the path ΓY, the two first modes are even, so they are obtained when the combination My 0My p/2is used, and the third odd mode when Ey 0Ey p/2is applied. It becomes apparent that the use of symmetry planes provides a good physical insight into the nature of the modes that propagate in the structure. The dispersion diagram along the contour of the IBZ is depicted in Fig. 10(c), which also This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. JIMENEZ-SUAREZ et al.: SYMMETRY-BASED MMTMM FOR ANALYSIS OF 2D-PERIODIC STRUCTURES 7 Fig. 10. (a) Square unit cell with a mirrored elliptical hole as motif with dimensions (in millimeters): p=3, h=1, dx=1.6, dy=2.6, and g=0.05. (b) First Brillouin and IBZ of the structure. (c) Dispersion diagram along the path ΓXMYΓ. TABLE III PHASE-SHIFT CONDITIONS FOR THE IBZ BOUNDARY IN FIG.10(B) exhibits a very good agreement between the CST results and the MMTMM with and without the use of symmetry planes. C. Glide-Symmetric Square Unit Cell With Rectangular Holes As a final example of a square unit cell, we study a periodic structure with higher symmetry; specifically, the square unit cell with rectangular glide-symmetric holes shown in the inset of Fig. 11(b). Observe that due to the glide symmetry, the plane z=0 no longer behaves as a mirror plane for the structure and, hence, only two symmetry planes perpendicular to the Cartesian axes can be implemented. The IBZ for this particular structure retains a square geometry, as shown in Fig. 10(b). Consequently, it can be analyzed by enforcing the same vertical symmetry planes as those used in the square unit cell with a mirrored elliptical hole. Thus, only the results for the phase shift and the attenuation constant along one of the IBZ paths are plotted in Fig. 11. Again, there is good agreement between the results obtained from CST and MMTMM, with and without the use of symmetry planes. Fig. 11. (a) Phase shift and (b) normalized attenuation constant along the path ΓX for the unit cell shown in (b). Dimensions (in millimeters) are p=3, a=1.3, b=2.1, h=1, and g=0.05. To understand the attenuation constants depicted in Fig. 11(b), it is convenient to first examine the phase shifts shown in Fig. 11(a). This figure reveals that the stopband is bounded by an even (first) mode and an odd (second) mode. Consequently, the attenuation constant for the even mode remains at 0 until the stopband, where it begins to increase. On the other hand, the odd mode is evanescent up to approximately 70 GHz, with its attenuation constant progressively decreasing to 0 at the cutofffrequency where propagation begins. These two attenuation constant curves do not couple each other since they correspond to modes of different natures. D. Rectangular Unit Cell With a Mirrored Circular Hole Next, we study a rectangular unit cell with a circular mirrored hole, as shown in Fig. 12(a). In this case, the IBZ is the rectangle shown in Fig. 12(b), which can be studied by imposing the same phase shift conditions as in the previous case shown in Table III. In the path ΓX, the first, second, and fourth modes, as well as the complex mode between the third and fourth modes, are calculated by applying the combination My 0My p/2. The third mode needs the combination Ey 0Ey p/2due to its odd nature. The combinations Ey 0My p/2and My 0Ey p/2are necessary to obtain all the modes present in the path XM. A similar procedure, but imposing symmetry planes perpendicular to x-direction, is used in paths MY and ΓY. The dispersion diagram computed for this structure is plotted in Fig. 12(c), again obtaining a very good agreement between CST and MMTMM with and without using symmetry planes. This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 8 IEEE TRANSACTIONS ON MICROWAVE THEORY AND TECHNIQUES Fig. 12. (a) Rectangular unit cell with a mirrored circular hole as motif. (b) First Brillouin and IBZ of this structure. (c) Dispersion diagram along the path ΓXMYΓ; dimensions (in millimeters): px=3, py=3.7, h=1, d=2.4, and g=0.05. V. HEXAGONAL UNIT CELL An additional and final application of the proposed methodology pertains to 2D-periodic structures with a hexagonal lattice. In contrast to square unit cells, the analysis of hexagonal cells with the MMTMM approach is more complex. In particular, some of the authors found numerical problems caused by matrices with high condition numbers and observed that the IBZ of hexagonal unit cells includes a boundary path that does not benefit from the linearization reported in [42]. To address the latter issue, a solution was proposed in [14], which involves solving the problematic path using a polynomial eigenvalue problem. However, this method still involves handling ill-conditioned matrix equations, which are susceptible to numerical instability and may require switching from double-to-quad precision during the MMTMM implementation. Moreover, this task often necessitates the use of specialized libraries, which complicates the implementation to preserve the method’s accuracy and efficiency. We demonstrate here that the use of symmetry planes resolves this issue, as it allows the challenging path to be addressed by a 1D-periodic eigenvalue problem. Employing a rectangular or square supercell facilitates the computation of the scattering parameters in this context. This approach is necessary because the CST-FDS with a hexahedral mesh requires that the waveguide ports be aligned with Cartesian coordinates, a requirement that cannot accommodate hexagonal layouts. Furthermore, a rectangular supercell enables us to preserve phase-shift independence in the xand y-directions, as enforced by the symmetry planes’ arrangement. The PC shown in Fig. 13(a) will be studied as an example of this type of 2D-periodic hexagonal structure. Its lattice, Fig. 13. (a) Hexagonal unit cell with a mirrored hole in each vertex. (b) Physical space, the bottom and top plate are overlapped due to the mirror symmetry. (c) 3-D model of the supercell highlighted in (b). (d) First BZ of the PC (colored blue) and the supercell (colored pink) illustrated in (a) and (c). including the hole arrangement, is illustrated in Fig. 13(b), where the PC is colored orange, and the supercell (SC) is colored light blue. A 3-D illustration of the SC used by the MMTMM to obtain the results is shown in Fig. 13(c). As elaborated in Section IV, the implementation of symmetry planes imposes a specific phase delay that facilitates the derivation of dispersion diagrams along the boundary paths of the PC/SC IBZ. However, in the current situation, this has important consequences for the calculation of the path ΓMK through the study of the SC. The SC BZ differs in size and shape with respect to the PC BZ, as illustrated in Fig. 13(d). Therefore, symmetry planes must be applied carefully to align the phase delay within the SC with the appropriate path in the PC IBZ. In particular, Table IV shows all the phase-shift relations to study the path ΓMK using the PC or SC edges. The path ΓM is studied in the SC by fixing φx=0 by means of the combination Mx 0Mx px/2. Here, we find the same situation as in the case of the square unit cell with circular mirrored hole discussed in Section IV-A, where multiple modal solutions appear due to compression of the SC BZ to half the size of the PC BZ in the y-direction (that is, a modal solution of the PC appears as two different solutions of the SC). In order to compute the odd modes, the combination Ex 0Ex px/2 is applied. As also reported in [14], the paths MK and ΓK can be calculated using the path ΓK1due to the symmetry of the unit cell. The path ΓK1can be analyzed in the SC by fixing φy=0 and imposing the combination My 0My py/2. A subsequent postprocessing discussed in Appendix Bis applied to rebuild the dispersion diagram and determine which modes of the SC path ΓK1are associated with the path MK or ΓK of the PC. Again, the odd modes are calculated using the This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. JIMENEZ-SUAREZ et al.: SYMMETRY-BASED MMTMM FOR ANALYSIS OF 2D-PERIODIC STRUCTURES 9 TABLE IV PHASE-SHIFT CONDITIONS FOR THE IBZ BOUNDARY IN FIG.13(D) Fig. 14. Dispersion diagram computed from the unit cell illustrated in Fig. 13(a). (a) Phase shift and (b) attenuation constant along the path ΓMK; dimensions (in millimeters) for the PC: a=2.9, r=0.27a,h=0.45a, and g=0.05. For SC px=a and py=√3a. combination Ey 0Ey py/2. Finally, it is found that the dispersion diagram along the PC IBZ edge is completely analyzed using only the combinations Mν 0Mν pν/2and Eν 0Eν pν/2in the SC, where νrepresents the xor y-direction. The results of the dispersion diagrams obtained following the methodology explained above are shown in Fig. 14. In this case, the HFSS ES was chosen to calculate the dispersion diagram of the hexagonal PC [13],[14]. A good agreement is observed from Fig. 14 for the results of different methods. The attenuation constant is only compared with different MMTMM approaches, as HFSS ES does not provide direct information about it. For effective convergence of the solution with the MMTMM w/o SP, 12 modes were utilized to solve the eigenvalue problem. In contrast, only two modes were required utilizing symmetry planes to calculate the scattering parameters and solve the eigenvalue problem. As expected, the use of symmetry planes significantly reduces the condition number. In Table V, a comparison of computing times using different methods is given, highlighting the benefits of the TABLE V COMPUTING TIME FOR 100 SAMPLES IN FIG.14 approach presented in this work. The use of symmetry planes in hexagonal unit cells is also found to be essential for diminishing numerical noise and computational effort, as solving the corresponding eigenvalue problem requires handling a larger structure and a greater number of modes to ensure the solution’s convergence. VI. CONCLUSION In this work, an efficient and robust approach to compute the dispersion diagram of 2D-periodic structures is presented by taking advantage of the existence of symmetry planes within their corresponding unit cells. The use of these symmetry planes allows the analysis of 2D-periodic structures to be reduced to a number of 1D-periodic cases. The analysis of 1D-periodic cases is advantageous, as the solutions are accurate, require less computational cost, and are less sensitive to numerical noise. This approach has been tested by computing the phase and attenuation constant of several examples of periodic structures with different lattices, symmetries, and motifs, demonstrating the versatility and applicability of the approach. In addition, the use of symmetry planes in the solution of the MMTMM provides a deeper understanding of the study of the BZ and the relationship between the phase shift and the path on its IBZ. APPENDIX A Typically, examining boundary paths within the IBZ of a 2D-periodic structure—periodic along xand y-axes—involves setting a phase shift in one direction and then analyzing wave propagation in the perpendicular direction. This phase-shift condition can be implemented during the postprocessing stage of the eigenvalue problem or by employing symmetry planes. In this work, these symmetry planes are treated as PEC/PMC walls, with their corresponding boundary conditions ˆ n×E=0 and ˆ n·E=0, respectively. The given boundary conditions determine a particular distribution of the electric field (E-field) for the modes. This distribution imposes a fixed phase shift condition in the direction orthogonal to the symmetry plane. Fig. 15 illustrates various examples of phase-shift conditions determined by the combination with the symmetry plane. In addition, all combinations of symmetry planes used in this work and their associated phase shift condition are tabulated in Table VI. Note that in the study of 2D-periodic structures with periodicity along xand y-directions, the required phase shift condition must also be fixed along these axes. This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination.