## Abstract

Multilayer metasurfaces have attracted much attention for broad bandwidth and high transmission efficiency. However, the electromagnetic (EM) coupling will be significant and greatly affect the transmission performance when the inter-layer distance is in subwavelength scale. Here, we present a generally equivalent circuit approach (ECA) that can be used to investigate the transmission performance and the near-field coupling mechanism in multilayer metasurfaces. Based on the proposed ECA, the analytical expression of transmission coefficient is obtained. The accuracy of the circuit model is further verified by several typical multilayer metasurfaces in ultrawideband from 0 to 30 GHz. The proposed method is valid for an arbitrary multilayer metasurface system, regardless of the layer number and the metallic elements.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Metasurfaces are two-dimensional artificially engineered materials formed by periodic/aperiodic arrangement of subwavelength unit cells. By rationally designing the geometry and arrangement of unit cells [1–3], metasurfaces can freely manipulate the electromagnetic (EM) waves like three-dimensional metamaterials. Attributing to the unprecedented ability to control EM waves, metasurfaces have received extensive applications from microwave to visible frequencies, such as EM absorption [4], polarization controls [5–7], and so on. In these applications, the transmission metasurface devices are an important branch. Due to the limitation of physical mechanism, single-layer metasurfaces have low transmission efficiencies and inadequate capabilities of phase control, thus their practical applications are hindered.

Multilayer metasurfaces possess many advantages of high transmission efficiency, wide bandwidth, and 360° phase control [8], which attract great research interest in transmission devices. Usually, ECA is an effective method to analyze and predict the EM behaviors in multilayer metasurface systems. In initial works, metasurfaces are modeled as impedance sheets, and dielectric plates are regarded as transmission lines. Based on these assumptions, multilayer metasurfaces are equivalent to a circuit network cascaded by impedance sheets and transmission lines. To characterize the transmission performance of equivalent circuit, precisely achieving the impedance sheet of each metasurface is critical. For metasurface with simple geometries, its impedance sheet can be expressed in a closed-form expression, while for complex metasurface, it is extracted usually by full-wave simulations or by curve-fitting procedures [9–11]. The existing ECA significantly simplifies the design of multilayer metasurface. When the inter-layer distance is in subwavelength scale, however, high-order harmonics will pass through the dielectric plate and affect the adjacent metasurfaces, resulting in low accuracy of ECA.

In order to address the forgoing issue, an ECA combined with Floquet’s theorem was proposed [12,13]. In this method, the periodic scattering problem is treated as the discontinuity problem in a general waveguide, where fundamental harmonic is considered as the only transmission mode, and all high-order harmonics are ignored. By investigating the electric and magnetic field distributions in this discontinuous region, one can obtain the equivalent impedance or admittance of metasurface. The achieved impedance and admittance are then characterized as the generalized impedance matrix (GIM) and admittance matrix (GAM) to further study the transmission of interested structures. The validity of this method has been verified by various single-layer structures, such as gratings [14], slot arrays [15], patch arrays [16,17], and frequency selective surfaces [18]. However, this method is only used to study single-layer metasurface and not applied to multilayer structure. Recently, Maria et al. proposed a $\pi \textrm{ - }shaped$ circuit to investigate multilayer structure with near-field coupling [19], but it is suitable to slot-based 2D periodic arrays. More recently, Olk et al. analyzed the electromagnetic coupling of multilayer metasurface by introducing ideal transformers to equivalent circuit [20], in which the circuit parameters such as the self impedance and mutual impedance are obtained from simulations, so that the physical insight into near-field coupling cannot be clearly provided.

In this paper, a novel ECA based on circuit theory and Floquet’s theorem is proposed to investigate the near-field coupling for multilayer metasurface system. In this ECA, the circuit elements are composed of self-impedance (z* _{ii}*), mutual-impedance (z

*), transmission line, and current-controlled voltage source. The self-impedance is employed to denote the surface impedance of each metasurface layer, and the mutual-impedance is used to quantify the near-field coupling. Furthermore, the self- and mutual-impedances are analytically expressed via Floquet harmonics of surface current. Consequently, the transmission and near-field coupling characteristics of multilayer system are obtained from the proposed ECA as long as we accurately calculate the surface current distributions of unit cell in each metasurface layer. The advantage of this model is that the analytical expressions of surface impedance, coupling impedance, and transimission coefficient, are all expressed by the surface current of unit cell. Therefore, the proposed ECA is not only used to study the transmission of arbitrary multilayer metasurface but also able to analyze the near-field coupling mechanism. Additionally, the ECA can also provide guidance to design multi-layer metasurface device.*

_{ij}## 2. Equivalent circuit approach

We first consider a bilayer metasurface separated by a dielectric plate, which is excited by a normal-incidence plane wave, as presented in Fig. 1. In this model, each metallic layer is presented by surface current. According to [21], the induced current density vector (${\overrightarrow J _p}$) in unit cell is denoted as

*p*denotes the

*p*th metasurface layer (here

*p*= 1, 2),

*lp*presents the

*l*th scatter in the

*p*th metasurface layer, ${c_{lp}}(\omega )$ is frequency-dependent complex amplitude coefficient, and ${\overrightarrow J _{lp}}(\rho )$ is the spatial current profile that is independent of $\omega$. The ${\overrightarrow J _p}$ can be derived from a 2-D propagation problem of a hollow metallic waveguide whose cross section is just equal to the periods of metasurface along

*x*- and

*y*-directions. For most scatters, the current profile is determined by the fundamental resonance when the wavelength of incident wave is far below the first resonant wavelength. In this case, the currents at high-order resonances are very weak, and hence they are negligible. As a result, the current profile at low frequency is considered as a static distribution of the first resonance. However, when the operation frequency is higher than the first resonant frequency, the currents at high-order resonances should be considered to improve the calculation accuracy [18]. By taking into account the current distributions of high-order resonances, the operation frequency of equivalent circuit can extend to the frequency range over the first resonance. Actually, solving the current distributions of these resonant modes is an eigenvalue problem solved by using the method of moments [22,23] or the eigenmode solver in most commercial electromagnetic simulators, such as CST Microwave Studio or High Frequency Structure Simulator (HFSS).

According to Floquet’s theorem, we decompose the bilayer system into three regions, as shown in Fig. 1(a). The EM fields in the three ranges of z < 0, 0 < z < *t*, and z > *t*, are respectively denoted by incident, scattering or transmission fields. Hence the total tangential EM fields in region 0 (z < 0) are the sum of incident and scattering fields, which are denoted as

The fields in region 1 (0 < z < *t*) are the sum of scattering fields and transmission fields, which are denoted as

The total fields in region 2 (z > *t*) only include transmission fields, which are

All the above parameters are calculated by the follow expressions.

*nm*th harmonic in medium

*i*(

*m*,

*n*= 0, ${\pm} $ 1, ${\pm} $ 2,…), ${\overrightarrow e _{q = TE}}$ and ${\overrightarrow e _{q = TM}}$ are the unit vectors of TE and TM harmonics, respectively, and $Y_{qmn}^i$ is the admittance of Floquet harmonic in medium

*i*.

The fields in regions 0 and 1 should satisfy the following boundary conditions at *z*_{1 }= 0

*a*

_{1}is the area of unit cell in 1# layer.

From Eq. (10), one further gets

Repeating the above process and using the boundary conditions at z_{2} = *t*, the following equation is obtained.

*a*

_{2}is the area of unit cell in 2# layer.

From Eqs. (11) and (12), the related variables can be written as

Then, the total electric field at the plane of *z*_{1} = 0 is denoted as

Supposing the bilayer metasurface is lossless, the current and the total electric field in unit cell (*u*.*c*.) should satisfy the integrated equation at the plane of z_{1} = 0.

By solving Eq. (15), one obtains the following equation.

The Eq. (16) can be rewritten as

where the elements are denoted asAt the plane of z_{2} = *t*, we repeat the aforementioned processes from Eqs. (1) to (15), and obtain the follow equation

*i*(

_{j}*j*= 1, 2) has the same form as that in Eq. (17), and the impedances are denoted as

It is noted that the primes in the summation in (18) and (20) exclude the fundamental TM_{00}/TE_{00} harmonic, namely *m*, *n* $\ne $ 0. By combining Eqs. (18) and (20), one obtains the impedance matrix, showing the near-field relationship between the equivalent voltages (*V*_{1} and *V*_{2}) and the equivalent currents (*i*_{1} and *i*_{2}) on each metasurface layer. This near-field relationship can be further denoted as

*(*

_{ii}*i*= 1, 2) are the self-impedance of the

*i*th metasurface layer, and the off-diagonal elements are the transfer impedances (Z

_{12}and Z

_{21}) indicating the EM coupling between 1# and 2# metasurface layer. For a practical bilayer structure, the incident plane wave is equivalent to a total voltage source (

*V*

_{1}), and a voltage drop (

*V*

_{2}) is formed when electromagnetic waves propagate in the bilayer structure. According to [20], the two-port impedance matrix in Eq. (21) can be described as an equivalent circuit presented in Fig. 1(b). The current-controlled voltage source, which is the product of mutual-impedance (z

*) and current passing through each metasurface, is introduced to quantify the near-field coupling between the two metasurface layers.*

_{ij}As for the dielectric plate, it can be regarded as a transmission line whose admittance matrix *Y _{tl}* is defined as [20]

*t*is the thickness of dielectric plate. As a result, the total admittance matrix of the bilayer system is obtained as

Then the transmission coefficient (S_{21}) of the bi-layer coupled system is achieved from the total admittance matrix

*Y*

_{0}= 1/Z

_{0}. Z

_{0}is the wave impedance in free space. This equation can accurately analyze the transmission of bilayer metasurface. Certainly, the coupling impedances (admittance) Z

_{12}(Y

_{12}) and Z

_{21}(Y

_{21}) can be used to investigate the near-field coupling property. In the following contents, for simplification, we only investigate the transmission coefficient. The investigations of coupling impedance will be mentioned in our future researches. By analyzing the coupling impedance, we will get a clear insight into near-field coupling mechanism, which can guide us to design multi-layer metasurface device with high performance.

## 3. Verification and discussion

#### 3.1 Single-layer metasurface

If we let the parameters of Z_{22}, Z_{12}, and Z_{21} be zero, the bilayer structure presented in Fig. 1(a) is degraded into a single-layer metasurface. Therefore, we first use single-layer metasurface to verify the correctness of the above theory. From Eq. (1), the surface current of single-layer metasurface is expressed as

In single-layer structure, the inter-layer coupling disappears, and therefore the current-controlled voltage sources in circuit model should be deleted. As a result, the equivalent circuit model is reduced to a parallel self-impedance, as presented in Fig. 2(a). According to the aforementioned theory, the self-impedance of single-layer metasurface is written as

Finally, the transmission coefficient (S_{21}) of single-layer metasurface is simplified as

_{0}is the wave impedance in free space.

Here we study several single-layer metasurfaces, whose unit cells are respectively rectangular patch, L-shaped resonator, I-shaped structure, and cross wire, as examples to verify the validity of the proposed method. Their geometric parameters are presented in caption of Fig. 3. In order to get the transmission from Eq. (27), accurately calculating the surface current distributions is critical. Here, the current profiles on metasurface are calculated by using commercial software CST Microwave Studio (CST-2016). As presented in Fig. 2(b), a unit cell surrounded by a rectangular box is used. Two longitudinal planes perpendicular to the z-direction are set as magnetic boundaries [see the orange areas in Fig. 2(a)], and the other four planes vertical to the *x*- and *y*- directions are defined as electric boundaries. By using eigenmode solver, several resonant modes can be calculated, among which the current of the first resonant mode is ${\vec{J}_{1}}$. Substituting the Fourier transforms of ${\vec{J}_{1}}$ into Eq. (26), we obtain the self-impedance (Z_{11}). Then the S parameter is calculated from Eq. (27). As presented in Fig. 3, both magnitude and phase of transmission coefficient calculated from the proposed ECA show perfect consistency with the simulations in a wideband from 0 to 30 GHz, demonstrating that the circuit model has a high accuracy in a broadband.

#### 3.2 Bilayer metasurfaces and inter-layer coupling

Now, we consider the bilayer metasurfaces formed by two identical patch arrays, whose unit cell is presented in Fig. 4(a). From Eq. (1), the current profiles of each metasurface are respectively written as

It is noted that ${\overrightarrow J _1}$ is equal to ${\overrightarrow J _2}$ because the two metasurfaces have identical unit cells. Like the single-layer metasurface [see Fig. 2(a)], ${\overrightarrow J _1}$ and ${\overrightarrow J _2}$ are respectively simulated by CST. Substituting the simulated ${\overrightarrow J _1}$ and ${\overrightarrow J _2}$ into Eq. (20), the self- and mutual-impedances (Z_{11} and Z_{12}) are achieved, and then the S parameters of the bilayer metasurface are calculated from Eq. (24). Figures 4(b)–4(d) present the calculated S_{21} parameters and simulations. For comparison, we also give the S_{21} parameters without considering inter-layer coupling, namely Z_{12 }= 0. It is found that the results calculated from the proposed ECA (Z_{12}$\ne $ 0) are good agreement with the simulations. While for Z_{12}=0, the calculated results [see the blue dashed lines] show significant difference with that for Z_{12}$\ne $ 0 at high frequency domain. This phenomenon clearly shows that the near-field coupling can significantly affect the transmission of the bi-layer metasurface. Moreover, the near-field coupling becomes stronger with decreasing inter-layer spacing. As we know, the inter-layer coupling is produced by high-order Floquet’s harmonics that decay exponentially in the dielectric plate. When the dielectric plate thickness becomes thin (*t* = 1 mm), more high-order harmonics can pass through the dielectric layer from one metasurface to the other, resulting in strong inter-layer coupling. In this case, an obvious difference between Z_{12 }= 0 and Z_{12} $\ne \; $ 0 appears. For a large inter-layer spacing (*t* = 4 and 8 mm), however, most harmonics decay before they reach the other metasurface, which causes very weak near-field coupling. Hence, the difference between Z_{12 }= 0 and Z_{12} $\ne \; $ 0 becomes smaller, as presented in Figs. 4(c)–4(d).

In addition, the proposed circuit model is also adaptive to the bilayer system consisting of two different metasurfaces. To demonstrate this property, we take a typical bilayer structure [see the insert in Fig. 5(a)] as an example to investigate its transmission. The unit cells of the two metasurfaces are cross wire and cutting line, respectively. Attributing to their different configurations, the surface currents at each metasurface are also diverse, which is described as

By simulating the two current distributions, respectively, the transmission coefficient of the bilayer metasurface is then calculated from Eq. (24), as presented in Fig. 5(a). The perfect consistency between the ECA results and the simulations are again observed.

If the metasurface is composed of complicated cells including multiple scatters [see the insert in Fig. 5(b)], the current profiles should be written as

#### 3.3 Tri-layer metasurfaces

Similar to the above analysis procedure, a tri-layer metasurface is further investigated. Figure 6(a) presents the transmission-reflection processes in a tri-layer metasurface. In this system, if the inter-layer spacing is small enough, the near-field coupling appears not only between two adjacent metasurface layers but also between 1# and 3# metasurface layers. Based on these considerations, the equivalent circuit presented in Fig. 6(b) is extracted. Here, the current-controlled voltage sources (*Z*_{13}*i*_{3} and *Z*_{31}*i*_{1}) are employed to describe the near-field coupling between 1# and 3# metasurface. Different from Fig. 1(b), the port 2, that is V_{2} (I_{2}), is deliberately introduced to facilitate analysis [20].

For convenience, a tri-layer metasurface consisting of three identical rectangular patch arrays [see Fig. 7(a)] is adopted as an example to demonstrate the availability of the proposed ECA. In this case, the surface currents of 1#, 2#, and 3# metasurface layesr are expressed as

Reusing (15) at the position of each metasurface, the self-impedance and mutual-impedance are obtained as following.

Writing (33) in the matrix form, we obtain

The admittance matrix of dielectric plates in this tri-layer system is

Then the near-field relationship between equivalent voltages (*V _{i}*) and equivalent currents (

*I*) on each metasurface [see Fig. 6(b)] is denoted as the following impedance matrix.

_{i}Actually, the multi-layer system is a two-port network [see Fig. 6(a)]. To reduce the 3-port network into a two-port network, one imposes the condition *I*_{2} = 0. Finally, Eq. (14) is simplified as

Figures 7(b)–7(d) present the calculated results from ECA and from simulations. As presented in Figs. 7(b) and 7(c), when we deliberately neglect all inter-layer couplings [${Z_{12}}({{Z_{21}}} )= {Z_{13}}({{Z_{31}}} )\textrm{ = }0$, see Fig. 7(b)] or neglect the near-field coupling between 1# and 3# metasurface layers [${Z_{13}} = 0,{Z_{12}} \ne 0$, see Fig. 7(c)], the results from ECA show a significant difference with simulations. While when all near-field couplings are considered [${Z_{12}} \ne 0,{Z_{13}} \ne 0$ see Fig. 7(d)], the calculated results are good agreement with simulations, demonstrating that the proposed ECA also present a high accuracy in tri-layer metasurface system. Actually, our proposed method is adaptive to multilayer system with arbitrary layer number. When the number of the layer is greater than 3, the calculation process may become more tedious.

## 4. Conclusions

We show that the multi-layer metasurface systems separated by dielectric plates exist strong near-field coupling when the inter-layer distance is in subwavelength scale. The existence of near-field coupling leads to significant inaccuracy of existing ECA, resulting in poor performance prediction of multilayer metasurface. To relieve this issue, we introduce a novel ECA that takes into account all the inter-layer coupling. In the proposed ECA, the analytical expressions of self-impedance, coupling impedance, and transmission coefficient are derived. Furthermore, these analytical expressions are expressed by surface currents of metasurface. Therefore, for an arbitrary multi-layer system, its electric performance can be effectively analyzed by calculating the surface current distributions of each metasurface layer, which greatly saves the computation time for multilayer metasurface design. For verifying the accuracy of our method, the EM behaviors of several typical metasurfaces including single-layer, bilayer and tri-layer structures are analyzed using the proposed ECA. The calculated results show a good consistence with simulations in ultra-wideband from 0 to 30 GHz.

## Funding

Natural Science Foundation of Guangxi Province (2019GXNSFDA245011, 2021GXNSFDA220003); National Natural Science Foundation of China (61761010, 62071133).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **H. T. Chen, A. J. Taylor, and N. F. Yu, “A review of metasurfaces: physics and applications,” Rep. Prog. Phys. **79**(7), 076401 (2016). [CrossRef]

**2. **H. H. Hsiao, C. H. Chu, and D. P. Tsai, “Fundamentals and Applications of Metasurfaces,” Small Methods **1**(4), 1600064 (2017). [CrossRef]

**3. **C. U. Hail, A. U. Michel, D. Poulikakos, and H. Eghildi, “Optical Metasurfaces: Evolving from Passive to Adaptive,” Adv. Opt. Mater. **7**(14), 1801786 (2019). [CrossRef]

**4. **H. Wakatsuchi, S. Kim, J. J. Rushton, and D. F. Sievenpiper, “Waveform-dependent absorbing metasurfaces,” Phys. Rev. Lett. **111**(24), 245501 (2013). [CrossRef]

**5. **X. Gao, X. Han, W. P. Cao, H. O. Li, H. F. Ma, and T. J. Cui, “Ultrawideband and High-Efficiency Linear Plarization Converter Based on Double V-Shaped Metasurface,” IEEE Trans. Antennas Propag. **63**(8), 3522–3530 (2015). [CrossRef]

**6. **X. Gao, W. L. Yang, H. F. Ma, Q. Cheng, X. H. Yu, and T. J. Cui, “A reconfigurable broadband polarization converter based on an active metasurface,” IEEE Trans. Antennas Propag. **66**(11), 6086–6095 (2018). [CrossRef]

**7. **J. B. Mueller, N. A. Rubin, R. C. Devlin, B. Groever, and F. Capasso, “Metasurface polarization optics: independent phase control of arbitrary orthogonal states of polarization,” Phys. Rev. Lett. **118**(11), 113901 (2017). [CrossRef]

**8. **A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Planar photonics with metasurfaces,” Science **339**(6125), 1232009 (2013). [CrossRef]

**9. **M. Hosseini and S. V. Hum, “A Circuit-Driven Design Methodology for a Circular Polarizer Based on Modiﬁed Jerusalem Cross Grids,” IEEE Trans. Antennas Propag. **66**(12), 7138–7147 (2018). [CrossRef]

**10. **M. Hosseini and S. V. Hum, “A Semi-Analytical Approach to Designing High Transparency Low-Profile Circular Polarizers,” IEEE Trans. Antennas Propag. **65**(10), 5322–5331 (2017). [CrossRef]

**11. **F. Costa, A. Monorchio, and G. Manara, “Efficient Analysis of Frequency-Selective Surfaces by a Simple Equivalent-Circuit Model,” IEEE Antennas Propag. Mag. **54**(4), 35–48 (2012). [CrossRef]

**12. **R. Rodríguez-Berral, F. Mesa, and F. Medina, “Analytical multimodal network approach for 2-D arrays of planar patches/apertures embedded in a layered medium,” IEEE Trans. Antennas Propag. **63**(5), 1969–1984 (2015). [CrossRef]

**13. **F. Mesa, R. Rodriguez-Berral, and F. Medina, “Unlocking Complexity Using the ECA: The Equivalent Circuit Model as An Efficient and Physically Insightful Tool for Microwave Engineering,” IEEE Microwave **19**(4), 44–65 (2018). [CrossRef]

**14. **R. Rodríguez-Berral, C. Molero, F. Medina, and F. Mesa, “Analytical wideband model for strip/slit gratings loaded with dielectric slabs,” IEEE Trans. Microwave Theory Tech. **60**(12), 3908–3918 (2012). [CrossRef]

**15. **C. S. R. Kaipa, A. B. Yakovlev, F. Medina, and F. Mesa, “Transmission through stacked 2-D periodic distributions of square conducting patches,” J. Appl. Phys. **112**(3), 033101 (2012). [CrossRef]

**16. **R. Rodríguez-Berral, F. Mesa, F. Medina, and M. García-Vigueras, “Analytical circuit model for dipole frequency-selective surfaces,” 2013 IEEE MTT-S International Microwave Symposium Digest (MTT)(2013), pp. 1–4.

**17. **S. Monni, G. Gerini, A. Neto, and A. G. Tijhuis, “Multi-mode equivalent networks for the design and analysis of frequency selective surfaces,” IEEE Trans. Antennas Propag. **55**(10), 2824–2835 (2007). [CrossRef]

**18. **F. Mesa, M. Garcia-Vigueras, F. Medina, R. Rodriguez-Berral, and J. R. Mosig, “Circuit-Model Analysis of Frequency Selective Surfaces With Scatterers of Arbitrary Geometry,” Antennas Wirel. Propag. Lett. **14**, 135–138 (2015). [CrossRef]

**19. **M. D. Astorino, F. Frezza, and N. Tedeschi, “Equivalent-circuit model for stacked slot-based 2D periodic arrays of arbitrary geometry for broadband analysis,” J. Appl. Phys. **123**(10), 103106 (2018). [CrossRef]

**20. **A. E. Olk and D. A. Powell, “Accurate metasurface synthesis incorporating near-field coupling effects,” Phys. Rev. Appl. **11**(6), 064007 (2019). [CrossRef]

**21. **F. Mesa, R. Rodríguez-Berral, M. García-Vigueras, F. Medina, and J. R. Mosig, “Simplified modal expansion to analyze frequency selective surfaces: An equivalent circuit approach,” IEEE Trans. Antennas Propag. **64**(3), 1106–1111 (2016). [CrossRef]

**22. **S. R. Rengarajan, “Choice of basis functions for accurate characterization of infinite array of microstrip reflectarray elements,” Antennas Wirel. Propag. Lett. **4**, 47–50 (2005). [CrossRef]

**23. **N. Don, A. Kirilenko, and A Poyedinchuk, “Mode basis computation for the waveguides of arbitrary cross-section,” 15th International Conference on Microwaves, Rader and Wireless Communications (IEEE Cat. No.04EX824)(2004), pp. 594–596.