Accurate metasurface synthesis incorporating near-field coupling effects

TL;DR

This paper introduces an extended equivalent-circuit model for metasurface synthesis that accurately accounts for near-field coupling effects, improving design precision especially at higher frequencies.

Contribution

The authors develop a novel extension of the impedance sheet model to include near-field coupling, enabling more accurate metasurface synthesis.

Findings

Extended model predicts scattering parameters accurately

Synthesized metasurface achieves near-theoretical efficiency

Improved design algorithm compensates for near-field perturbations

Abstract

One of the most promising metasurface architectures for the microwave and terahertz frequency ranges consists of three patterned metallic layers separated by dielectrics. Such metasurfaces are well suited to planar fabrication techniques and their synthesis is facilitated by modelling them as impedance sheets separated by transmission lines. We show that this model can be significantly inaccurate in some cases, due to near-field coupling between metallic layers. This problem is particularly severe for higher frequency designs, where fabrication tolerances prevent the patterns from being highly-subwavelength in size. Since the near-field coupling is difficult to describe analytically, correcting for it in a design typically requires numerical optimization. We propose an extension of the widely used equivalent-circuit model to incorporate near-field coupling and show that the extended model can predict the scattering parameters of a metasurface accurately. Based on our extended model, we introduce an improved metasurface synthesis algorithm that gives physical insight to the problem and efficiently compensates for the perturbations induced by near-field coupling. Using the proposed algorithm, a Huygens metasurface for beam refraction is synthesized showing a performance close to the theoretical efficiency limit despite the presence of strong near-field coupling.

Tables

Table 1. Table 1: Geometric parameters of the cells in the refracting metasurface.

Geometry Parameters [$\mu$m]
cell	1				2				3
parameter	$S_1$	$G_1$	$S_2$	$G_2$	$S_1$	$G_1$	$S_2$	$G_2$	$S_1$	$G_1$	$S_2$	$G_2$
before	300	126	400	149	300	204	172	220	300	175	400	138
after	300	114	400	139	300	210	172	218	300	172	400	137

Equations

$$\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right]=\left[\begin{array}{ccc}{Z}_1& Z_{mut}^{12}& Z_{mut}^{13}\\ Z_{mut}^{21}& Z_2& Z_{mut}^{23}\\ Z_{mut}^{31}& Z_{mut}^{32}& Z_3\end{array}\right]\left[\begin{array}{c}{i}_{z,1}\\ i_{z,2}\\ i_{z,3}\end{array}\right].$$

$$\left[\begin{array}{c}{i}_{in,1}\\ i_{in,2}\\ i_{in,3}\end{array}\right]={\mathbf{Y}}_{\mathbf{tot}}\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right].$$

$$\left[\begin{array}{c}{i}_{tl,1}\\ i_{tl,2}\\ i_{tl,3}\end{array}\right]={\mathbf{Y}}_{\mathbf{ff}}\left[\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right],$$

$${\mathbf{Y}}_{\mathbf{tot}}={\mathbf{Y}}_{\mathbf{nf}}+{\mathbf{Y}}_{\mathbf{ff}}.$$

$$\left[\begin{array}{c}{i}_{in,1}\\ i_{in,3}\end{array}\right]={\mathbf{Y}}_{\mathbf{tot}}^{\prime }\left[\begin{array}{c}{v}_1\\ v_3\end{array}\right],$$

$$\text{min}(|S_{21}^{\mathrm{re}}-S_{21}^{\mathrm{des}}{|}^2).$$

$${\mathbf{Z}}_{\mathbf{nf},\mathbf{p}}={\mathbf{Z}}_{\mathbf{nf},\mathbf{p}-1}+\gamma \text{diag}(Z_{1,p-1},Z_{2,p-1},Z_{3,p-1})$$

$${\mathbf{A}}_{\mathbf{Zn}}=\left[\begin{array}{cc}1& 0\\ 1/Z_n& 1\end{array}\right]$$

$${\mathbf{A}}_{\mathbf{tl}}=\left[\begin{array}{cc}\cos (\beta t)& j\eta \sin (\beta t)\\ j1/\eta \sin (\beta t)& \cos (\beta t)\end{array}\right].$$

$${\mathbf{Y}}_{\mathbf{ff}}=\left[\begin{array}{ccc}\frac{-i}{\eta \tan \beta t}& \frac{i}{\eta \sin \beta t}& 0\\ \frac{i}{\eta \sin \beta t}& \frac{-2i}{\eta \tan \beta t}& \frac{i}{\eta \sin \beta t}\\ 0& \frac{i}{\eta \sin \beta t}& \frac{-i}{\eta \tan \beta t}\end{array}\right].$$

$${\mathbf{Y}}_{\mathbf{tot}}^{\prime }=\left(\begin{array}{cc}{Y}_{\mathrm{tot},11}-\frac{Y_{\mathrm{tot},21}}{Y_{\mathrm{tot},22}}{Y}_{\mathrm{tot},12}& Y_{\mathrm{tot},13}-\frac{Y_{\mathrm{tot},23}}{Y_{\mathrm{tot},22}}{Y}_{\mathrm{tot},12}\\ Y_{\mathrm{tot},31}-\frac{Y_{\mathrm{tot},21}}{Y_{\mathrm{tot},22}}{Y}_{\mathrm{tot},32}& Y_{\mathrm{tot},33}-\frac{Y_{\mathrm{tot},23}}{Y_{\mathrm{tot},22}}{Y}_{\mathrm{tot},32}\end{array}\right).$$

$${\mathbf{Y}}_{\mathbf{ff}}^2=\left(\begin{array}{cc}\frac{-j}{\eta \tan \beta t}& \frac{j}{\eta \sin \beta t}\\ \frac{j}{\eta \sin \beta t}& \frac{-j}{\eta \tan \beta t}\end{array}\right).$$

Full text

I Introduction

Metasurfaces are thin sheets of subwavelength resonators which have emerged as a versatile platform for wavefront manipulation, addressing applications from microwave to visible frequencies Yu and Capasso (2014); Ding et al. (2018); Chen et al. (2018). Where high transmission efficiency is required, Huygens metasurfaces are widely used as they feature equal electric and magnetic dipole responses to suppress spurious reflections Pfeiffer and Grbic (2013a); Monticone et al. (2013); Decker et al. (2015). Different metasurface architectures and corresponding synthesis methods have been proposed within the last two decades Ding et al. (2018). For microwave to terahertz frequencies, most reported work is based on planar stacks of three patterned metallic layers, separated by dielectric substrates, as this architecture is compatible with planar circuit manufacturing techniques Chen et al. (2018); Epstein and Eleftheriades (2016a); Pfeiffer and Grbic (2014). In addition to Huygens metasurfaces, the same approach has been used to realize bianisotropic metasurfaces, which are more efficient for very large refraction angles Wong et al. (2016); Cole et al. (2018); Fathnan and Powell (2018).

In synthesizing metasurfaces, it is often necessary to assume that meta-atoms are small compared to the wavelength. This has motivated the development of deeply subwavelength meta-atoms Wong et al. (2014); Epstein and Eleftheriades (2016a), with feature sizes on the order of $\lambda/200$. For metasurfaces operating at W-band frequencies (75 - 110 GHz) and higher, the resulting feature sizes are too small for conventional printed circuit board fabrication techniques, so simpler geometries such as dogbone Lavigne et al. (2018); Capolino et al. (2013); Rabinovich and Epstein (2018) must be used. This frequency range is of interest for applications such as wireless communication Rappaport et al. and radar Patole et al. ; Hasch et al. , and metasurfaces should enable devices such as efficient transmitarray antennas Pfeiffer and Grbic (2013b); Jiang et al. . Despite the high potential for applications, to the best of our knowledge, there is only one work reporting transmissive refracting metasurfaces operating at W-band frequencies Pfeiffer and Grbic (2013b). Efficient modelling of near-field interaction would enable metasurfaces to be synthesized with simpler geometries, which are feasible to fabricate for this frequency range.

In modelling multilayered printed circuit metasurfaces, each metallic layer is represented by an equivalent surface impedance and the dielectric layers in between are treated as transmission line sections (See Fig. 1) Epstein and Eleftheriades (2016b); Pfeiffer and Grbic (2013b); Monticone et al. (2013). The advantage of this model is that it yields expressions for the required surface impedance of each metallic layer, greatly simplifying the process of designing the metasurface. The transmission line model accounts for the propagation of the fundamental Floquet harmonic within a periodic metasurface, neglecting all higher-order harmonics, since they are evanescent. Neglecting evanescent waves is accurate for structures with large separation relative to the wavelength, however in metasurfaces it is usually required that the separation between the metallic layers is small. Thus near-field coupling between the layers is neglected. As we show here, this approximation is significantly inaccurate in cases of practical interest, and can lead to performance degradation. Therefore, accurate modelling and synthesis of transmissive metasurfaces must take account of near-field effects between metallic layers. We note that the near-field coupling effects that are investigated here are not related to the evanescent waves (or auxiliary fields) deliberately introduced in Ref. Epstein and Eleftheriades, .

Early works on near-field interaction in metamaterials aimed to model the propagation of magneto-inductive waves in arrays of split ring resonators Syms et al. (2005), to analyze tunability through modifying lattice parameters Powell et al. (2010) and to analyze optical activity of twisted dimers Liu et al. (2007, 2012). Furthermore, several recent works propose analytical design approaches for metagratings based on rigorous descriptions of the mutual interaction of individual particles Chalabi et al. (2017); Ra’di et al. (2017); Rabinovich and Epstein (2018); Popov et al. (2018); Epstein and Rabinovich (2018). While these works on near-field coupling in metagratings provide full insight into the interaction of individual scatterers and demonstrate highly efficient devices, they have been largely restricted to reflective operation at large angles. More recently, near-field effects within and between cells of a metasurface have been investigated in gap plasmon structures Deshpande et al. (2018); Sharac et al. . For other important metasurface geometries, especially those for operation in transmission, near-field coupling is still largely unexplored and achieving the designed transmission response often requires black-box optimization Lavigne et al. (2018); Cole et al. (2018) or ad hoc iterative design methods Pfeiffer and Grbic (2013b).

In this paper, we demonstrate how the model of cascaded impedance sheets can be significantly inaccurate due to near-field coupling between metallic layers. We propose an extension of this model to incorporate near-field coupling and show that the extended model can predict the scattering parameters of metasurfaces accurately. This proposed coupling model does not require the explicit calculation of electromagnetic interaction integrals as in previous works Powell et al. (2010); Liu et al. (2012). The inaccuracies caused by near-field coupling can hamper the synthesis of transmissive multilayer metasurfaces. Therefore, we introduce an improved synthesis algorithm based on the proposed model, which provides good physical insight and reduced computational effort compared to black-box optimization. The algorithm is applied to a homogeneous Huygens metasurface with dogbone resonators, then to an inhomogeneous metasurface exhibiting anomalous refraction.

II Design of metasurfaces with three metallic layers

The model of cascaded impedance sheets is widely used in the literature for the analysis of transmissive metasurfaces Monticone et al. (2013); Epstein and Eleftheriades (2016b); Pfeiffer and Grbic (2013b). With this model, the properties of each layer are specified based on the scattering properties of the multi-layer stack, allowing the synthesis to be reduced to readily solvable problems. However, the model does not account for near-field coupling between layers, which can lead to significant inaccuracy. In the following, we briefly recapitulate the model and show that neglecting near-field coupling can lead to significant inaccuracy in the design of a homogeneous Huygens metasurface.

Each metallic layer is represented by a shunt impedance $Z_{n}$ and the dielectric layers are represented by a transmission line with propagation constant $\beta$, wave impedance $\eta=\sqrt{\mu/\epsilon}$ and thickness $t$ Tretyakov (2003) as depicted in Figure 1(a). Utilizing a circuit analogy, currents are equivalent to the magnetic field and voltages are equivalent to the electric field Tretyakov (2003). The transimssion line model accounts for far-field (or first order Floquet) coupling between the metallic-layers. It is difficult to calculate a meaningful propagation angle within a single meta-atom, therefore within each cell propagation normal to the surface is assumed Epstein and Eleftheriades (2016b). The corresponding circuit for the three layer system is shown in Figure 1(b). Using this circuit, the network parameters of each cell can be determined by cascading the individual building blocks in terms of transmission (ABCD) matrices $\mathbf{A_{tot}}=\mathbf{A_{Z1}}\mathbf{A_{tl}}\mathbf{A_{Z2}}\mathbf{A_{tl}}\mathbf{A_{Z3}}$, where the transmission matrix of a shunt impedance $\mathbf{A_{Zn}}$ and a transmission line $\mathbf{A_{tl}}$ are defined in Appendix A. The transmission matrix of the complete cell $\mathbf{A_{tot}}$ is transformed into a scattering matrix $\mathbf{S_{tot}}$. In Ref. Epstein and Eleftheriades (2016b) it was shown how the designed transmission response of the cell leads to specification for each of the three sheet impedances shown in Fig. 1. We perform numerical simulations of the dogbone structure to create a lookup table relating its geometric parameters to the sheet impedance $Z_{n}$, which is shown in Fig. 2. We note that the relationship between impedance and geometric parameters differs between the inner and outer layers due to the difference in surrounding dielectric medium.

To apply the model, we synthesize one unit cell of a Huygens metasurface for operation in the millimeter wave band at $f_{0}=$ 80 GHz with the metallic layers having the dogbone geometry depicted in Figure 1(c). With the smallest metallic feature size of 0.1 mm, a substrate thickness of $t=254\,\upmu$m and a permittivity of $\epsilon_{r}$=3.0, this unit cell design is compatible with standard printed circuit board fabrication while operating at W-band frequencies (similar to Ref. Pfeiffer and Grbic, 2013b). This results in relatively large resonators, with lateral unit cell dimensions of $a\gtrsim\lambda/5$. To allow comparison with theory, we assume lossless dielectric layers, and perfectly conducting metal patterns with zero thickness. Since Huygens metasurfaces are symmetrical, the first and third layer are identical, i.e. $Z_{1}=Z_{3}$, $G_{1}=G_{3}$ and $S_{1}=S_{3}$.

The metasurface is designed to reach near unity transmission amplitude and a transmission phase response of $\phi_{\mathrm{des}}=22.9^{\circ}$ at frequency $f_{0}$. In Figure 3, the transmission response of the model (black dashed line) is compared to results from full wave simulation (blue solid line). The model predicts that transmission is close to unity amplitude and meets the expected phase $\phi_{\mathrm{des}}$ at $f_{0}$. However, full wave simulation reveals that the transmission maximum and the corresponding phase are significantly shifted to 81.9 GHz. We attribute this discrepancy to near-field coupling (or equivalently higher order Floquet modes), which is to be expected, given that $\beta t\ll 1$. This clearly demonstrates that this widely-used model can fail for some geometries of practical interest, and consequently it is not precise enough for the synthesis of metasurfaces having high efficiency.

III Modelling near-field coupling between layers

Here we extend the model of metasurfaces with three metallic layers to incorporate near-field coupling between the layers. We note that the presented model is quite general and could also be applied to bianisotropic metasurfaces, or to designs containing four or more layers Abdo-Sanchez et al. (2018); Liu et al. (2016); Pan et al. (2018); Jiang et al. (2018).

III.1 Near-field impedance matrix

We extend the circuit model shown in Figure 1(b) by accounting for the mutual impedance between each of the metallic layers. As shown in Fig. 4, this is represented within the equivalent circuit model by current controlled voltage sources. To facilitate our analysis, we include a virtual port at the location of the centre metallic layer. Consistent with the existing model, we assume that propagation within a single meta-atom is normal to the surface.

We introduce the near-field impedance matrix $\mathbf{Z_{nf}}$ which gives the near-field relationship between the equivalent voltages $v_{n}$ and currents $i_{z,n}$ on each of the metallic layers

[TABLE]

The diagonal terms $Z_{n}$ are the self impedances of layer $n$ and the off diagonal elements are the transfer impedances $Z_{mut}^{lm}$ of the corresponding current controlled voltage sources, hence they represent coefficients that quantify the near-field coupling between layer $l$ and $m$. We consider only reciprocal structures, therefore $\mathbf{Z_{nf}}$ is a symmetric matrix with $Z_{mut}^{lm}=Z_{mut}^{ml}$. Also note that the use of a scalar current $i_{z,n}$ to represent the current flow on a metallic layer is valid so long as only one mode can be excited on it.

To relate $\mathbf{Z_{nf}}$ to the scattering parameters of the entire cell, we first convert it to admittance parameters $\mathbf{Y_{nf}}=\mathbf{Z_{nf}}^{-1}$. The admittance matrix of the complete cell $\mathbf{Y_{tot}}$ relates the equivalent current $i_{in,n}$ flowing into to each port to the equivalent voltages

[TABLE]

Admittance matrix $\mathbf{Y_{ff}}$ represents the far-field contribution of the transmission lines through their equivalent current $i_{tl,n}$.

[TABLE]

where we introduce the total transmission line current at the central node as $i_{tl,2}=i_{tl,2a}+i_{tl,2b}$. An analytical expression for $\mathbf{Y_{ff}}$ is given in Appendix A.

Using the relation $i_{in,n}=i_{z,n}+i_{tl,n}$ and Eq. (1) to (3), we can write

[TABLE]

Using these admittance matrices, near-field and far field contributions to the inter-layer coupling simply sum up and can be linked to the network parameters of the complete cell system $\mathbf{Y_{tot}}$.

As no current can flow through the virtual port 2 of Fig. 4, we impose the condition $i_{in,2}=0$. This leads to a reduced admittance matrix $\mathbf{Y_{tot}^{\prime}}$ for the complete cell, defined as

[TABLE]

with its relationship to $\mathbf{Y_{tot}}$ shown in Appendix A.

Given a near-field impedance matrix $\mathbf{Z_{nf}}$, the reduced two port admittance matrix $\mathbf{Y_{tot}^{\prime}}$ can be determined using equation Eqs. (4) and (11). This is then converted to a scattering matrix $\mathbf{S_{tot}}$, to obtain transmission and reflection properties of the metasurface. Note that if the near-field coupling coefficients $Z_{mut}^{lm}$ are set to zero, this formalism is equivalent to the existing model outlined in Section II.

III.2 Determination of the near-field impedance terms

To apply the formalism developed in Subsection III.1 to metasurface analysis, it is necessary to find the values of elements of $\mathbf{Z_{nf}}$. Since $\mathbf{Z_{nf}}$ is a 3x3 matrix and $\mathbf{S_{tot}}$ is a 2x2 matrix, a single far-field scattering simulation yields insufficient information to resolve $\mathbf{Z_{nf}}$. Therefore, we propose a technique to determine the elements of $\mathbf{Z_{nf}}$ successively from full wave simulations of reduced systems containing only two layers. A two layer formalism analogous to Section III.1 is applied to these full-wave scattering parameters, which allows the unambiguous extraction of individual elements of $\mathbf{Z_{nf}}$. This two layer formalism is outlined in Appendix A.

A typical configuration for a full wave simulation of the unit cell of a metasurface with three layers is shown in Fig. 5 (a). The unit cell is centered at $z=0$, periodic boundary conditions apply to all sides and Floquet ports are located at a large distance $d_{p}$ to the top and bottom of the structure, where the influence of high order Floquet modes is negligible. The reference plane of these ports however needs to be moved back to the top and bottom to $z=\pm t$, respectively, to correct the reflection and transmission phase Epstein and Eleftheriades (2016a).

By removing the central metallic layer, as shown in Fig. 5(b), the two layer formalism can be applied to determine the self impedance of layers one and three, $Z_{1}$ and $Z_{3}$, and their mutual impedance $Z_{mut}^{13}$. To determine $Z_{1}$, $Z_{2}$ and $Z_{mut}^{12}$, the configuration from Fig. 5(c) is utilized. Here, however, the reference plane of Port 2 needs to be moved to the central metal layer, which can be done by inverting the ABCD matrix of the dielectric layer Pozar (2011). Similarly, $Z_{2}$, $Z_{3}$ and $Z_{mut}^{23}$ can be obtained from the configuration from Fig. 5 (d). Using this successive approach, the full coupling matrix $\mathbf{Z_{nf}}$ is determined and the three layer system is fully characterized.

We apply this technique to the metasurface considered in Section II to determine the near-field impedance matrix $\mathbf{Z_{nf}}$. As we consider a Huygens metasurface with $Z_{1}=Z_{3}$ and $Z_{mut}^{lm}=Z_{mut}^{ml}$, only two simulations are required to determine the near-field impedance matrix $\mathbf{Z_{nf}}$, since the configurations in Figure 5(c) and (d) are equivalent. The resulting self impedances $Z_{n}$ and mutual impedances $Z_{mut}^{lm}$ are shown in Figure 6. As we consider lossless structures, these functions are purely imaginary to within numerical accuracy ($|\mathrm{Re}(Z)|/\eta_{0}<10^{-5}$). At the design frequency, the mutual impedance $Z_{\mathrm{mut}}^{12}$ and $Z_{mut}^{13}$ are clearly smaller than the self impedance $Z_{1}$ and $Z_{2}$, however they are not negligible. The mutual and self impedance terms have the form of a series LC circuit. Using the calculated near-field impedance matrix, the metasurface transmission is plotted in Figure 3 (red dashed line), showing good agreement with the full wave simulation. This result confirms the accuracy of our near-field coupling model, demonstrating why the conventional metasurface synthesis procedure fails in this case.

IV Metasurface synthesis incorporating near-field coupling

IV.1 Refinement algorithm for metasurface synthesis

Having characterized how near-field coupling can render existing metasurface design techniques inaccurate, we present an approach to include it within the metasurface synthesis. Since the near-field coupling cannot be evaluated until a geometry of the metal layer has been specified, this necessarily involves the use of an iterative algorithm. The goal of the algorithm is to find the geometrical parameters of the unit cell that match the target scattering parameters. In the case of Huygens metasurfaces, it is sufficient to minimize the difference between the realized transmission $S_{21}^{\mathrm{re}}$ and designed transmission $S_{21}^{\mathrm{des}}$

[TABLE]

The realized transmission $S_{21}^{\mathrm{re}}$ is calculated from our near-field coupling model with the matrix $\mathbf{Z_{nf}}$. A fundamental assumption of the algorithm is that for small changes of self impedances $Z_{n}$, the geometric change is small, and the mutual coupling does not change strongly. This allows us to correct only the diagonal elements of the coupling matrix with a factor $\gamma$ at each iteration

[TABLE]

where the iteration index is denoted with $p$. The factor $\gamma$ is chosen by minimizing the term $|S_{12}^{\mathrm{re}}(\mathbf{Z_{nf,p}})-S_{12}^{\mathrm{des}}|^{2}$ numerically in each iteration.

Based on the corrected self impedances $Z_{n,p}$, the geometry of the unit cell is updated with new geometrical parameters $G_{n}$ and $S_{n}$ using the lookup table in Figure 2. If a further iteration is necessary, the mutual impedances are recalculated with full wave simulations of the new geometry. Note that Eq. (7) could be modified to adjust outer layer impedances separately from the inner layer ones, however this was found empirically to offer no improvement in convergence. The optimization process is outlined in full in Algorithm 1.

The algorithm terminates when the relative change of the self impedances from one iteration to the next is smaller than $\gamma_{c}$. This convergence parameter was chosen as $\gamma_{c}=$0.085 for the examples shown in this work. The absolute value of the resulting correction factor for the self impedances $|\gamma|$ ranged between 0.038 and 0.25. The ability of this algorithm to correct the transmission response of a cell is shown in Fig. 7. The red dashed line shows the initial response according to the design procedure outlined in Section II. The black dash-dot curve shows the result of the model after the optimization, and the blue lines show the corresponding results of full-wave simulation. It can be seen that the strong perturbation of the operating frequency has been corrected, and the metasurface now achieves complete transmission and the designed phase value at the operating frequency.

IV.2 Application to metasurface demonstrating anomalous refraction

We demonstrate the further utility of our iterative algorithm by applying it to an inhomogeneous Huygens metasurface demonstrating anomalous refraction of a plane wave. The incident angle was chosen $\theta_{\mathrm{in}}=\,0$ and the refraction angle $\theta_{\mathrm{out}}=\,55^{\circ}$, leading to a supercell width $d=4.56\,$mm. Each cell is based on the dogbone geometry shown in Fig. 1(a). To minimize the influence of coupling to non-identical neighboring cells which is not described in our model, the lateral spacing $a_{x}$ was chosen as large as possible. With the supercell width used here, there are three propagating Floquet modes in transmission, $T_{0}$ and $T_{\pm 1}$, and three in reflection $R_{0}$ and $R_{\pm 1}$. To control these 6 degrees of freedom, a minimum of three cells per supercell is required, each having an engineered electric and magnetic response.

The initial sheet impedances $Z_{n,p=0}$ and the required transmission response for each cell $S_{21}^{\mathrm{des}}$ were specified according to the design procedure in Epstein and Eleftheriades (2016a). With these target values, impedance matching to the refracted wave is stipulated. The geometry of each cell was synthesized using the algorithm from Section IV.1. The self impedance $Z_{n,p}$ of each cell $n$ at iteration $p$ is shown in Figure 8. In the case of the first cell, having identical geometry to that was analyzed in Section II, three iterations were necessary to obtain the required scattering properties. For the other two cells, one iteration was enough, due to the lesser influence of near-field coupling in these cases. The resulting geometrical parameters of the supercell are given in Table 1.

We note that an alternative approach to correct for near-field coupling is to use black-box optimization methods, which are built in to many commercial simulation packages. Compared to these methods, our proposed algorithm has several advantages. Unlike these black-box optimization methods, our approach gives clear physical insight into the near-field coupling which causes the metasurface design to be non-optimal. In particular, with the quantities $Z_{nm}^{mut}$, near-field coupling effects can be quantified and monitored from iteration to iteration. Additionally, our algorithm is much more computationally efficient, as is demonstrated by the results in Appendix B, which demonstrates a speed up factor of 3.6 for the design considered. For metasurface designs with a large number of cells, this speed up would also be a significant benefit. Furthermore, using the proposed algorithm, the influence of near-field coupling effects can be corrected without any prior knowledge. Applying black-box optimization methods requires guessing the required geometric parameter changes and setting suitable parameter boundaries.

While we considered only a single polarization in this work, we expect that our method could be adapted to tensorial impedances Selvanayagam and Eleftheriades (2016) used to model polarization rotating metasurfaces. We note that for chiral metasurfaces based on coupled high-symmetry layers Hannam et al. (2014), polarization rotation only occurs when near-field coupling is strong, thus near-field coupling cannot be neglected in equivalent circuit modelling.

The refraction performance of the synthesized metasurface is shown in Figure 9, where $T_{1}$, the transmission coefficient into the designed Floquet mode, should be maximized, and other terms minimized. The performance of the initial geometry (a) is compared with the performance of the geometry obtained by our iterative optimization process (b). For the initial metasurface, the energy refracted into the desired Floquet mode is as low as 56 % at the design frequency $f_{0}=80\,$GHz. The maximum refraction efficiency is 78 %, but it is significantly shifted to 81.4 GHz. The optimized metasurface on the other hand refracts 82 % into the desired Floquet mode $T_{1}$ at the design frequency. In this case, the maximum refraction efficiency is 90 % at 79.7 GHz. We note that the maximum efficiency of the optimized metasurface is close to the theoretical efficiency limit Epstein and Eleftheriades (2016a) for Huygens metasurfaces, which is 92 % for a refraction angle of $55\,^{\circ}$.

V Conclusion

We showed that metasurfaces based on patterned metallic layers separated by dielectric layers can exhibit strong inter-layer near-field coupling. This can lead to significant inaccuracy of existing models used for metasurface synthesis, resulting in poor performance of the metasurfaces. We developed an improved model which accounts for near-field coupling and showed that it closely matches the results of full-wave simulation. This is of particular significance for higher frequency designs, where fabrication tolerances mean that feature sizes are not highly sub-wavelength, leading to strong near-field coupling.

To mitigate the performance degradation induced by near-field coupling, we introduced an iterative algorithm for accurate metasurface synthesis which takes it into account. The presented algorithm gives more physical insight into the problem than black-box optimization, and reduces computation time significantly. Additionally, it is very generic and can be applied to multilayer metasurfaces in various frequency ranges. As this work introduces an efficient synthesis method for such geometries, it helps to relax fabrication requirements and pave avenues for new metasurface geometries.

Acknowledgements.

This work was financially supported by the Australian Research Council (Linkage Project LP160100253) and the Luxembourg Ministry of the Economy (grant CVN 18/17/RED).

Appendix A Two and three-port matrix representations

The ABCD matrices of a shunt impedance $\mathbf{A_{Zn}}$ and transmission line $\mathbf{A_{tl}}$ Pozar (2011) are given by

[TABLE]

and

[TABLE]

The admittance matrix $\mathbf{Y_{ff}}$ for the three layer system can be determined analyzing the circuit from Figure 4 (b) without shunt elements and dependent sources utilizing the definition of an admittance matrix element ${Y_{ff,lm}=\frac{i_{tl,l}}{v_{m}}\left.\right|_{v_{k\neq m}=0}}$ and states

[TABLE]

To reduce the three port equivalent circuit shown in Fig. 4 to a two-port network compatible with circuit simulations, we impose the condition $i_{in,2}=0$. This yields the following expression for the two port matrix $\mathbf{Y_{tot}^{\prime}}$ in terms of elements of the three port matrix $\mathbf{Y_{tot}}$

[TABLE]

Equations (1) to (4) can be defined analogously for the two layer system using the 2x2 matrices $\mathbf{Y_{tot}^{2}}$,$\mathbf{Y_{nf}^{2}}$ and $\mathbf{Y_{ff}^{2}}$. In the corresponding circuit model in Figure 10, the admittance matrix for the transmission line is

[TABLE]

With known admittance matrix $\mathbf{Y_{tot}^{2}}$, we can determine $\mathbf{Z_{nf}^{2}}=\mathbf{Y_{nf}^{2}}^{-1}$ without any ambiguity. The coefficients of $\mathbf{Z_{nf}^{2}}$ for configuration (b) to (d) in Figure 5 are then used subsequently to build the 3x3 matrix $\mathbf{Z_{nf}}$. Note, that in configuration (b), the thickness is $2t$.

Appendix B Computational efficiency

To evaluate the computation time of the algorithm proposed in this work, it was compared with a black-box trusted region framework algorithm implemented in the commercial package CST Microwave Studio Studio (2018). The merit function for the optimization is $f_{m}=|S_{12}^{\mathrm{re}}-S_{12}^{\mathrm{des}}|^{2}$, as previously. By far the most significant contribution to the computation time comes from full wave simulations. Therefore, we compare here the number of full wave simulations $n_{fw}$ that is required until the respective optimization algorithm converges for both of these methods. Since two full wave simulations are required to determine the near-field impedance matrix $\mathbf{Z_{nf}}$, the proposed model based optimization requires $n_{fw}=2p_{t}$ when terminating at iteration $p_{t}$. The black-box optimization on the other hand requires $n_{fw}=p_{t}+1$.

In Figure 11, the computational effort for the two algorithms is compared in terms of $n_{fw}$, for each of the three cells in our anomalously refracting metasurface. Our proposed optimization technique was on average 3.6 times faster than the black-box optimization. Additionally, the trust region optimization requires some prior knowledge to set the boundaries of the parameter space. Here, we chose $\pm 25\%$ on the initial parameter set.