Theoretical and experimental exploration of finite sample size effects on the propagation of surface waves supported by slot arrays

PHYSICAL REVIEW B 95, 245425 (2017) Theoretical and experimental exploration of finite sample size effects on the propagation of surface waves supported by slot arrays Miguel Camacho,1,*Rafael R. Boix,2Francisco Medina,2Alastair P. Hibbins,1and J. Roy Sambles1 1Electromagnetic and Acoustic Materials Group, Department of Physics and Astronomy, University of Exeter, Stocker Road, Exeter EX4 4QL, United Kingdom 2Department of Electronics and Electromagnetics, College of Physics, Av. Reina Mercedes s/n, 41012 Seville, Spain (Received 12 February 2017; revised manuscript received 5 June 2017; published 22 June 2017) The propagation of surface waves supported by a finite array of slots perforated on a zero thickness perfect electrically conducting screen is studied both experimentally and theoretically. To generate numerical results, the integral equation satisfied by the electric field in the slots is efficiently solved by means of Galerkin’s method, treating the metal as perfectly conducting. The finite size of the array along the direction of propagation creates a family of states of higher momentum and lower amplitude than the single mode for the corresponding infinite array. These modes are spaced in momentum with a periodicity inversely proportional to the length of the array. In addition, the finite width in the transverse direction produces a set of higher frequency modes due to this additional quantization. Both effects arising from finite sample length and width are explained by the theoretical model and validated experimentally. DOI: 10.1103/PhysRevB.95.245425 The existence of strongly localized surface modes supported by corrugated metallic surfaces at frequencies varying from terahertz [1]tomicrowaves[2] is well established. At optical frequencies, due to the Drude-like properties of the metal, these exist even on flat surfaces, being labeled surface plasmons polaritons (SPPs) [3,4]. However, following the discovery of extraordinary optical transmission (EOT) [5] and the later theoretical calculations, these surface waves were linked to a massively enhanced transmission through subwavelength holes [6,7] due to the diffractive coupling of radiative modes with the SPPs. This theory was, however, essentially developed for optical frequencies relying on the presence of SPPs. The subsequent discovery of EOT at frequencies at which metals behave as near perfect electric conductors (PEC), such as microwaves [8,9], launched the search for a more general theory. Many workers in the field made analogies of the surface waves supported on a patterned PEC with the aforementioned SPP, and due to the similarities with the SPP dispersion, the modes became commonly known as the “spoof surface plasmons” [10–12]. It should be noted, however, that there has been a wealth of studies of surface waves on metals at low frequencies, dating back to the 1940s [13,14]. This connection between the surface waves supported by hole arrays and the transmission through them has given rise to much research [15,16] in which the transmission through both infinite and finite arrays of holes and slits has been studied often using a modal matching approach [17,18]. However, there is no corresponding study on the effects of finite sample size on the propagation of surface waves themselves. At microwave frequencies, these surface modes have been explored particularly using metasurfaces [19], subwavelengthstructured metallic surfaces that yield novel properties such as negative refraction [20], band gaps [21], and dispersive behavior leading to the design of surface wave lenses and antennas [22–25]. With the development of miniaturization *[email protected] techniques, some of these applications have also been extended to optical frequencies [26,27]. For the design of such surface wave controlling structures, periodic boundary conditions are assumed therefore reducing the analysis to a single unit cell. In reality, however, infinite arrays do not exist, leading to the necessity of characterizing a finite system assuming that it has the same characteristics of an infinite array. In reality, when experimentally measuring the dispersion relation of the surface waves supported by an array of slots, effects due to the finite dimensions along the different directions of a rectangular array cannot be avoided. It is just these effects that are presented here and numerically reproduced by using a method of moments (MoM) technique, which is based on a previous study of some of the authors [28]. To be able to compare the behavior of infinite and finite arrays this method is first applied to the analysis of infinite periodic arrays, and the results obtained compared to those from a commercial software based on finite element method. Then, the results for finite arrays with MoM are presented and compared to the experimental results, confirming the effects associated with the finite sample size. I. INFINITE PERIODIC ARRAY Let us start by considering the propagation of surface waves supported by an infinite periodic array of slots. In Fig. 1,the unit cell of dimensions a×bis shown, containing a negligible thickness PEC screen perforated with a slot of dimensions ls×ws. This lossless approximation is accurate in the microwave regime up to millimeter-wave frequencies. We first derive the integral equation that the electric field in the slot needs to satisfy, for which a time dependence given by ejωt will be assumed and suppressed from this point forward. Let us define the in-plane component of the electric field as Et(x,y,z =0) =Ex(x,y,z =0)ˆ x+Ey(x,y,z =0)ˆ y.The electric field at any point of the surface will produce a differential current at any other point, which can be calculated in the form of a dyadic Green’s function. However, the total electric current obtained by the sum of all those differential 2469-9950/2017/95(24)/245425(7) 245425-1 ©2017 American Physical Society CAMACHO, BOIX, MEDINA, HIBBINS, AND SAMBLES PHYSICAL REVIEW B 95, 245425 (2017) FIG. 1. Infinite periodic array of slots perforated in a negligible thickness perfectly conducting screen. The periodicities of the array along the xand ydirections are aand b, respectively, and lsand ws are the dimensions of the slots. currents needs to vanish where there is no metal, i.e., at the surface occupied by the slots. This reasoning can be expressed mathematically as the following integral equation ∞ −∞ ∞ −∞ GM(x−x,y −y)·Et(x,y,z =0) dxdy =0(x,y)∈slots.(1) Here GM(x,y) is the dyadic Green’s function relating the surface electric current density on the conducting screen and the tangential electric field in the slots. The double integral has to be extended to all the slots in the infinite periodic structure. The dyadic Green’s function can be obtained as GM(x,y)=k2 0+∂2 ∂y2g(x,y)−∂2g(x,y) ∂x∂y −∂2g(x,y) ∂x∂y k2 0+∂2 ∂x2g(x,y)(2) where g(x,y)=− je−jk0√x2+y2 πk0Z0x2+y2,(3) where k0and Z0represent the free space wave vector and impedance, respectively. Due to the periodicity of the structure, the limits of this integral can be reduced, applying Bloch’s theorem, to a single unit cell by defining a periodic Green’s function given by Gper M(x,y)=+∞  m,n=−∞ GM(x−ma,y −nb)ej(kx0x+ky0y).(4) In this equation the phase factor accounts for the propagation of a surface wave mode along the periodically perforated screen and kx0and ky0are the components of the momentum along the xand ydirections, respectively (kx0awould represent the phase shift between the centers of two adjacent cells in the x direction). When (4) is substituted in (1), an integral equation can be obtained given by a 0b 0 Gper M(x−x,y −y)·Et(x,y,z =0) dxdy =0(x,y)∈δ00,(5) where δ00 is the slot {(a−ws)/2<x<(a+ws)/2; (b−ls)/2<y<(b+ls)/2}of the periodic cell C00 covering the rectangular domain {0⩽x⩽a;0⩽y⩽b}. To convert this integral equation into a system of linear equations we make use of Garlerkin’s version of the MoM. First, the electric field at the surface of the slots is approximately expanded as a linear combination of known basis functions, this is, Et(x,y,z =0) ≈ Nb  j=1 e∞,j bj(x,y)(x,y)∈δ00,(6) where e∞,j are unknown constant coefficients. Note that, although in principle only an infinite number of basis functions would be an exact solution, very good convergence can be obtained when these are chosen adequately to the geometry of the problem (in this case, these will be given by Chebyshev polynomials multiplied by the edge behavior of the electric field for each polarization) [28,29]. By substituting (6)in(5) and using each of the basis functions as weighting and then projecting that expression into the domain C00 the following system of equations is obtained Nb  j=1 ij e∞,j =0(i=1,...,N b)(7) where ij =a 0b 0 b∗ i(x,y)·a 0b 0 Gper M(x−x,y −y) ·bj(x,y)dxdydxdy. (8) These matrix elements can be calculated very efficiently by making use of the properties of the Fourier transform, which lead to quasianalytical expressions, as shown elsewhere [28]. With the addition of the in-plane momentum components (kx0and ky0), the calculation of the periodic dyadic Green’s function needs to be generalized but can still be efficiently calculated by means of Ewald’s method [30,31]. Once the matrix ij is calculated, the dispersion relation is obtained by imposing the determinant of ij to be equal to zero. Figure 2 shows the dispersion curve along the xdirection (ky0=0) obtained from (7) and (8) compared to the results predicted by a commercial finite element method (FEM) solver [32] for different lengths of the slot ls(varying from 6 to 9 mm). The good agreement shown validates both the method and the chosen set of basis functions. It can be observed that, as the length of the slot increases, the wave vector between the surface mode and the radiative zone (the zone for which 2πf > ckx0being fthe frequency and cthe speed of light) increases, leading to a more rapid exponential decay of the fields along the zdirection. It is interesting to note that, when the fields in the slots are excited by incident plane waves, the resonance (associated with maximum transmission for 245425-2 THEORETICAL AND EXPERIMENTAL EXPLORATION OF . . . PHYSICAL REVIEW B 95, 245425 (2017) FIG. 2. Predicted dispersion obtained using the presented MoM implementation (lines) compared to FEM models (markers) for different slot lengths ls. The gray area corresponds to the nonbound regime (above the light line). The values for the other dimensions are the same for every curve, with ws/a =0.05 and a=b=10 mm. dipole/slot elements) [33] is located near the frequency at which the length of the slot corresponds to half the wavelength. However, here we find that the surface mode is supported at frequencies well below that resonance. For the case of ls=9 mm, the expected transmission peak would be near 16.7 GHz, while we find that the bound surface mode only exists below 12 GHz. We have found numerically, however, that the same set of basis functions can reliably explain the behavior of slot elements for both radiative (as shown elsewhere [28]) and bound regimes (with weights e∞,j that will vary from one problem to the other). II. TRUNCATED PERIODIC ARRAY Let us now consider the propagation of surface modes along a finite array of slots. For this purpose, we will consider an infinite (negligible thickness) metal layer perforated with a truncated periodic array of M=Nx×Nyslots with constant separation given by aand bin the xand ydirections, respectively, as that shown in Fig. 3. Due to the lack of strict periodicity, Bloch’s theorem is no longer applicable (which means that the electric field distribution varies both in amplitude and phase from one cell to another) and therefore the study cannot be reduced to a single cell by assuming a phase shift in the fields between two adjacent unit cells. Therefore, no periodic Green’s function can be defined, and we need to use a different approach to study the propagation of surface waves. From a physical point of view this means that an excitation needs to be included into our model. Here we replicate one method used experimentally, by using a localized current density near the central element of the array. This can be easily implemented as an excitation term when deriving the integral equation that the electric field needs to satisfy on the perforated surface. Following the reasoning used for the infinite array, the electric field on the surface needs to satisfy the following M FIG. 3. Truncated periodic array of slots perforated in an infinite negligible thickness perfectly conducting screen. The periodicity of the array along the xand ydirections are aand b, respectively, and lsand wsare the dimensions of the slots. coupled integral equations Jas(x,y)+ M  j=1ηj GM(x−x,y −y) ·Esc t(x,y,z =0) dxdy=0(x,y)∈ηi (i=1,...,M),(9) where GM(x,y) is the dyadic Green’s function defined in (2) and (3). Since in the case of a truncated array of slots Esc t(x,y,z =0) is not a periodic function of xand y,theset of integral equations (9) cannot be reduced to one single unit cell as for the infinite array as shown in (5). Note that in this case, Esc t(x,y,z =0) represents only the field scattered by the slots when illuminated with an excitation in the form of a position-dependent surface current given by Jas(x,y)inthe absence of slots. To determine the value of Esc t(x,y,z =0) in ηj(j=1,...,M), we will apply the MoM to each of the M integral equations of (9). That is, the tangential electric field in the jth slot will be approximately expressed in terms of the proposed basis functions (electric field profiles) djl(x,y)= djl,x(x,y)ˆ x+djl,y(x,y)ˆ y(l=1,...,N b) as shown below Esc t(x,y,z =0) ≈ Nb  l=1 ejldjl(x,y)(x,y)∈ηj.(10) In this case, a nonhomogeneous linear system of equations can be derived, given by M  j=1 Nb  l=1 kl ij ejl =pik (i=1,...,M;k=1,...,N b),(11) 245425-3 CAMACHO, BOIX, MEDINA, HIBBINS, AND SAMBLES PHYSICAL REVIEW B 95, 245425 (2017) where kl ij =ηi d∗ ik(x,y)·ηj GM(x−x,y −y) ·djl(x,y)dxdydxdy (i,j =1,...,M;k,l =1,...,N b) (12) and where pik =− ηi d∗ ik(x,y)·Jas(x,y)dxdy (i=1,...,M;k=1,...,N b).(13) In order to model the excitation of a very localized source, we have made use of a delta-gap excitation [34], which replicates an infinitely thin wire connecting the two parallel vertical edges of the slot. For symmetry reasons it has been chosen to be located in the middle xplane of the slot, but it could be placed elsewhere. Mathematically, this can be easily introduced into the model as Jas(x,y)∝δ(y−ycc)for (xcc −ws/2) <x<(xcc +ws/2) and zero otherwise (where ycc and xcc are the yand xcoordinates of the center of the central slot, respectively). In a very similar manner to the case of an infinite array, the matrix elements given by (12) can be efficiently calculated when rearranged in terms of cross correlations of the basis functions as explained in detail elsewhere [28]. Once these coefficients are calculated, the electric field is completely determined on the surface. Following the experimental techniques used for the study of surface waves [20,35], the dispersion diagram is calculated by computing the Fourier transform of the field distribution into the spectral domain (which is continuous given the lack of periodicity). Taking advantage of the analytical basis functions used to expand the electric field on the slots, we can compute the analytical Fourier transform. This makes it possible to obtain as small a resolution as needed in the Fourier transformed fields, in contrast to the fast Fourier transform algorithm used experimentally. To test the convergence of our MoM implementation for a finite array by comparison to the solution obtained for an infinite array of slots (shown in Fig. 2), we have initially modeled an array of 49 ×49 cells (this gives good definition of the modes in kspace in comparison to that obtained for smaller arrays within a reasonable computational time) with ls=9mm,ws=0.5mm,a=b=10 mm as previously studied for the infinite case. The obtained dispersion diagram is shown in Fig. 4, where the brightest feature can be directly compared to the dispersion of the bound mode of the infinite array identified in Fig. 2(purple line). This is not, however, the only mode supported by the truncated array: A family of bound modes can be found at higher frequencies within the light cone above the limit of the infinite array dispersion. In addition, each of the modes appears to be periodically repeated in kx(this can be easily seen outside of the cone formed by the two light lines). In the following we will explore the origin of these effects and how the dimensions of the array can be optimized to reduce them. Although the numerical results for an array of 49 ×49 elements show good agreement with those of an infinite array, the FIG. 4. Electric field Fourier transform predicted by MoM for a range of frequencies for an array of Nx=Ny=49 slots for which ls/a =0.9, ws/a =0.05, and a=b=10 mm. It is represented linearly in color as the square root of the absolute value of the Fourier amplitude. The white superposed squares correspond to the dispersion diagram obtained for an infinite periodic array with the same unit cell. validation of the truncation effects is most practically achieved experimentally, as full-wave numerical solvers require a significant computational effort to solve large but finite arrays. In addition, the large number of higher order modes observed (Fig. 4) would also make experimental validation challenging since their position and number will also depend on the symmetry of the sample. For these reasons, we have studied the dependence of these higher order modes on the size of the array in the xand ydirections. We discover that the number of higher frequency modes is only affected by the length of the array in the direction perpendicular to that of the propagation of the surface wave (Nybas we are studying propagation in x), with the number of modes increasing with the number of elements. In contrast, the spacing in kxbetween two of the translated modes is only affected by the length of the array along the propagation direction (Nxain our case), the separation between the modes being reduced as the length increases. In order to simplify our analysis, we have chosen to experimentally study the case where Nx=39 and Ny=5. Figure 5shows the results obtained by means of MoM, where FIG. 5. Electric field Fourier transform predicted by MoM in a range of frequencies for an array of Nx=39 and Ny=5 slots for which ls/a =0.9, ws/a =0.05, and a=b=10 mm. 245425-4 THEORETICAL AND EXPERIMENTAL EXPLORATION OF . . . PHYSICAL REVIEW B 95, 245425 (2017) FIG. 6. Experimental electric field Fourier transform in a range of frequencies for an array of Nx=39 and Ny=5 slots for which ls/a =0.9, ws/a =0.05, and a=b=10 mm. three modes are clearly distinguished and the same periodic repetition of the modes along the kxaxis as for the case of Nx=Ny=49 is present. The experimental sample has been prepared by making use of a very thin plastic layer (50 μm with r=2.8) coated with 17 μm thick copper, in which slots have been etched. By making use of a pair of probe antennas (made of a coaxial cable with the inner conductor exposed), the electric field is excited at the center of the array and scanned at a constant height across the surface [20,35]. This electric field distribution is then Fourier transformed from the spatial domain into the spectral domain making use of a fast Fourier transform algorithm for each frequency, obtaining the dispersion diagram shown in Fig. 6. A small frequency redshift is present with respect to the numerically obtained dispersion as a consequence of the introduction of the thin dielectric layer. This shift, however, does not change any of the physics involved, and the same truncation effects can be observed in both figures. Note that in the experimental data, there is a decrease in Fourier amplitude as kxincreases. This is because the source antenna does not provide an ideal wave-vector (momentum) spectrum. By investigating this simpler case (with 39 ×5 slots), and thereby reducing the complexity of the mode structure compared to the case of 49 ×49 slots, we are now in a position to better understand the origin of these additional modes. In Fig. 7the electric field magnitude and phase distribution at the surface of the array have been represented. The first and third color maps correspond to the magnitude and phase of the lowest frequency mode (excited at 12.19 GHz, corresponding to kxa/π =1) and show how the fields in every slot in each column responds in phase to the excitation while two consecutive rows are out of phase, in concordance with the mode being at the Brillouin zone boundary. In contrast, the second and fourth color maps show the electric field magnitude and phase distribution for the second lowest frequency mode (excited at 13.54 GHz). In this case we observe that two of the rows are not excited and that the slots on the top and bottom of each column are excited out of phase with respect to the slot in the central row. This shows that although we are studying propagation along the xdirection, the absence of infinite periodicity allows for the excitation of modes with a FIG. 7. MoM predictions of the electric field distribution for the two lowest frequency modes supported by a finite array of slots with Nx=39 and Ny=5 slots for which ls/a =0.9, ws/a =0.05, and a=b=10 mm excited at the central slot with a delta gap excitation. From the top, first and third figures correspond to the magnitude and phase of the electric field at 12.19 GHz while the second and fourth correspond to the magnitude and phase at a frequency of 13.54 GHz, respectively. structure in the direction perpendicular to that of propagation. These occur at higher frequencies than for the most similar mode to that supported by infinite arrays. Although not shown here, the third mode that appears in Fig. 5corresponds to the excitation of all five rows with each element being out of phase with respect to its closest neighbors, leading to the highest frequency of the modes shown. This can be easily predicted using simple reasoning: Each of the higher order modes will require an extra (transverse) momentum ktrans yto excite such phase differences between the rows, which needs to be of the order of ktrans y≈π/(2b) in the case of the second lowest frequency mode. This predicts, at the Brillouin zone boundary along the xdirection, a frequency separation of ≈1.36 GHz which is close to the 1.67 GHz separation between the two lowest order modes found in Fig. 5. This discrepancy can be explained by edge effects for the top and bottom rows, at which the excited fields of these higher order resonances decay over a short distance due to the presence of the semi-infinite conducting plane. We finally analyze the phenomenon of the multiple mode dispersion curves that are translated in kxevident in Figs. 4, 5, and 6. The evolution of the spacing between such repeated modes has been analyzed for different array sizes and is only affected by the length of the array along the direction of propagation (Nxain this case). In Fig. 8, data is plotted having been extracted from different dispersion diagrams along the Brillouin zone boundary (kxπ/a =1) for different values of Nxwhen keeping Ny=5. It can be seen how the frequency separation of two consecutive modes reduces with the length of the array (Nxa), with corresponding narrowing of the original mode (infinite array). This separation in frequency is associated with a change in wavelength λ ≈Nxaand a change in wave vector kx≈π/Nxa. This quantization is associated with the boundary condition imposed by the ends of the array that dictates that the mode has zero electric field amplitude. This effect cannot be found, however, for the third 245425-5 CAMACHO, BOIX, MEDINA, HIBBINS, AND SAMBLES PHYSICAL REVIEW B 95, 245425 (2017) FIG. 8. Normalized amplitude of the Fourier transform of the MoM prediction of the electric field distribution at the Brillouin zone boundary along the kxdirection (kx=π/a) for finite arrays of slots of different lengths when excited at the central slot. All three cases have been calculated for Ny=5, ls/a =0.9, ws/a =0.05, and a= b=10 mm. nonbound mode, whose amplitude decays with distance from the center to negligible values at the edges of the sample and hence little scattering occurs. Consequently, no such repeated modes can be observed in Fig. 5above 14 GHz. Note that only modes with nonzero amplitude in the center of the array can be excited since this is where the source is positioned. From the results presented here, one could think of the modes supported by the finite array as very similar to those of a two-dimensional cavity. However, the fact that it is an open (infinite) system containing a finite array means that its field distribution can only be decomposed in terms of a continuous momentum spectrum. This can be understood from the need of a Fourier transform (instead of Fourier series) to find the modes supported by the structure from the decomposition of the electric field on the surface. This continuous momentum spectrum allows the scattering of the surface mode (even if it was the only component excited by an ideal source) at the end of the array into both radiative and nonradiative components. Although some of the latter may excite the mode associated with the infinite periodic array, it will also excite components not allowed there, as we could see in Figs. 4,5, and 6. III. CONCLUSION In conclusion, the propagation of surface waves supported by both infinite and finite periodic arrays of slots has been investigated. To do this we have made use of an efficient and accurate approach based on the method of moments to solve the integral equation satisfied by the component of electric field tangential to the surface of the rectangular holes for both infinite and finite problems. Predictions from this model have been verified by means of FEM modeling in the case of infinite periodic arrays for different geometries and experimentally for finite arrays. For finite systems, two additional effects have been observed, separately related to the size of the array along the propagation direction and the transverse direction, respectively. By using the field distributions predicted by the MoM, the existence of higher frequency modes not present in the infinite case has been explained in terms of quantization arising from the finite width of the system which can quasiquantitatively predict this frequency splitting. In addition, the periodic repetition of the modes with decreasing amplitude along the kxaxis has been explained by the backscattering of the surface waves by the ends of the array, with spacing in the wave vector controlled by the length of the array. ACKNOWLEDGMENTS The authors wish to acknowledge financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the United Kingdom, via the EPSRC Centre for Doctoral Training in Metamaterials (Grant No. EP/L015331/1). R.R.B. and F.M. acknowledge financial support from the Spanish MINECO (Grant No. TEC2013-41913-P) and the Spanish Junta de Andalucía (Grant No. P12-TIC-1435). All data created during this research are openly available from the University of Exeter’s institutional repository at https://ore.exeter.ac.uk/. [1] R. Ulrich and M. Tacke, Appl. Phys. Lett. 22,251 (1973). [2] L. Goldstone and A. Oliner, IRE Trans. Antennas Propag. 7, 274 (1959). [3] R. H. Ritchie, Phys. Rev. 106,874 (1957). [4] W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature (London) 424,824 (2003). [5] T. W. Ebbesen, H. J. Lezec, H. F. Ghaemi, T. Thio, and P. A. Wolff, Nature (London) 391,667 (1998). [6] L. Martín-Moreno, F. J. García-Vidal, H. J. Lezec, K. M. Pellerin, T. Thio, J. B. Pendry, and T. W. Ebbesen, Phys. Rev. Lett. 86,1114 (2001). [7] C. Genet and T. W. Ebbesen, Nature (London) 445,39 (2007). [8] M. Beruete, M. Sorolla, I. Campillo, J. S. Dolado, L. MartínMoreno, J. Bravo-Abad, and F. J. García-Vidal, Opt. Lett. 29, 2500 (2004). [9] F. Medina, J. A. Ruiz-Cruz, F. Mesa, J. M. Rebollar, J. R. Montejo-Garai, and R. Marqués, Appl. Phys. Lett. 95,071102 (2009). [10] J. B. Pendry, L. Martín-Moreno, and F. J. Garcia-Vidal, Science 305,847 (2004). [11] A. Hibbins, B. Evans, and J. Sambles, Science 308,670 (2005). [12] F. J. García de Abajo and J. J. Sáenz, Phys. Rev. Lett. 95,233901 (2005). 245425-6 THEORETICAL AND EXPERIMENTAL EXPLORATION OF . . . PHYSICAL REVIEW B 95, 245425 (2017) [13] C. Cutler, Proceedings of IEEE Antennas and Propagation Society International Symposium and URSI National Radio Science Meeting, Vol. 3 (IEEE, Seattle, WA, 1994), p. 1456. [14] R. Collin, Field Theory of Guided Waves (IEEE Press, Piscataway, New Jersey, 1991). [15] F. J. Garcia-Vidal, L. Martin-Moreno, T. W. Ebbesen, and L. Kuipers, Rev. Mod. Phys. 82,729 (2010). [16] F. J. Garcia De Abajo, Rev. Mod. Phys. 79,1267 (2007). [17] J. Bravo-Abad, A. Degiron, F. Przybilla, C. Genet, F. J. GarcíaVidal, L. Martín-Moreno, and T. W. Ebbesen, Nat. Phys. 2,120 (2006). [18] A. I. Fernández-Domínguez, F. J. García-Vidal, and L. MartínMoreno, Phys.Rev.B76,235430 (2007). [19] C. L. Holloway, E. F. Kuester, J. A. Gordon, J. O’Hara, J. Booth, and D. R. Smith, IEEE Antennas Propag. Mag. 54,10 (2012). [20] J. A. Dockrey, S. A. R. Horsley, I. R. Hooper, J. R. Sambles, and A. P. Hibbins, Sci. Rep. 6,22018 (2016). [21] D. Sievenpiper, L. Zhang, R. F. Jimenez Broas, N. G. Alexöpolous, and E. Yablonovitch, IEEE Trans. Microw. Theory Techn. 47,2059 (1999). [22] J. A. Dockrey, M. J. Lockyear, S. J. Berry, S. A. R. Horsley, J. R. Sambles, and A. P. Hibbins, Phys.Rev.B87,125137 (2013). [23] C. Pfeiffer and A. Grbic, IEEE Trans. Antennas Propag. 58, 3055 (2010). [24] M. Bosiljevac, M. Casaletti, F. Caminita, Z. Sipus, and S. Maci, IEEE Trans. Antennas Propag. 60,4065 (2012). [25] S. Maci, G. Minatti, M. Casaletti, and M. Bosiljevac, IEEE Antennas Wireless Propag. Lett. 10,1499 (2011). [26] N. Yu and F. Capasso, Nat. Mater. 13,139 (2014). [27] N. Yu and F. Capasso, J. Lightwave Technol. 33,2344 (2015). [28] M. Camacho, R. R. Boix, and F. Medina, Phys. Rev. E 93, 063312 (2016). [29] S. R. Rengarajan, IEEE Antennas Wireless Propag. Lett. 4,47 (2005). [30] P. P. Ewald, Ann. Phys. 369,253 (1921). [31] K. E. Jordan, G. R. Richter, and P. Sheng, J. Comput. Phys. 63, 222 (1986). [32] COMSOL AB, COMSOL Multiphysics R. [33] B. A. Munk, Frequency Selective Surfaces: Theory and Design (Wiley, Hoboken, NJ, 2000). [34] A. Neto and S. Maci, IEEE Trans. Antennas Propag. 51,1572 (2003). [35] I. R. Hooper, B. Tremain, J. A. Dockrey, and A. P. Hibbins, Sci. Rep. 4,7495 (2014). 245425-7