Far-Field and Near-Field Physics of Extraordinary THz Transmitting Hole-Array Antennas

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 9, SEPTEMBER 2019 6029 Far-Field and Near-Field Physics of Extraordinary THz Transmitting Hole-Array Antennas Miguel Camacho ,Student Member, IEEE, Rafael R. Boix ,Member, IEEE, Sergei A. Kuznetsov ,Member, IEEE, Miguel Beruete , and Miguel Navarro-Cía ,Senior Member, IEEE Abstract— Despite three decades of effort, predicting accurately extraordinary transmission through subwavelength hole arrays has proven challenging. The lack of quantitative design and modeling capability to take into account the inherent complexity of high frequency instrumentation has prevented the development of practical high-performance components based on this phenomenon. This paper resorts to the Method of Moments to provide not only such missing quantitative prediction but also a theoretical framework to understand and shed more light on the far-field and near-field physics of the extraordinary terahertz (THz) transmission through subwavelength hole arrays under different illumination and detection conditions. An excellent agreement between the numerical and experimental results with various illumination and detection setups is obtained, demonstrating the suitability of this computationally efficient modeling tool to predict the response of extraordinary transmission structures in practical situations. Index Terms— Extraodinary transmission, frequency selective surface, Method of Moments, quasi-optics, terahertz (THz), time-domain spectrometer. Manuscript received October 25, 2018; revised April 30, 2019; accepted May 11, 2019. Date of publication June 19, 2019; date of current version September 4, 2019. The work of M. Camacho was supported by the Engineering and Physical Sciences Research Council (EPSRC), U.K., through the EPSRC Centre for Doctoral Training in Metamaterials under Grant EP/L015331/1. The work of S. A. Kuznetsov was supported by the Russian Foundation for Basic Research (RFBR) under Grant 18-29-20066. The work of M. Beruete was supported by the Ministerio de Ciencia, Innovación y Universidades under Grant RTI2018-094475-B-I00. The work of M. Navarro-Cía was supported in part by the Royal Society under Grant RSG/R1/180040 and in part by the University of Birmingham (Birmingham Fellowship). (Corresponding authors: Miguel Beruete; Miguel Navarro-Cía.) M. Camacho was with the Department of Physics and Astronomy, University of Exeter, Exeter EX4 4PY, U.K. He is now with the Department of Electrical and Systems Engineering, University of Pennsylvania, Philadelphia, PA 19104 USA (e-mail: [email protected]). R. R. Boix is with the Department of Electronics and Electromagnetism, College of Physics, University of Seville, 41012 Seville, Spain (e-mail: [email protected]). S. A. Kuznetsov is with the Rzhanov Institute of Semiconductor Physics SB RAS, Novosibirsk Branch “TDIAM”, Novosibirsk 630090, Russia, and also with the Physics Department, Novosibirsk State University, Novosibirsk 630090, Russia (e-mail: sakuznets[email protected]u.ru). M. Beruete is with the Antennas Group-TERALAB, Universidad Pública de Navarra, 31006 Pamplona, Spain, and also with the Institute of Smart Cities, Universidad Pública de Navarra, 31006 Pamplona, Spain (e-mail: [email protected]). M. Navarro-Cía is with the School of Physics and Astronomy, University of Birmingham, Birmingham B15 2TT, U.K. (e-mail: m.nav[email protected].uk). Color versions of one or more of the figures in this article are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TAP.2019.2920262 I. INTRODUCTION THE discovery of extraordinary optical transmission (EOT) in the late 1990s [1]–[3] stimulated the research on a type of frequency selective surfaces that had been little studied until then. It was found that, thanks to periodicity, subwavelength hole arrays present a very narrow passband at frequencies slightly below the first onset of diffraction. This feature attracted much attention for its promising filtering and sensing applications at optical frequencies [4]. A few years later, this finding was explained in terms of the coupling of the incident wave to surface plasmon polaritons supported by the metal–air interface [5]. This phenomenon was later found also at microwave frequencies [6], where the role of the surface plasmons was taken by leaky waves, supported by the hole array [7]. Due to the deep connection between the extraordinary transmission and periodicity, practical quasi-optical components such as filters and wave plates [8], [9] usually consist of large, but finite, arrays for which the number of holes required for the transmission peak to appear remains as an open question [10]. Recent studies have approached this problem both experimentally [11], [12] and theoretically [13]–[15], assuming plane wave illumination that fails to capture the subtle underlying mechanisms of practical scenarios. For instance, in truncated periodic arrays, plane waves couple to nonradiating and radiating modes through edge diffraction. However, at large angles of incidence, the coupling efficiency is low and, therefore, does not introduce noticeable effects on the far-field pattern [16]. As shown in [17], the electric field distribution of subwavelength hole arrays operating at the EOT frequency changes sharply when they are excited with a localized source, a fact that is not observed in the theoretical works that use a plane wave illumination or in the experimental works done in the far-field range. The coupling between localized sources and leaky waves has been extensively exploited for the design of Bull’s- Eye [18] and metasurface antennas [19]. However, more effort is needed to exploit the leaky-wave coupling mechanism and improve current applications of EOT in sensing, color filtering, etc. [4] or for the design of extraordinarily transmitting antennas [20]. Given that EOT is an enabling technology for applications in the whole electromagnetic spectrum, one misses a numerically efficient framework that can unlock its This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ 6030 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 9, SEPTEMBER 2019 Fig. 1. Schematic of the freestanding subwavelength hole array along with the incident Gaussian beam. Lattice period a=470 μm and hole side s=230 μm. full potential. In this paper, this issue will be addressed first theoretically with the help of a Method of Moments approach for the analysis of large but finite arrays of holes under different illuminations. This comprehensive study, not carried out before, will cover a wide range of experimental scenarios seen in today’s EOT-based applications and will provide the missing basic design, modeling, and optimization tools. Second, a numerical study considering several illumination and detection cases will be presented and the results will be compared with quasi-optical experimental measurements with the aim to provide physically appealing intuition for the design of these devices in terms of classic leaky-wave theory. Such systematic study has been disregarded by literature but has strong implications for terahertz (THz) free-space metrology. II. METHOD OF MOMENTS ANALYSIS Let us consider an array of Nxby Nysquare holes with lattice period aand lateral hole dimension sperforated on a perfectly conducting screen of negligible thickness, as shown in Fig. 1. This is a good approximation at THz frequencies [21] as shown later in Section IV. The impinging wave is a Gaussian beam whose electric field on the plane z =0isassumedtobeoftheformEinc = E0e−[(x−xc)2+(y−yc)2]/w2 0ˆx,where(xc,yc, z =0)represents the coordinates of the center of the array and w0represents the beam-waist, i.e., beam radius [22]. Then, the electric field on the holes satisfies the electric field integral equation given by [15] 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)(1) where Esc t(x,y, z =0)represents the tangential component of unknown electric field distribution on the surface of the array, consisting of a rectangular grid with M=NxNysquare cells, each of them containing a square hole, which occupies a surface denoted by ηiand GM(x,y)is the dyadic Green’s function 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) with g(x,y)=− je−jk0√x2+y2 πk0Z0x2+y2.(3) At normal incidence, the current excited on the metal by the impinging wave in the absence of holes is given by Jas(x,y)=(2Einc/Z0),whereZ0=(μ0/0)1/2is the free space impedance. The integral equation in (1) can be solved approximately by expanding the unknown electric field on the holes in terms of a set of Nbbasis functions, djl(x,y), such that Esc t(x,y, z =0)≈ Nb  l=1 ejldjl(x,y)(x,y)∈ηj(4) where ηjrepresents the jth hole of the array and where ejl are unknown coefficients. Then, by introducing (4) into (1) and projecting the resulting expression on the same set of basis functions (Galerkin’s version of MoM), one can derive a system of linear equations for the unknown amplitudes ejl given by M  j=1 Nb  l=1 kl ijejl =pik (i=1,...,M;k=1,...,Nb)(5) 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,...,Nb)(6) and where pik =−ηi d∗ ik(x,y)Jas(x,y)dxdy (i=1,...,M;k=1,...,Nb). (7) As shown in [15], the matrix coefficients kl ij can be efficiently and accurately calculated in the spatial domain in terms of cross-correlations between basis functions and their divergences. When Chebyshev polynomials weighted by the appropriate edge behavior of the electric field for each of the polarizations are chosen as basis functions (see [23, (30) and (31)] for the definition of these basis functions), these cross-correlations can be calculated in a closed form [15], [24], reducing the number of numerical integrations required to two. Note that these can also be efficiently calculated using Ma–Rohklin–Wandzura quadratures [25]. CAMACHO et al.: FAR-FIELD AND NEAR-FIELD PHYSICS OF EXTRAORDINARY THz TRANSMITTING HOLE-ARRAY ANTENNAS 6031 For the calculation of the coefficients pik in (7), we have made use of the Gauss-Chebyshev quadratures [26]. With this method, we have been able to efficiently simulate the transmission through large arrays of up to 2500 holes, which would require an enormous computational effort using commercial electromagnetic simulators. Once the system of equations is solved, the values of ejl are retrieved, and from them, one can obtain the far-field distribution. This is done through the 2-D continuous Fourier transform, which can be calculated analytically thanks to the choice of basis functions with analytical Fourier transform as  Esc t(kx,ky, z =0)= M  j=1 Nb  l=1 ejl djl(kx,ky). (8) Using this spectral decomposition, one can calculate the electric field at any given position of space as Esc t(x,y, z ) =∞ −∞ ∞ −∞ Esc t(kx,ky, z =0)e−j(kxx+ky+k zz )dkxdky(9) where k z =k2 0−k2 x−k2 y.(10) If the point of observation is in the far-field region, the previous expression can be obtained in an asymptotic form, as shown in [27, Sec. 4.1], given by Esc(r,θ,φ)  z >0 =jk0e−jk0r 2πr × Esc x(kx=k0sin θcos φ,ky=k0sin θsin φ, z =0)cos ϕ + Esc y(kx=k0sin θcos φ,ky=k0sin θsin φ, z =0)sin ϕˆ θ + Esc y(kx=k0sin θcos φ,ky=k0sin θsin φ, z =0)cos ϕ − Esc x(kx=k0sin θcos φ,ky=k0sin θsin φ, z =0)sin ϕ ×cos θˆϕ.(11) In that case, the transmission coefficient can then be computed numerically as T=2 πE0w2 0π/2 02π 0|Esc(r,θ,φ)|2r2sin θdφdθ. (12) As discussed later in detail, in most quasi-optical experimental setups, it is not possible to locate the receiving antenna in the far-field range, and one needs to resort to the numerical evaluation of (9). III. NUMERICAL RESULTS Let us first explore the dependence of the EOT phenomena on both the number of holes and size of the Gaussian beam (represented by its beam-waist w0). In Fig. 2(a)–(c), we have plotted the transmission coefficient [as calculated in (12)] for increasing beam-waists of 1, 2, and 10 times the periodicity, respectively. For each of these cases, we study the evolution of the array response with the number of elements on the x-direction, as the direction set by the impinging polarization plays the main role in the appearance of the EOT phenomena for highly symmetric unit cells. This is due to the transverse-magnetic nature of the surface waves supported by a freestanding metal-connected structure as the one studied here, which is incompatible with the transverse-electric nature of the diffracted orders that propagate on the surface of the array along the y-direction at their onset [7], [28]. This mismatch provokes the disappearance of Wood’s anomaly, which is normally connected to the constructive interference of the diffracted orders at their onset. For this reason, the number of elements along the y-direction does not play an important role [29] and is kept constant but large enough for the array to cover the widest beam spot. In Fig. 2, all three cases show very well-defined EOT peaks with as few as three holes, although the number of apertures required to obtain saturation of the transmission peak varies depending on the size of the spot. In particular, the spot sizes of a,2a, and 10arequire a minimum of 3, 7, and 21 elements, respectively, confirming that the number of nonilluminated holes play a fundamental role in the EOT resonance due to the excitation of a leaky wave that runs along the surface, as outlined in [12]. The leaky modes (the m=−1 space harmonic of the surface wave, supported by the hole array propagating away from the illumination spot) allow for the exploration of nonilluminated holes, which help reproduce the behavior of a large periodic array even with a confined illumination. In this process, the effect introduced by the edges is greatly diminished due to the exponential decay of these leaky modes due to the gradual radiation of energy. Otherwise, two wave mechanisms will happen at the edge [16], [30]: (1) diffraction and (2) excitation of back-propagating modes that can be either nonradiating like surface waves (e.g., creeping waves and Norton waves [14], [31]), and all the other higher order Floquet wave modes (i.e., space harmonics), or radiating like the back-propagating leaky wave m=−1 mode. Only the latter could significantly contribute to radiation at broadside. Although the excitation of these back-propagating m=−1 leaky modes could be accomplished through the continuous spectrum of the diffracted fields on the edges, the efficiency of the process is low unless the edge is engineered as in [32] and [33]. Fig. 2 shows that modifying the spectrum of the illumination is a far more efficient approach when the edge is not engineered. In addition, on the nonilluminated half-space one could consider the hole array as an antenna, whose size will depend on the size of the illumination. As we will discuss later, changing the size of the beam spot will have a large effect on the distance at which the antenna operates in far-field conditions. This miniaturization of the EOT array, however, cannot be done at zero cost. As shown by Fig. 2, reducing the beam waist diminishes the transmissivity from the total transmission found in periodic arrays under normal plane wave incidence due to the wider spectrum of the Gaussian beams. We find that the values of the maximum transmission correspond to −2, −1.7, and −0.2 dB for beam-waists of a,2a,and 10a, respectively. 6032 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 9, SEPTEMBER 2019 Fig. 2. Transmission coefficient for varying number of square holes Nx under Gaussian beam illumination with beam-waist (a) w0=a,(b)w0=2a, and (c) w0=10aand for Ny=21. The dashed line marks the onset of the first Wood’s anomaly for the nontruncated 2-D periodic array at normal plane wave illumination. The presence of leaky waves has a strong effect on the radiating properties of the array, as explained in the following. Fig. 3(a)–(c) represents the evolution of the radiation diagram on the E-plane (x− z in our particular case) when the number of elements along the x-direction is increased analogously to Fig. 2. To discern better the effect of the presence of the leaky waves, let us first focus on Fig. 3(a) where the beam-waist corresponds to the periodicity. In this case, in which a single unit cell is being illuminated, we can see a great reduction in the power radiated into nonbroadside directions with the introduction of two lateral holes. One can see that, when the number of apertures increases, the maximum of the radiation is not at broadside, which is consistent with the Fig. 3. Radiation pattern (E-plane) of the extraordinary THz transmitting antenna for varying number of square holes Nxunder Gaussian beam primary feeding with beam-waist (a) w0=a,(b)w0=2a,and(c)w0=10aand for Ny=21. These have been obtained at the frequencies of the maximum transmission shown in Fig. 2. leaky mode (the m=−1 space harmonic of the complex wave supported by the hole array) presenting a nonzero wavevector component along the x-direction, as we will later confirm. In fact, strictly speaking, the absolute broadside radiation is not possible because of the open stopband [34]. Additionally, a noticeable amount of energy is radiated into large angles. This can be associated with the fundamental m=0 space harmonic of the complex wave, which runs very close to the light line [35] (i.e., it is almost like a plane wave at grazing angle). As the number of apertures is increased, these two contributions become clearer, because the contribution from the diffracted field by the edges of the array is reduced. CAMACHO et al.: FAR-FIELD AND NEAR-FIELD PHYSICS OF EXTRAORDINARY THz TRANSMITTING HOLE-ARRAY ANTENNAS 6033 Fig. 4. Normalized surface scattered electric field distribution along the E-plane of the truncated subwavelength hole array for varying beam-waist. Green: w0=2a.Red:w0=5a.Blue:w0=10a. Horizontal dashed line: normalized electric field amplitude of the impinging Gaussian beam on the central hole. Comparing Fig. 3(a)–(c), the main differences are found in the convergence of the energy distribution. In Fig. 3(b), for instance, one observes a characteristic pattern in the region of 30◦–50◦that arises from the interference among the increasing number of elements illuminated by the beam spot, a feature that was not present in Fig. 3(a). However, the regions of small and large angles are dominated by the leaky waves (m=−1andm=0 space harmonic, respectively) and, therefore, are less affected by the change in the number of illuminated holes. Finally, in Fig. 3(c), in which the largest beam spot is considered, the radiation pattern is dominated by the interference between the holes directly illuminated, with little energy going into the excitation of leaky modes, due to the narrow spectrum of the beam. Most of our experimental studies, until now devoted to EOT [8], [20], [36]–[39] and those reported by others at optics and infrared (see [5] and references therein), have been dealing with this situation. The radiation patterns shown in Fig. 3 allow us to study the spectral distribution of the electric field on the surface of the hole array, leading to different species of radiating contributions. The presence of these species can be easily noticeable from the electrical field distribution itself, as shown in Fig. 4, where the electric field at the center of the hole has been represented for the central row of an array with Nx=107 and for beam-waists equivalent to 2a,5a, and 10a, respectively, at the frequency of the EOT peak. The curves have been normalized to the value of w0=10a. Additionally, the dashed lines represent the electrical field amplitude of the Gaussian beam on the central hole, normalized again to the maximum of w0=10a. Within a local periodicity assumption, one would expect the maximum amplitude of the electrical field distribution to grow proportionally to that of the local driving field. However, this approximation does not take into account that changing the beam spot strongly modifies the spectral distribution of energy on the surface, as we have seen in the radiation patterns. This means that, as the beam spot becomes narrower, more energy is being captured by the leaky modes, and there is less energy available for each Fig. 5. Complex kxspace showing the complex zeros of the determinant of the matrix kl per for an infinitely periodic array at (a) 0.58 THz and (b) 0.76 THz. The vertical dashed lines at (a) kxa/π =±1.821 and (b) kxa/π =±2.4 represent the normalized value of the free space wavevector k0and thus, the limits of the radiation region for each frequency. hole to resonate with the locally impressed field. Therefore, we find that the field on the central hole is reduced by a 10% and 27% with respect to what one would obtain through a local periodicity assumption when the array is illuminated with beam-waists of 5aand 2a, respectively. To confirm the existence of these leaky waves, we consider an infinite array along both xand y-directions and calculate the complex propagation constants of the modes traveling along x in the absence of excitation. More specifically, the procedure consists in finding the complex zeros of the determinant of the system of equations that has the same set of basis functions and the 2-D periodic Green’s function, following the Method of Moments approach presented in [40], at a given frequency. In Fig. 5, we show the colormap of the value of the determinant of the matrix kl per (representing the coefficient matrix of the system of linear equations) in logarithmic scale at two different frequencies. Fig. 5(a) and (b) have been calculated at the EOT transmission peak (0.58 THz) and at a frequency above the first Wood’s anomaly (0.76 THz), respectively. Let us focus on Fig. 5(a). In that, one can clearly see the presence of 8 zeroes, half of them with positive imaginary part of their wavevector and half of them with negative imaginary part. Thanks to the symmetry of the unit cell, the position of these zeroes is symmetric. To understand their origin, one needs to remember that, due to the periodicity of the infinite array, the reciprocal space of the xvariable can be obtained by periodically repeating the Brillouin zone defined by the region limited by kx∈[−π/a,π/a]. In addition, hole arrays are known to support complex waves, which, for the case of holes that are small compared to the wavelength and to the periodicity, present a value of their wavevector that is just slightly larger than the free space wavevector, |ksw x|≈ k0[7], [40]. Hence, when k0approaches the Brillouin zone boundary, thanks to the periodicity of the spectrum, some of the harmonics with wavevector km x=ksw x+2πm/amay enter the visible region [41] and become a leaky mode. According to this reasoning, we can identify the zeroes with wavevector 6034 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 9, SEPTEMBER 2019 close to k0as the zeroth harmonics of the complex wave supported by the hole array. Then, as the array is being excited at a frequency below the first onset of diffraction, the other zeroes correspond to the m=−1 space harmonics of those complex waves that propagate away from the illumination spot. This zero appears in the negative Re(kx)-space. The symmetric zeroes linked to the m=1 also appear in the Brillouin zone, but they are only physically meaningful when the array supports leaky waves propagating back to the illumination spot (e.g., in truncated arrays) [41]. The bright spots correspond to the wavevector of the diffracted orders, which lead to poles in the determinant of the system of equations [42]. The exact positions of these zeroes are kx,−1=±k0(0.1±j0.029)and kx,0=±k0(1.0±j0.029). In Fig. 5(b), one can see the presence of 12 zeroes. Following the aforementioned rationale, the zeroth harmonics can be identified as those zeroes with a real part of kxslightly larger than k0, as the phase velocity of these waves reduces when the frequency approaches the λ/2 resonance of a single hole [40]. Then, using the formula for the space harmonics and what we learned from Fig. 5(a), the four zeroes with the smallest value of the real part of kxcan be identified as the m=−1 space harmonics, which, above the first Wood’s anomaly, shift away from the origin of the complex kxplane as frequency increases [35]. Furthermore, the value of k0surpasses the limits of the first Brillouin zone, and thus, it does not contribute to the far-field radiation. However, m=−2 space harmonics enter the radiation region [35] and contribute to the far-field. The exact positions of the 12 zeroes are kx,−1=±k0(0.31 ±0.028),kx,−2=±k0(0.53 ±0.028), and kx,0=±k0(1.15 ±0.028). Finally, from each of the aforementioned pairs of leaky modes of the hole array, only those that decay in amplitude along the direction on which their energy propagates will be physical solutions. In particular, only the zeroes of the determinant with opposite signs of real and imaginary parts are physical solutions. This understanding of the resonant modes of a nontruncated array of holes confirms our hypothesis for the radiation diagrams shown in Fig. 3, as we can associate the maximum radiation in a direction close (but not equal) to broadside with the energy carried by the m=−1 space harmonics of the complex wave supported by the array and the large energy radiated into large angles as the energy carried by the m=0 harmonic of the aforementioned complex wave. In addition, the radiation diagrams present the contribution from the nontruncated array at directions close to broadside, as well as diffraction from the edges of the array. IV. MEASUREMENTS AND DISCUSSION To corroborate the previous numerical study, we performed several experimental tests in the THz range. Quasi-optical systems are the preferred setups for measurements at THz frequencies [22]. However, they rarely operate in the far-field range, and thus, comparison between measurements and theoretical or numerical results should be done with caution. Given the leaky-wave nature of the extraordinary transmission resonance, one should expect a strong influence of Fig. 6. Measured transmission coefficient as a function of angle of emission (E-plane) and frequency for a truncated (a) 7 ×7 and (b) 107 ×107 subwavelength hole array under collimated beam illumination and detection (setup 1). (c) On-axis experimental and modeled transmission coefficient for a truncated (a) 7×7 and (b) 107×107 subwavelength hole array. Inset: schema of the experimental setup. the quasi-optical system in the angle-resolved emission pattern from the subwavelength hole array. The importance of the illumination has been recently reported for the on-axis detection [12]. Here, we focus on the angle-resolved emission pattern as a function of the incident beam-waist in the light of the Method of Moments’ findings. To do a fair comparison between the numerical and experimental results in this section, (9) is numerically evaluated at a distance congruent with the measurements. Two different substrate-free truncated subwavelength hole arrays (7 ×7 and 107 ×107) are investigated with three CAMACHO et al.: FAR-FIELD AND NEAR-FIELD PHYSICS OF EXTRAORDINARY THz TRANSMITTING HOLE-ARRAY ANTENNAS 6035 Fig. 7. Measured transmission coefficient as a function of angle of emission (E-plane) and frequency for a truncated (a) 7 ×7 and (b) 107 ×107 subwavelength hole array under focused beam illumination and detection (setup 2). (c) On-axis transmission coefficient for a truncated (a) 7 ×7and (b) 107 ×107 subwavelength hole array. Inset: schematic of the experimental setup. different quasi-optical systems whose schema can be found in Figs. 6–8 as insets. The arrays were fabricated with the electroplating technology described in [43]. Technical details of the experimental settings and the quasi-optical systems used in this work can be found in the Appendix. From Section III and the observation in Bull’s Eye antennas supporting EOT [35], one should expect the m=−2 diffraction order to emerge in angle-frequency maps with a frequency-dependent angular response. Color maps of Figs. 6–8 confirm the existence of the m=−2 space harmonic of the complex wave supported by the array for all cases above Fig. 8. Measured transmission coefficient as a function of angle of emission (E-plane) and frequency for a truncated (a) 7 ×7 and (b) 107 ×107 subwavelength hole array under focused beam illumination and no collimating lens in detection (setup 3). (c) On-axis experimental and modeled transmission coefficient for a truncated (a) 7 ×7 and (b) 107 ×107 subwavelength hole array. Inset: schematic of the experimental setup. the EOT frequency and at large radiation angles as previously predicted by the Method of Moments in Fig. 5(b). From previous observations in Bull’s Eye antennas [35], one should expect dispersion in the m=−1 diffracted order responsible for the EOT peak. Such characteristic frequency-dependent angular response is only noticeable unambiguously in setup 3 for 107 ×107 [see Fig. 8(b)], which is, arguably, the closer analog among the three setups to the Bull’s Eye antenna. This case shows around the EOT frequency a clear avoided crossing at broadside of the m=−1 space harmonics [41] each supported by each half of the subwavelength hole array, 6036 IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. 67, NO. 9, SEPTEMBER 2019 Fig. 9. Comparison between measured and numerically computed angle-resolved transmission coefficient (E-plane) at the resonance frequency for truncated (7×7 and 107×107) subwavelength hole arrays. (a) Collimated beam setup (w0=10a). (b) Focused beam setup (w0=3a). (c) Focused beam illumination (w0=3a)+collimated beam detection setup. while in the 7 ×7 counterpart [Fig. 8(a)] such expected avoided crossing is masked by the Fresnel diffraction. This avoiding of crossing produces a mode splitting at normal incidence and leads to the existence of a maximum transmission below the first onset of diffraction. This phenomenon was more clearly observed in plane wave simulations for single- (see [44, Fig. 4]) and multi-layered structures, i.e., fishnet metamaterial (see [45, Fig. 3]). A close look at setup 2 results [Fig. 7(a) and (b)] may also reveal the avoided crossing, but such response disappears completely in setup 1 results [Fig. 6(a) and (b)] because of the nonlocalized nature of the collimated illumination, as the radiation pattern is dominated by the interference between the holes directly illuminated. For validation purposes, we have plotted the normalized numerical on-axis transmission coefficient in Figs. 6(c)–8(c) for frequencies around the EOT peak, at which the beam-waist was experimentally determined to be approximately 10aand 3a, respectively. We observe an excellent agreement between the numerical and measured transmission coefficient in the three experimental configurations. Contrary to the experiment, only 21 rows parallel to the impinging electric field were included in the simulations, as we know from the results in Fig. 5 and from [46] and [42], that the leaky waves responsible for the EOT propagate along such direction. Such approximation was verified for the cases shown in Figs. 6(c) and 8(c). This allowed us to reduce the computational burden of the problem, leading to a code that calculates any of the cases in Figs. 6(c) and 8(c) in less than 1 h in a personal computer with one processor Intel Core i7 @ 3.6 GHz with 32 GB of RAM, whereas the same result requires 76 h using a full-wave commercial simulator running on a HP Z820 with two processors Intel Xeon E5-2660 @ 2.2 GHz with 128 GB of RAM [12]. Fig. 9 shows the measured and numerically calculated angle-resolved data at the extraordinary transmission peak for the six different scenarios. One can see a good agreement between the measured and numerical results. Transmission is not maximum on-axis but slightly off-axis (1◦)forsetup1 [Fig. 9(a)] and setup 3 [Fig. 9(c)]. This feature is well captured by the Method of Moments. This response agrees with the leaky-wave formalism and the fact that an open stopband at broadside should be expected for symmetric EOT structures as those reported here. This key feature was not observed in the past because of the experimental limitations [46]. Remarkably, at such angle of radiation, the measured transmission is indeed above 0 dB for the 107 ×107 sample in both setups. This gain stems from the larger effective (antenna) area that the leaky wave produces at the exit interface of the subwavelength hole array compared to the beam spot illuminating the subwavelength hole array and the negligible absorption of the high-quality freestanding samples fabricated by electroplating technology. V. CONCLUSION In this paper, the physics of extraordinary THz transmission plate antennas has been explored. Thanks to an efficient implementation of the Method of Moments for the analysis of the transmission through large square arrays of holes under Gaussian beam illumination, we have elucidated the main mechanisms that control both the excitation of surface modes and their subsequent radiation properties. We have found that the size of the beam spot controls the amplitude of the variety of leaky modes supported by the array and, therefore, allows us to manipulate the far-field energy distribution on the other side of the plate. This intuitive mechanism has been corroborated by studying the electrical field distribution both on the surface of the plate and in the Fresnel and far-field regions. The last two have also been experimentally explored in detail both to validate our theoretical analysis and to provide CAMACHO et al.: FAR-FIELD AND NEAR-FIELD PHYSICS OF EXTRAORDINARY THz TRANSMITTING HOLE-ARRAY ANTENNAS 6037 guidance on the practical limitations presented by the traditional quasi-optical measurements at THz frequencies. These limitations are related to the large differences in the position of the Fresnel-zone boundary introduced by a tunable beam size, due to the change in the radiating area of the antenna. These different regimes can, if not dealt with properly, introduce large discrepancies between theoretical predictions and experimental measurements, which are analyzed here with the help of numerical near-field calculations. APPENDIX The samples were characterized with the all fibercoupled THz time-domain spectrometer TERA K15 from Menlo Systems in a quasi-optical configuration without purging. The lock-in constant was set to 300 ms and the total temporal length of the recorded waveforms was at least 208 ps to have a spectral resolution of 4.8 GHz in the worst case. Three different sets of optics were used to realize three different configurations [see insets in (c) of Figs. 6–8]. In all cases, the detector unit was placed on a free rotation arm, whose pivot point coincides with the sample position. The distance between the detector unit and the sample was approximately 110 mm. Setup 1 dealt with collimated beam illumination. At the sample position, the frequency-dependent beam-waist of the THz beam was estimated to be 5 mm at 0.6 THz. Setup 2 and setup 3 used focused beam illumination. The beam-waist was estimated to be 1 mm at 0.6 THz in these cases. Unlike setup 3, setup 2 has a TPX planoconvex lens (effective local length 54 mm) on the free rotation arm of the quasi-optics to collimate radiation emerging from the samples. This collimating lens effectively increases the numerical aperture of the detection. An angular step of 1◦was used until 20◦and then an angular step of 2.5◦was applied in the measurements. All configurations were measured twice in different days to check that the results were repeatable. Calibration was done by comparing the measurements with the aligned configuration (i.e., emitter and detector on-axis) without the sample on the sample holder. REFERENCES [1] E. Betzig, A. Lewis, A. Harootunian, M. Isaacson, and E. Kratschmer, “Near field scanning optical microscopy (NSOM): Development and biophysical applications,” Biophys. J., vol. 49, no. 1, pp. 269–279, 1986. [2] R. E. Betzig, “Nondestructive optical imaging of surfaces with 500 angstrom resolution,” Ph.D. dissertation, School Appl. Eng. Phys., Cornell Univ., Ithaca, CA, USA, 1988. [3] T. W. Ebbesen, H. J. Lezec, H. F. Ghaemi, T. Thio, and P. A. Wolf, “Extraordinary optical transmission through sub-wavelength hole arrays,” Nature, vol. 391, pp. 667–669, Feb. 1998. [4] S. G. Rodrigo, F. de Leon-Perez, and L. Martin-Moreno, “Extraordinary optical transmission: Fundamentals and applications,” Proc. IEEE, vol. 104, no. 12, pp. 2288–2306, Dec. 2016. [5] F. J. Garcia-Vidal, L. Martin-Moreno, T. W. Ebbesen, and L. Kuipers, “Light passing through subwavelength apertures,” Rev. Mod. Phys., vol. 82, pp. 729–787, Apr. 2010. [6] M. Beruete et al., “Enhanced millimeter wave transmission through quasioptical subwavelength perforated plates,” IEEE Trans. Antennas Propag., vol. 53, no. 6, pp. 1897–1903, Jun. 2005. [7] V. Lomakin and E. Michielssen, “Transmission of transient plane waves through perfect electrically conducting plates perforated by periodic arrays of subwavelength holes,” IEEE Trans. Antennas Propag., vol. 54, no. 3, pp. 970–984, Mar. 2006. [8] M. Beruete, M. Navarro-Cía, I. Campillo, P. Goy, and M. Sorolla, “Quasioptical polarizer based on self-complementary sub-wavelength hole arrays,” IEEE Microw. Wireless Compon. Lett., vol. 17, no. 12, pp. 834–836, Dec. 2007. [9] V. Torres, N. Sánchez, D. Etayo, R. Ortuño, M. Navarro- Cía, A. Martínez, and M. Beruete, “Compact dual-band terahertz quarter-wave plate metasurface,” IEEE Photon. Technol. Lett., vol. 26, no. 16, pp. 1679–1682, Aug. 15, 2014. [10] A. Bitzer and M. Walther, “Terahertz near-field imaging of metallic subwavelength holes and hole arrays,” Appl. Phys. Lett., vol. 92, no. 23, p. 231101, Jun. 2008. [11] F. Przybilla et al., “Efficiency and finite size effects in enhanced transmission through subwavelength apertures,” Opt. Express, vol. 16, no. 13, p. 9571, Jun. 2008. [12] M. Navarro-Cía, V. Pacheco-Peña, S. A. Kuznetsov, and M. Beruete, “Extraordinary THz transmission with a small beam spot: The leaky wave mechanism,” Adv. Opt. Mater., vol. 6, no. 8, Apr. 2018, Art. no. 1701312. [13] J. Bravo-Abad, A. Degiron, F. Przybilla, C. Genet, F. J. García-Vidal, L. Martín-Moreno, and T. W. Ebbesen, “How light emerges from an illuminated array of subwavelength holes,” Nature Phys., vol. 2, no. 2, pp. 120–123, Feb. 2006. [14] H. Liu and P. Lalanne, “Microscopic theory of the extraordinary optical transmission,” Nature, vol. 452, no. 7188, pp. 728–731, 2008. [15] M. Camacho, R. R. Boix, and F. Medina, “Computationally efficient analysis of extraordinary optical transmission through infinite and truncated subwavelength hole arrays,” Phys. Rev. E, Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top., vol. 93, no. 6, p. 063312, Jun. 2016. [16] M. Camacho, A. P. Hibbins, F. Capolino, and M. Albani, “Diffraction by a truncated planar array of dipoles: A Wiener–Hopf approach,” Wave Motion, vol. 89, pp. 28–42, Jun. 2019. [17] M. Camacho, R. R. Boix, and F. Medina, “Comparative study between resonant transmission and extraordinary transmission in truncated periodic arrays of slots,” in Proc. IEEE MTT-S Int. Conf. Numer. Electromagn. Multiphysics Model. Optim. RF, Microw., Terahertz Appl. (NEMO), May 2017, pp. 257–259. [18] M. Beruete, U. Beaskoetxea, and T. Akalin, “Flat corrugated and Bull’seye antennas,” in Aperture Antennas for Millimeter and Sub-Millimeter Wave Applications. Signals and Communication Technology.Cham, Switzerland: Springer, 2018, pp. 111–141. [19] G. Minatti et al., “Metasurface antennas,” in Aperture Antennas for Millimeter and Sub-Millimeter Wave Applications. Cham, Switzerland: Springer, 2018, pp. 289–333. [20] M. Beruete, I. Campillo, J. E. Rodríguez-Seco, E. Perea, M. Navarro-Cía, I. J. Núñez-Manrique, and M. Sorolla, “Enhanced gain by double-periodic stacked subwavelength hole array,” IEEE Microw. Wireless Compon. Lett., vol. 17, no. 12, pp. 831–833, Dec. 2007. [21] M. Naftaly, Terahertz Metrology. Boston, MA, USA: Artech House, 2015. [22] P. F. Goldsmith, Quasioptical Systems: Gaussian Beam Quasioptical Propogation and Applications. Piscataway, NJ, USA: IEEE Press, 1998. [23] R. Florencio, R. R. Boix, and J. A. Encinar, “Fast and accurate MoM analysis of periodic arrays of multilayered stacked rectangular patches with application to the design of reflectarray antennas,” IEEE Trans. Antennas Propag., vol. 63, no. 6, pp. 2558–2571, Jun. 2015. [24] R. Florencio, R. R. Boix, and J. A. Encinar, “Enhanced MoM analysis of the scattering by periodic strip gratings in multilayered substrates,” IEEE Trans. Antennas Propag., vol. 61, no. 10, pp. 5088–5099, Oct. 2013. [25] J. Ma, V. Rokhlin, and S. Wandzura, “Generalized Gaussian quadrature rules for systems of arbitrary functions,” SIAM J. Numer. Anal., vol. 33, no. 3, pp. 971–996, Jun. 1996. [26] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions. New York, NY, USA: Dover, 1970. [27] R. E. Collin, Antennas and Radiowave Propagation.NewYork,NY, USA: McGraw-Hill, 1985. [28] M. Beruete, M. Navarro-Cía, and M. S. Ayza, “Understanding anomalous extraordinary transmission from equivalent circuit and grounded slab concepts,” IEEE Trans. Microw. Theory Techn., vol. 59, no. 9, pp. 2180–2188, Sep. 2011. [29] J. Bravo-Abad, F. J. García-Vidal, and L. Martín-Moreno, “Resonant transmission of light through finite chains of subwavelength holes in a metallic film,” Phys.Rev.Lett., vol. 93, no. 22, Nov. 2004. [30] G. Minatti, E. Martini, and S. Maci, “Efficiency of metasurface antennas,” IEEE Trans. Antennas Propag., vol. 65, no. 4, pp. 1532–1541, Apr. 2017.