Multi-modal Scattering and Propagation Through Several Close Periodic Grids

The presented method describes reflection and transmission of electromagnetic waves at multiple closely stacked metal grids using multi-modal S-parameter propagation. For that, the periodic surface admittance curve of every metal layer was determined through a novel semi-analytic approach that uses simple Fourier transformations instead of solving an integral equation. The modal components of this surface admittance were used to express the generalized scattering matrices of the individual grids. By applying multi-modal propagation techniques to the resulting scattering parameters, it was possible to model the electromagnetic interactions within a multilayer stack of periodic impedance sheets. Resulting reflection and transmission parameters perfectly matched the corresponding full-wave simulations even above the grating lobe regime. In the end, the universal mode propagation method enabled modelling of a layer stack, connected to lumped components.


I. INTRODUCTION
T HE widely studied problem of modelling periodically structured metal sheets experienced an upsurge in the past decades due to modern investigation of metamaterials and frequency selective surfaces. This is thanks to the promising application fields these impedance layers open up, starting from simple filters [1] over polarisers [2], perfect absorbers [3], reflect arrays [4], [5], thin lenses [6], and even active structures [5], [7], to only name a few. Their capability of manipulating electromagnetic waves applies to bands from microwave up to optical frequencies and is therewith of interest in various engineering fields as well as in natural sciences [8]- [11]. Nowadays, commercial electromagnetic simulators are used to accurately extract the scattering properties of one-or twodimensionally periodic layers. These tools provide universal structure investigation possibilities via finite elements or finite difference methods and often serve as reference for analytical modelling approaches. However, depending on the structure Manuscript received August 17th, 2021; revised January 17th, 2022; accepted January 22nd, 2022. Date of publication mm dd, yyyy; date of current version January 22nd, 2022. This work has been jointly supported by Infineon Technologies Linz and by Silicon Austria Labs (SAL), owned by the Republic of Austria, the Styrian Business Promotion Agency (SFG), the federal state of Carinthia, the Upper Austrian Research (UAR), and the Austrian Association for the Electric and Electronics Industry (FEEI).
Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. complexity, a full-wave solver requires significant computation power to thoroughly calculate the desired magnitudes and it does not provide much insight into the physical nature of the problem. In contrast to this, a quasi-analytic solution often applies integral equations in combination with spectral method of moments techniques [12]- [14]. The resulting generalized scattering, impedance or admittance matrices describe then electromagnetic wave scattering at a grid [15]- [18]. Interestingly, a similar behaviour was investigated years before at waveguide discontinuities [19]- [21]. The equivalent multimodal networks describing these discontinuities were applied to layer systems incorporating dielectric slabs and later lumped elements [16], [17], [22]- [24]. Full-analytic modelling of periodic metal sheets typically results in scalar, frequency dependant shunt admittances which are calculated based on Bloch-Floquet theory [25]- [27]. These equivalent admittances can be integrated into rather simple equivalent circuits, which give frequency responses accurate enough for many applications [28], [29]. However, if it comes to multiple stacked grids, there occur evanescent (nonpropagating) waves at and through the sheets. These interactions certainly influence the system's frequency response, but cannot be modelled through scalar admittances. To fully consider interactions between two or more close periodic metal sheets, one needs to solve the single layer integral equations and perform multi-modal propagation through generalized scattering matrices [18], [30], [31]. Thereby, the equivalent multi-port networks, corresponding to the metal sheets, interconnect all considered propagating and evanescent modes. A combined scattering matrix of the whole layer system, incorporates then multi-modal scattering between the various layers and contains information about reflected and transmitted propagating modes. It is further possible to define more complex equivalent networks that directly model multimodal propagation and therewith evanescent interactions between grids [32], [33]. However, these networks require careful consideration of grid alignment and layer stack, which may complicate structure investigation. In contrast to the previous works, the here presented method uses existing analytic solutions of the electric field or surface current from e.g. [14], [26] to calculate the admittance curve within a grid period. Thereby, solely Fourier transforms and characteristic waveguide admittances at both sides of the metal layer need to be applied to the analytic expressions. The generalized scattering matrices are subsequently extracted from the modal components of the resulting surface admittance [34]. It needs to be mentioned that the utilization of analytic  field or current density solutions as well as the extraction of the generalized scattering matrices only apply to 1-D periodic grids. An extension to arbitrary 2-D periodic structures, that are treated in e.g. [18], [31], would still require the Method of Moments technique. However, for structures like slits and strip arrays, the presented method spares the solution of integral equations, on the one hand. On the other hand, it does not require partly complex equivalent networks, as multimodal propagation is described through generalized scattering parameters from modal surface admittances. Fig. 1 depicts an arbitrary 1-D periodic grid on a dielectric substrate with the periodic part of the electric field E P and surface current density J P depending on the incident wave polarization. A fundamental approximation of the proposed method, which in the end enables the application of predefined expressions for field and current, is the frequencyindependence of the shape of these quantities. According to [22], [24], [27], this is valid until high order harmonics of the current in a scatterer or the field in an aperture are excited by the incident wave. Metasurfaces and frequency selective surfaces are typically designed with periodicities D smaller than the wavelength, thus the approximation above is applicable in most cases. The frequency dependence of the presented scattering parameters comes merely from the considered characteristic waveguide admittances.
November 22nd, 2021 II. THEORY In equivalent circuit theory, scattering of electromagnetic waves at periodic metal strips is usually depicted as reflection and transmission at shunt impedances or admittances. Single surface admittances are calculated by investigating induced currents and reflected or transmitted fields at the metal surface. According to Floquet theory, the linear differential equations which describe electromagnetic quantities at x-periodic structures, result in surface tangential fields E S (x) or currents J S (x) that can be expressed as periodic functions modulated by a complex number q (see e.g. [25]) as In fact, q corresponds to the incident wave vector component parallel to the surface k ∥,i = |k ∥,i | = |k i | sin(θ), which is shown in Fig. 2. The periodicity of the field E P (x) = E P (x + nD), with n ∈ Z and D as size of a grid unit cell, allows to expand it into its modal components with well defined wave vectors k ∥,n = n 2π D : Obviously, this expansion can also be applied to the surface magnetic field H S (x) with the induced sheet current density as the field difference between both sides of the metal layer In the following, these principles are used to generate expressions for modal surface admittances of periodic metal sheets.

A. Surface Admittance Extraction
Generally, the total electric field E(x) in an aperture consists of the parallel component of the incident (external) field E i and a local field generated by the surface E S (x). The same accounts for the total magnetic fields H I (x) and H II (x) at both sides of the surface. The ratio between the parallel components of current density and electric field is then defined as surface admittance This leads to a periodic function for Y S (x), which can again be expanded into modal components. One could think, it is essential to know the electric field E(x) and current density J(x) to determine the surface admittance. However, the modal components of the magnetic field are directly related to the electric field components through well known characteristic waveguide admittances for TE polarization kiεr η0βn for TM polarization.
Here, η 0 = µ 0 /ε 0 is the vacuum wave impedance, k i = ω 0 √ ε 0 µ 0 is the magnitude of the incident vacuum wave vector, and is the imaginary part of the propagation constant (phase constant). Obviously, the permittivity ε r and β n may differ at the two sides of the surface. With the waveguide admittances above it is possible to derive an expression for the H-field from the E-field or vice versa and calculate Y S (x) and its modal components. Using (2) and k n = k ∥,i + k ∥,n , the expansion of the total field reads A similar expression can be found in [26], where the field components are normalized by E i . As E 0 is the amplitude of the excited ground mode, the ratio between E 0 and E i equals the reflection coefficient In contrast to the tangential electric field, which only exists in the aperture of the grid and therewith is continuous through it, the tangential magnetic field differs at both layer sides due to excited current densities. Furthermore, H I (x) and H II (x) depend on the dielectric surrounding at the respective side of the surface. Expressing the modal components of the magnetic field through the electric field components according to (4) gives With all electric field components E n and E i known, the magnetic fields were determined and the surface admittance could be calculated. However, the zero order mode component of the Fourier transform of E(x) equals the sum E i + E 0 . To separate E 0 from E i + E 0 one needs to know R. According to propagation theory, the reflection coefficient of a shunt admittance between two different characteristic admittances is Here, Y eq is the sheet equivalent admittance from [26], which is defined by a multi-modal equivalent circuit through parallel input admittances It is mentionable that in this work Y eq is calculated merely from the electric field's Floquet components and the characteristic waveguide admittances. This is achieved, as the influences of dielectrics with finite thickness or neighbouring grid structures are treated separately in the multi-modal propagation calculations. Now that R is known, one can express the surface current density through (8), (9), and (7) as Once the modal admittances and the electric or magnetic field are known, the surface admittance is calculable.
In principle the procedure of calculating J(x) can be broken down to a discrete Fourier transform (DFT) of the electric field, a scalar multiplication of the coefficients E n with the characteristic admittances Y n , and an inverse DFT of the current density components J n , as illustrated in Fig. 3. This can be simplified by applying the convolution theorem. The multiplication in Fourier domain becomes then a convolution in real space is the inverse DFT of the characteristic admittance series. Rewriting (12) with the identity above makes Obviously, both expressions for J(x) give the same results. In some cases, it is favourable to calculate the electric field in the aperture from the surface current density of the strip. This becomes important at grids excited by TM polarized waves, as here the curve of the surface current density appears continuous, whereas the electric field in the aperture becomes discontinuous at the grid edges (see e.g. Fig. 6b). The current density is then expanded by Certainly, it is possible to perform the same extraction as above the other way round, meaning the functional direction depicted in Fig. 3 turns. The major difference comes then from the calculation of the sheet equivalent admittance, as the corresponding equivalent circuit changes. The shunt admittance Y eq is now represented by a series connection of the modal input admittances Using this in (10), the reflection coefficient of the single layer is determined and the electric field in the aperture reads or with the convolution theorem applied Both approaches above provide straight forward surface admittance calculation, independent on whether the electric field or the surface current density is known. Furthermore, one can see from (6) and (12) or from (15) and (17) that the amplitude of the a priori known field or current density does not influence Y S (x). A scalar factor on E(x) or J(x) would apply to every modal component and cancel in (11) or (16) and later in (3). Therefore, only the shape of the aperture field or the surface current needs to be known. Related analytical expressions are found in e.g. [14], whereas the field profile can also be extracted numerically from full-wave simulations. The normalized aperture field and current density are approximated as frequency independent. This is valid for slits and strips that are small compared to the vacuum wavelength [26].

B. Multi-modal Propagation
To investigate scattering of an incident wave at serried impedance sheets into high-order modes, one needs to expand the single surface admittances into its modal components The modal admittances g n directly relate the tangential field components to the surface current components as in [34] through and indicate how E(x) and J(x) are linked in a corresponding equivalent circuit (see e.g. [16], [20]). Apparently, the periodic sheet interconnects discrete modes and distributes power depending on the properties of the impedance layer, where the modal admittances can be seen as coupling coefficients. In Fig. 4, a network resulting from the distributive nature of two periodic grids combined by a dielectric is shown. As high-order modes are increasingly attenuated, the amount of practically needed g n as well as the sum in (20) is finite. That allows to practically express the surface current components J n through a multiplication of an admittance matrix Y S (Toeplitz matrix), containing all modal admittances, with a column vector consisting of modal where N is the highest considered mode number [34]. Subsequently, the column vectors in (21) will be written as In fact, there are now 2N +1 independent modes transmitted and reflected at the metal grid, which makes every layer a 2 × (2N + 1)-multi-port. The propagation of an incident wave through several of such multi-ports above can be described by multi-modal transmission parameters where A, B, C, and D are (2N + 1) × (2N + 1) matrices. Consequently, an impedance sheet is modelled by transmission parameters of the shunt surface admittance with where I is the identity matrix of rank 2N + 1. Intermediate dielectric slabs of thickness d are represented by a lossless multi-modal transmission line where Sin and Cos are diagonal matrices containing sin(β n d) and cos(β n d), respectively. The dielectric admittance matrix Y d is also diagonal with its non-zero elements being the material specific characteristic waveguide admittances from (4). In all sub-matrices of (24), n indicates the row and column index which goes from −N to N .
Several circuit elements are finally combined by consecutive multiplication of the single transfer matrices where T tot represents the whole layer stack. Reflection and transmission properties for TE and TM polarization can now be extracted from T tot , by relating incident, reflected and transmitted modes as in [34]. However, depending on the desired model accuracy and the dimensions of the structure, the calculations above require large mode numbers. This can become an issue as the maximum value entries in T d increase exponentially with N , leading to floating point precision errors in the final transmission matrix.

C. Numerical Stability Improvement
For a certain high mode number m ≫ D λ0 √ ε r , the imaginary part of the propagation constant from (5) is approximated by With this, entries of the matrices Sin and Cos in (24) become and Both hyperbolic expressions above exceed the precision limit of a double floating-point number (2 53 ≈ 9 · 10 16 ) already at These high numbers lead to a singular admittance matrix [19], [21] and subsequently to instabilities at the conversion of transfer into scattering parameters, which finally limits the number of calculable modes at certain slab thickness d. Furthermore, a stack of dielectric layers, implying a multiplication of matrices with high value entries, additionally reduces the maximum mode number due to catastrophic cancellation (see e.g. [35]). In the following, the stability problem emerging from consecutive multiplication of transfer parameters and the one coming from single thick dielectrics will be solved.
To avoid the increase of high numbers through matrix multiplication in (25), single equivalent circuit elements have to be combined as scattering parameters [19]. These are defined by where ⃗ E ± = 1 2 ( ⃗ E ± Z 0 ⃗ H) are the incident and reflected waves at the respective side (I or II) of the multi-port. The diagonal impedance matrix Z 0 = Y 0 −1 is the inverse of the vacuum admittance matrix which consists of the free-space characteristic waveguide admittances. By transforming (22) into Z-parameters one can express the incident and reflected waves from (28), as in [36], through Z and ⃗ H: The scattering parameters result from eliminating ⃗ H in (30), which finally gives (31) Multiple elements of the equivalent circuit are now combined stepwise by evaluating the total S-parameters S (1+2) of two individual consecutive scattering matrices S (1) and S (2) (see e.g. [19], [37]): The formulas above are applied until the resulting Sparameters represent the whole layer stack. Reflection and transmission parameters correspond then to S 11 and S 21 , with S 11 = S 11 [0, 0] and S 21 = S 21 [0, 0] as the reflection and transmission coefficient of the propagating ground mode. In some cases, already a single dielectric layer in combination with high required mode number leads to floating point numbers exceeding the double precision limit. To circumvent this problem, the dielectric can be artificially separated into p equally thin layers such that the corresponding S-parameters can be calculated and combined through (32) afterwards. The number of separations needed per layer is according to (27)  In Fig. 5, the principle of converting transmission parameters into the total scattering matrix is depicted. Here, the evaluation of a double layer with e.g. D = 3.0 mm, d = 0.5 mm, and m = N = 40, requires the dielectric slab to be split up into two consecutive sheets, exhibiting half of the actual thickness. Once converted into scattering parameters, all elements are combined to a total matrix. Obviously, the two calculation steps above increase the amount of performed matrix multiplications and therewith the computation time. However, in most applications presented in this work, the reached calculation accuracy without stability improvement was not sufficient.

III. FULL-WAVE SIMULATION
Subsequently, resulting fields and frequency responses from the above presented approach will be compared to full-wave simulations. All numerical operations were performed on a notebook with an Intel(R) Core(TM) i7-7820HQ CPU and 32 GB RAM, where the simulation software Ansys HFSS(TM) from Ansys Electronics 2020 R2 occupied three cores of the processing unit. The grids were modelled as infinite periodic arrays with square unit cells through Master/Slave boundaries and Floquet Ports. Metal structures within the unit cells consisted of two-dimensional perfect conductors and the dielectrics were considered to be lossless.
Even if the subsequently presented scattering parameters are shown over a broad frequency range, at most structures a single mesh with corresponding frequency sweep resulted in decent frequency responses. Thereby, meshing was performed at a quarter of the maximum frequency f SIM = f max /4. Only in few cases (Fig. 10), computationally expensive broadband meshing was needed to model the multiple resonances of more complicated structures. The mode calculator in HFSS evaluates the attenuation coefficients for the different Floquet modes according to f SIM . With this, the number of simulated modes M was chosen in order to only consider modes with an attenuation below 20 dB/m. The final frequency responses were evaluated by an interpolating frequency sweep with 201 linearly distributed sample points.

IV. INVESTIGATION OF SINGLE LAYER GRIDS
Before the scattering properties of multi-layer systems are shown, the calculated tangential fields and the surface admittance of single layers are presented. This gives an important insight into the relation between electromagnetic properties under TE and TM polarized excitation. Furthermore, not only the frequency response, but also the extracted fields can be compared to full-wave simulations. In the following, analytic expressions proportional to the aperture field or patch current density from [14] or [26] are applied to the procedure above. Depending on the polarization of the incoming wave, the corresponding formulas read where the expressions for E(x) apply to the grid aperture and those for J(x) are valid at the metal sections. For both polarization states, the field or current with the formula exhibiting a positive exponent (E TE and J TM ) is centred and preferably used in the subsequent calculations. This is important to facilitate and stabilize the later applied Fourier transformation as E TM (x) and J TE (x) are discontinuous at x = ± D−w 2 . Obviously, surface admittance extraction depends on the number of sampling points within a structure period, where an odd number of points is favourable to resolve the maximum in the center. In the following, the fields of the grid unit cells are discretised by 601 equally distributed samples. Fig. 6 depicts calculated and simulated quantities of a grid under TE or TM polarization. The electric field used for computations related to TE excitation corresponds to (34a) in the aperture and is considered to be nought at the metal parts (thick yellow lines in the plots 6a to 6d). To prevent a division by zero in (3), a comparably small number (e.g. 10 −10 ) is added to E(x), which is corrected for in the subsequent calculations. In Fig. 6a, E(x) is additionally scaled to match the field resulting in the full-wave simulation. This affected the resulting current density J(x), but not the surface admittance, shown in Fig. 6c. As expected, Y eq is in the range of Y S (x) averaged over the grid unit cell, which further corresponds to g 0 . The frequency responses are extracted from Y S (x) by applying the transformation from modal admittances to scattering parameters. In Fig. 6e, a considered mode number of N TE = 40 or above results in a frequency response of a single grid that matches the HFSS simulation very well, even above the grating lobe region where n ≥ 1. As anticipated, the grid behaves as shunt inductor under TE polarized excitation, resulting in a negative imaginary Y S (x) and consequently high-pass behaviour. The corresponding scalar impedance from e.g. [25] results in scattering parameters that fit the multi-modal approach and simulations at low frequencies.
However, close to the grating love regime, the single-mode solution becomes inaccurate even at single layer grids without evanescent mode propagation. Under TM polarized excitation, J(x) is used to extract electromagnetic quantities. In Fig. 6b, the current density corresponds to expression (34d) at the metal parts, whereas the aperture is considered current free. As above, an offset is added to J(x) to gain numerical stability, which is corrected for when calculating E(x). The extracted current density from the HFSS simulation matches the scaled analytic formula well. Fig. 6d depicts once more the rough conformity between the average of Y S (x) and Y eq . As expected, the grid reacts capacitively under TM excitation with positive imaginary surface admittance and low-pass behaviour. It is noticeable that proper frequency response extraction of grids under TM polarized excitation needed more modes than under TE polarization. Therefore, the frequency response in plot 6f was determined with a maximum mode number of N TM = 60. The mode numbers above are applied to all the subsequent frequency responses. Again, the reflection magnitude of a scalar, singlemode surface impedance fits the simulations well at low frequencies.
Obviously, the extraction of scattering parameters from Y S (x) of single grids does not require multi-modal propagation. Nevertheless, higher order modal admittances g n are needed to accurately compute frequency responses (see e.g. Fig. 8, 12, and 13 in [31]), as they influence the fundamental mode reflection and transmission. The here presented scattering parameter S 11 incorporates higher order modal admittances g n at the matrix operations (29) to (31).

V. MULTI-LAYER STRUCTURES
The actual investigation of evanescent mode propagation was performed on multiple serried layers. In the following, the extracted surface admittances of different grids as e.g. in Fig. 6 and corresponding modal admittances are used to derive the scattering at multi-layer systems. To see, whether the interaction between grids is modelled properly, two double-layers are investigated, which only differed by the lateral alignment of their metal parts. The bottom grid of the double layer in Fig. 7a is aligned with the top one, whereas it is shifted by half of the periodicity D in Fig. 7b. This displacement is modelled by shifting the calculated admittance of the bottom grid according to Y S,2 (x) = Y S,1 (x+ D/2), which affects the modal admittances and further the scattering parameters. By investigating both systems through  Fig. 7 under TM polarization. Plot (a) corresponds to the structure in Fig. 7a and (b) corresponds to Fig. 7b, respectively. a single-mode equivalent circuit, their frequency responses are equal, since the related shunt impedances of the layers are the same. However, one can imagine that a variation of the distances between consecutive strips changes existing and generates additional resonances. Exactly this can be seen when comparing Fig. 7c with 7d, where both multi-modal frequency responses match the corresponding simulations.  As a prove of concept, the two structures of Fig. 7 were also investigated under TM polarization. Fig. 8 depicts a significant difference in the frequency responses of the structures between 60 GHz and 90 GHz. The multi-modal calculations again match the related simulations over the full frequency range.
To understand the influence of certain high order modes on the resulting scattering parameters, one has to examine the transmission parameters in (24). These depend on the distance between the metal grids d and the mode propagation constant β n , which is further related by (5) to incidence angle, frequency, structure periodicity, and dielectric permittivity. In Table I, examples for propagation magnitudes at four different frequencies are presented. Here, the propagation constants of the first high order modes are either imaginary or real, resulting in an exponential decay or constant propagation of One can see that evanescent modes propagate between close layers even at low frequencies related to the grid periodicity D. However, with increasing grid distance d, the influence of these modes vanishes and a single mode approach could accurately model multilayer scattering. As a further prove of the proposed model's correctness, the structure in Fig. 7a was investigated under oblique incidence. Fig. 9 depicts frequency responses with θ = 30 • and θ = 60 • . One sees that the tilt of the beam under TE polarization results in an additional resonance around 62 GHz and simultaneously flattens the second minimum at 90 GHz, similar as in plot 7d.
For an oblique incident wave, the grid seems narrower with the second layer appearing laterally shifted. Under TM polarization, the two minima between 60 GHz and 80 GHz shift to lower frequencies and the reflection magnitude decreases with increasing incidence angle. It is mentionable that the reflection coefficient rises when tilting the beam further than θ = 60 • , which is an indication to be close to Brewster's angle. Again, the calculated frequency responses match the simulations. Under the presented incidence conditions, the expressions for the scalar impedances in (36) and (37) are independent of θ. Therefore, the frequency response calculated by the single-mode solution becomes increasingly inaccurate with steeper incidence.
Finally, the frequency response of the multi-layer structure shown in Fig. 10a is examined. This structure exhibits four interacting grids with strongly varying widths and dielectric slabs of different thickness and permittivity. As shown in Fig. 10b and 10c, the scattering behaviour induced by various mismatches is extracted properly for TE and TM wave incidence. When simulating the frequency responses of this multilayer with a single mesh operation, mismatches occurred at higher frequencies (f > 150 GHz). Therefore, broadband meshing was applied here, which, however, increased the former simulation time of around half an hour tremendously to more than 20 hours. It is mentionable that the defined mode numbers N TE and N TM did not have to be adjusted for any of the presented examples. This means, once the necessary number of required modes is reached, the surface admittance Y S (x) is properly represented by its generalized admittance matrix Y S . A further increase of N TE and N TM barely affected the coincidence between simulation and calculation, but increases the computation time.

VI. SURFACES WITH LUMPED ELEMENTS
Besides quick and accurate analysis of periodic multilayers, the presented method is also capable of investigating the effects of lumped circuit elements on such structures. This becomes important when designing active metasurfaces, where components on top or in between periodic grids can influence the system's scattering behaviour. Subsequently, the frequency selective surface shown in Fig. 11a with a lumped impedance Z c at the bottom layer will be investigated under TM polarisation. The impedance interconnects two neighbouring grid lines and is represented by a serial connection of a resistance R c , an inductance L c , and a capacitance C c . Depending on the considered mode number n, the effective impedance reads Interestingly, only the real part of the phase constant contributes to the reactive properties of the lumped component. The related imaginary part, which can be seen as damping factor in the time domain, must be neglected at evanescent modes. According to the multi-port system shown in Fig. 11b, the lumped component can then be described through the Zparameters where Z c is a diagonal matrix consisting of Z c [n, n] = Z c,n . The corresponding scattering parameters, following from (31), are then added to the combined scattering matrix of the whole layer stack. Obviously, this procedure is very efficient when it comes to investigations on a certain layer system loaded with a varying impedance. In that case, only the S-parameters of the lumped component need to be calculated repeatedly, whereas the multilayer parameters remain untouched.
To verify the correctness of the approach above, the presented structure was again modelled and evaluated in Ansys HFSS. A unit cell of the corresponding structure is shown in Fig.  11c, where a red strip between two grid lines represents the lumped component. As the strip exhibits an intrinsic, structural inductivity of L strip = 0.52 nH, the set inductance needs to be lowered accordingly. In contrast to the performed calculations, a complete full-wave simulation needs to be carried out after every change at the impedance boundary. Fig. 11d and 11e depict the good match between HFSS simulation and calculation. Here, the lumped impedance adds an additional resonance at lower frequencies (≈ 10 to 25 GHz), whereas the frequency response above 60 GHz remains almost unchanged, even if Z c is removed completely. The slight mismatch at the two frequency responses from 75 to 100 GHz might come from the additional scattering at the strip structure that is not considered in the quasi analytic calculations.

VII. COMPUTATION EFFICIENCY
In general, the full-wave simulations applying Master/Slave boundaries and Floquet ports provide a reliable reference for the electromagnetic quantities of the investigated structures. However, finding the optimal simulation settings, which grant time efficiency on the one hand and high accuracy over the whole frequency range on the other hand is not always straight forward. Ideally, the mode calculator as well as the meshing procedure are run at the highest occurring frequency, which indeed would result in a huge number of mesh elements and therewith would exacerbate mesh convergence. To be able to investigate the presented structures efficiently, without the use of a powerful simulation server, f SIM was lowered according to section III. The chosen simulation settings are a good compromise between computation effort and accuracy.
To depict the effectiveness of the proposed method, the approximate time taken for calculating or simulating the presented frequency responses are shown in Table II. As the simulated S-parameters, also the calculated responses consist of 201 sample points. One can see that the quasi-analytic solution always finishes within one minute even for the most complicated layer stacks, whereas the simulation time rises dramatically with structural complexity due to sophisticated meshing.

VIII. CONCLUSION
Summarizing, this paper presents an accurate and efficient approach to calculate multi-modal scattering parameters of several stacked grids by using predefined solutions of the electric aperture fields or the current density on the grids. Although these solutions only hold under certain circumstances as e.g. low frequencies or narrow grids/slits, the procedure works valid over a broad frequency range and for most grid widths. Additionally, its modularity allows a joint consideration of passive, wave scattering structures and lumped impedances of e.g. circuit elements. The quickness of the method is especially convenient when it comes to broadband optimization problems, as e.g. in [6], or when the continuous change of the frequency response with a design parameter, as e.g. a lumped impedance, is of interest. Both of these challenges are found in the field of active metamaterial and frequency selective surface design. However, since the surface admittances of grids under TE polarization (inductive) cannot be modally combined with grids under TM polarization (capacitive), it is not possible to design high order metasurface bandpasses efficiently. This limitation could be overcome by applying the presented method to grids exhibiting 2-D periodicity. To do that, surface and current distribution of such surface structures as well as the corresponding equivalent networks need to be investigated properly. Nevertheless, various combinations of grids and dielectric slabs can provide sufficiently broad pass-or stopbands with corresponding transmission or reflection phase shift.