Geometrical Dependence on the Onset of Surface Plasmon Polaritons in THz Grid Metasurfaces

The transmission response of metallo-dielectric grid metasurfaces is experimentally investigated through Terahertz Time Domain Spectroscopy and the corresponding effective dielectric function is retrieved. Using a lumped element model we can determine the dependence of the effective plasma frequency (the transition frequency) on the metasurface filling factor F. The change of the transition frequency vs. F spans over one order of magnitude and sets the threshold between the metamaterial (homogeneous) and the photonic crystal (diffraction-like) regime, ruling the onset of two different Surface Plasmon Polaritons, spoof and high order. Field symmetry and spatial extension of such excitations are investigated for the possible applications of THz grid metasurfaces in bio- and chemical sensing and sub-wavelength imaging.

Surface Plasmon Polaritons (SPPs) are collective excitations induced at the boundary between metallic and dielectric materials. The underlying physics attracted a huge interest in the last few years because of their high potential in real world applications 1 . The possibility to control at will the activation energy, the extension, and the propagation length of SPPs make them optimal candidates for sub-wavelength imaging 1,2 , lithography 3 , surface sensing 4 , wave-guiding 5 , local field emission 6,7 , and many other applications.
Furthermore, the recent development of THz technology, based on reliable, compact, and high-performance sources and detectors 8 , is giving a boost to the study of tailored 9 "cold" (in terms of energy) SPPs for the development of nearly zero refractive index metamaterials 10 and for the control of enhanced transmission phenomena [11][12][13][14][15][16][17] .
The surface mode phenomenology may be divided into two categories: the spoof and the tunneling SPPs. Spoof Surface Plasmon Polaritons (SSPPs) refer to modes that are generated and remain confined at the metal-dielectric interface (MDI) 18 , with an appropriate electromagnetic perturbation having a frequency close to the plasma frequency ω p . The latter terminology pertains instead to high order modes (HOSPPs) excited at MDI but that propagate across the sample thickness, favoring the photon tunneling 15 .
In ordinary metals the plasma frequency can be as high as 10 15 Hz 19 and is given by ω ε = ⁎ n e m / p 2 0 , where n is the charge carrier density, e is the elementary charge, ε 0 the vacuum permittivity and m* the effective carrier mass. Spoof SPPs are upper bounded by the plasma frequency and their dispersion has the same signature of a standard SPP 20 . To excite them at lower frequencies, it is possible to decrease ω p diluting the metal layer, with an operation of dimensional "squeezing" in which the unit cell volume Vf of a plane film is reduced to a value V u < V f 20 . In this way, the diluted metal has an effective charge carrier density n′ = nF, where F = V u /V f is the filling factor. Since the effective mass becomes related to F, ω p of a diluted metal relies on the geometrical parameters of the unit cell only 21 .
We study the THz response of metasurfaces intended as diluted metals obtained through an appropriate combination of square symmetry conducting patches of side d, deposited over a dielectric substrate. Their periodical distribution over the surface composes a grid with periodicity p sized the unit cell (grid metasurfaces, GMs). In particular, we investigate the dielectric function of two GMs consisting of a patterned layer of copper, 30 μm thick, deposited over a FR4 foil with a nominal thickness s = 160 μm. Samples are realized through a standard process readily available from PCB industry. The conductive patches have the same periodicity p = 600 μm and the same size d ≈ 0.5p, composing a standard grid (GR, Fig. 1(a)) or a chessboard (CB, Fig. 1(b)) geometry. The latter grid ScIeNtIFIc REPORTS | (2019) 9:924 | DOI:10.1038/s41598-018-36648-x medium can be obtained from the former one by removing the central metallic patch in the intersection of each perpendicular and parallel wire ( Fig. 1(c)), resulting in a conducting network of wires.
To distinguish a GM from a standard MDI realized with an unpatterned metal layer, we define the reduced plasma frequency of the system as a transition frequency ω r , depending on the substrate permittivity and F only [22][23][24] .
We clearly show how ω r provides the threshold for the onset of the high order SPPs 25 , typically observed in structured surfaces characterized by sub-wavelength holes 15 .
Time Domain Spectroscopy (TDS) is used to characterize in the transmission mode the response of the GMs and the onset of surface plasmon polaritons in the 0.1-1.0 THz range. Above ω r , tunnelling produced by the excitation of HOSPPs show up in the spectrum of both GR and CB samples as maxima in the amplitude. Experimental values of the transition frequency well agree with full wave electromagnetic simulations. Using a simple lumped element circuital model to describe the grid metasurfaces, we can obtain the analytical dependence of the transition frequency on F.
The analysis of the frequency spectrum allows to get a clear picture not only on the GM excited states, but also on the concept of dilution, somehow setting a threshold (ruled by ω r ) between the metamaterial/homogeneous and the (metallic) photonic crystal regime. In this respect, we show that homogeneity can be seen as a frequency (energy) related condition rather than being determined from the trivial comparison between unit cell periodicity and the wavelength of the impinging signal. Furthermore, full wave electromagnetic simulations of the transmission response allow to identify the characteristics of the E-field at the HOSPP resonances, providing interesting information on the spatial symmetry and poles of the excitations in both grid metasurfaces.

Results
In Fig. 2(a), a detail of the acquired time dependent electric field (in arbitrary units) of the beam passing through the two metasurfaces under study is shown and compared with the free space. Black curve is the reference (air) signal E ref (t), whereas the blue and red curves refer to the transmission across GR, E GR (t), and CB, E CB (t), respectively. For the sake of clarity, E GR (t) and E CB (t) are shown scaled by the arbitrary factors reported in the figure. We assume that the THz signal is propagating along the z-direction, whereas the grid medium extends over the x-y plane.
The observed time shift in the peak maxima can be used to roughly evaluate the effective sample refractive index. Indeed, we estimate on average a delay of about δt ~ 0.5 ps (only marginally dependent on the geometry) which corresponds to a mean refractive index ≅ + . , lower than the corresponding value for FR4 (n FR4 ~ 2.1). This phenomenon can be accounted for by considering the dielectric function of the system, that in standard plasmonic structures 3 is given by ε ε ε ε ε ). ε  FR4 and ε  m are the dielectric functions of the substrate and the metal layer respectively, the latter one being described through the standard Drude function ε ω ω ωω where ω r is the relaxation frequency. Setting for a lossless diluted metal ω τ = 0 and ω p = 1 THz, we can easily verify that the effective refractive index of the composite metal-FR4 slab is lower than the bare FR4, in qualitative agreement with the experimental results.
Using a standard transmission function theory 8 , one can easily calculate the dielectric function ε  of each sample from the Fourier transform of the time signal. The grid metasurfaces do not present symmetry breaking in the plane of unit cells, therefore we use a retrieval process adopted in previous papers on homogeneous samples 26,27 . Since this procedure cannot account for the presence of diffracting/resonating processes above ω r , the dielectric function ε  GM is obtained at its zero order. Results as a function of frequency are shown in terms of real ε GM ′ ( Fig. 2(b)) and imaginary ε GM ″ ( Fig. 2(c)) components. Red and blue curves in each graph refer to the behavior of ε  GM in CB and GR samples respectively. The dielectric function of the bare FR4 substrate is also reported (black curve). All curves are affected by an uncertainty lower than about 7%, mostly related to the indetermination in the substrate thickness value. Both the real and imaginary part of the dielectric response can be described dividing the frequency spectrum into three regions marked by the transition from the plasmonic regime and the onset of relevant diffraction process approximately at Full plasmonic behavior (ε GM ′ < 0) can be clearly observed in Fig. 2 and ω r , CB /2π ≈ 0.21 THz for GR and CB samples respectively.
Above the transition frequency the impinging radiation couples with higher order modes connected to the Bragg wave vectors (see eq. (12) in Methods), inducing HOSPP resonances in transmission that reflect in a dielectric function with dumped oscillations tending towards ε  FR4 . For f > f D , diffraction starts smearing collective phenomena and further geometrical resonances are barely observed.
Following previous studies 28,29 , the resonating behavior of ε ω  ( ) GM can be in principle modeled through an electrical circuit whose impedance can fit the entire electrodynamics of the grid metasurfaces. However, an in-depth discussion on this is beyond the scope of our paper, focusing instead on the transition from SSPPs to HOSPPs.
The onset of SSPPs is testified by the dispersion diagram ω(k) of the fundamental mode 20 . This is reported in Fig. 2(d) presenting the standard monotonic growth upper bounded by the transition frequency, which scales according to the filling factor of the GM.
To evaluate the dependence of the transition frequency ω r on F, we performed full wave electromagnetic simulations using commercial software. The filling factor can be related to the ratio η = d/p, yielding F GR = 2η − η 2 and F CB = 2η 2 − (2η − 1) 2 for the GR and CB structure respectively (see Methods). Simulation results are shown in Fig. 3, as dot-dashed blue and red curves for the GR and the CB grid metasurfaces respectively. In the case of the CB sample the grid structure is lost for F CB < 0.5, making the concept of plasma frequency meaningless. The function ω r (F) increases almost linearly as long as F < 0.8, to become exponential-like at higher values, so that its variation spans a frequency decade (from 0.1 to 1.1 THz). Analytically, the ω r dependence on F can be easily obtained using an approach based on lumped circuital elements 24 . The model calculates the grid medium dielectric function describing each unit cell as a RLC circuit, with parameters r, l and c representing the resistance, inductance and capacitance per unit length respectively. Specifically, while a simple wire medium can be described with a unit cell composed of a resistance in series with an inductance 23 , a grid metasurface can be modeled, at the lowest order, adding a capacitor in parallel to the basic RL unit cell ( Fig. 1(d)). As shown in 23 (and described in detail in Methods), the real part of the dielectric function for a generic grid medium can be expressed as where ω′ p 2 is the effective plasma frequency, ε ∞ is the asymptotic permittivity, ε d is the substrate dielectric constant and τ = l/r accounts for the relaxation time. The lumped element model defines the equivalence ω′ = l p t 1/ p 2 , where t is the medium effective thickness. Following the results reported in 21 one can express the geometry dependence of the lumped inductance as π where c 0 is the speed of light in vacuum. By imposing the condition ε′ = 0 GM and assuming a perfect conductor, from eq. (2) we can write the transition frequency as a function of F: In the above formula, t represents the effective optical thickness of the GM taken equal to the dielectric substrate thickness (t = s).
Eq. (3) is plotted in Fig. 3 as a black dash-dotted curve. The analytical expression nearly overlaps the dependence given by the numerical data, showing that one can express the transition of a GM from the plasmonic regime in a very simple and general way. Black spheres represent the experimental values obtained from the real part of the dielectric function (see Fig. 2(b)), in very good agreement with both the analytical and numerical curves. The ω r (f ) curve sets also the transition from the "plasmonic" (homogeneous) to the"photonic" (diffractive) regime (see diagram in Fig. 3)). Above ω r the GM assumes a dielectric behavior (ε′ > 0 GM ), and because of the periodicity specific resonating modes are geometrically selected, enabling to consider the system as composed of separated metallic unit cells 30 .
Transmission peaks above ω r are reported in Fig. 4(a,b). Black curves represent the full wave numerical simulations. In the graphs, HOSPP excitations are labeled as M i GR and M i CB respectively and are bounded at higher where ε ε ε ε = ′ ′ + n / eff md m d is the effective refractive index at the MDI. From eq. (4), it is possible to approximately estimate the value of n eff using the frequency of transmission maxima in the SPP spectrum. In Tables 1 and 2 all states excited in the CB and GR samples respectively, and corresponding to the peaks marked in Fig. 4(a,b), are listed according to the triad ( ) n m m , , eff max x y . The effective refractive index ranges in between 2 and 3, in agreement with the corresponding values found above for ε′ GM and reported in Fig. 2(b).
Although ε′ GM assumes positive values in correspondence of the HOSPP band, the plasmonic features of the patterned metallic layer is preserved. Indeed, by employing the standard dielectric function of the metal-dielectric boundary ε ε ε 4 , we can extract ε′ m using, in place of ε  FR4 and ε  GM , the function retrieved from the experimental data. The real part of the metal layer dielectric function is reported in the inset of Fig. 4,    , as clearly seen as minima in the transmission response approximately at f 1,0 = 0.5THz, f 1,1 = 0.7THz.
Further information on SPP dynamics can be obtained by performing full wave simulations of the electric field spatial distribution of each excited state. In the following figures, the impinging radiation propagates along the z-axis with the electric field oriented parallel to the y-axis.
We first report in Fig. 5 the field distribution of SSPPs for ω < ω r (f = 0.15 THz) for the GR (a) and CB (b) metasurface. We choose the E z component to highlight the dipolar symmetry of the field. In particular, the figure shows the electric field distribution in correspondence of the capacitive line where the GM displays a series of interrupted metallic patches.
For ω > ω r , high order SPPs come into play and show up in the excitation spectra of the GM samples. We have selected the three preeminent maxima in | | ∼ T . In particular, we refer to the peaks M GR 1 , M GR 2 , M GR 4 for the GR and M CB 1 , M CB 2 , M CB 3 for the CB metasurface. The analysis of field distribution of HOSPPs is operated by choosing two different observation planes. In Figs 6, 7 and 8 the z-y plane is selected, so that one can look at the (mostly dipolar) excitations when they propagate along the sample thickness. Instead, in Figs 9 and 10 the HOSPP modes are studied in the x-y plane, to highlight the onset of multipolar features at the MDI.
In Figs 6  . Most of the dipolar field builds up inside the substrate underneath the metallic holes, producing field symmetry above and beyond the conducting layer.
In Fig. 8 the E z intensity along the inductive line (see sketch in the figure) is reported for the GR metasurface. Here the electric field symmetry of HOSPPs differs from what shown along the capacitive line, mostly for the dipolar symmetries of M GR 2 , M GR 4 which appear reversed. The symmetries of dipoles distribution of HOSPP above and beyond the metallic layer directly depend on the different electromagnetic coupling at the sample interface (air and FR4 respectively).
Another interesting information regarding the field spatial distribution in correspondence of each HOSPP frequency is acquired by investigating E z on the metasurface plane, as reported in Figs 9 and 10 for the CB and GR samples respectively. The former always displays a dipolar distribution independently of the mode order. On the    character for the corresponding modes.
The onset of multipolar excitations in GMs corresponds to the behavior observed in metamaterials composed of separated unit cells [33][34][35] , confirming the observation of a diffraction-like response for ω > ω r .

Conclusion
The THz response of grid metasurfaces realized using different patterns of a single square metal patch has been presented. These systems are valuable for investigating the whole family of SPP excitations observable below and above the transition frequency ω r , that separates the plasmonic/metamaterial region from the photonic crystal regime. ω r sets the threshold between the energy necessary to excite either spoof SPPs, for ω < ω r , or high order SPPs, for ω > ω r . This transition frequency can be tuned and spans over one order of magnitude by only changing the filling factor of the GMs. Combining the experimental results and full wave electromagnetic simulations the dependence of the transition frequency on the metasurface filling factor F is retrieved. A simple lumped circuit is applied to model the GM unit cell, and from here an analytical expression of ω r (F) is obtained, depending on the geometrical characteristics of each grid and on the permittivity of the substrate only.
Since the fundamental SSPP mode behaves like in the case of an unpatterned MDI, ω r represents the reference energy determining the upper bound for the homogeneous response of the system. Therefore, it can be used as a quantitative observable alternative to the wavelength threshold condition λ ε ∼ p thr d below which the system discretizes.
The electric field at the cross section of the grid metasurfaces shows that the dipolar distribution of SSPPs is mostly localized at the metal/dielectric boundary, whereas HOSPPs basically present a spatial extension across the GM sample, being responsible for the electromagnetic tunneling and giving rise to enhanced transmission phenomena. Instead, the electric field at the MDI presents the onset of multipolar excitations in correspondence of transmission peaks.  Measurements. Grid metasurfaces are printed over 10 cm × 10 cm FR4 foils using standard electrochemical and galvanic processes. Samples are cut in shape of squares, having an area of 2 × 2 cm 2 . In Fig. 11 a picture of the samples is shown.
The metasurfaces are placed on an aluminum plate with circular holes 15 mm in diameter. Similarly to previous experiments 26,27 , the sample holder can be finely moved in and out the THz beam to acquire in the same measurement session both the signal transmitted through the sample, E s (t), and the reference one in air, E air (t). In order to remove unwanted water absorptions, measurements are carried out in a dry nitrogen environment with humidity level of the order of 0.1% or less. Spectral analysis of samples are performed calculating the transmitted signal =    T E E / s air , where  E is the fast Fourier transform of the time signal.
Electrical model of wire and grid media. The dielectric function of a bare metallic film, having thickness t, can be simply described sampling the film as composed by parallel array of wires (wire medium), whose electrical unit cell is constituted by a series of a resistance R and an inductance L, periodically distributed over a plane with step p. The total impedance of the structure matches with the impedance of the single unit cell ω = + z r i l u , where = r R p / and = l L p / (see Fig. 1(c)). Following the same approach described in 23 , the dielectric function ε  WM of a wire medium can be obtained by using ε where the conductivity σ  is related to the unit cell impedance ~σ = z pt 1/ u . After some algebra, the real and imaginary parts of ε  WM can be written in a Drude-like form as where ω′ = l pt 1/ r 2 is the diluted plasma frequency, τ = l r / is the relaxation time, and ε = ∞ 1 for standard metals. From eq. (5) it is easy to find the root frequency ω r satisfying ε ω Obviously, in a metal wire medium ω ω = ′ r p 2 2 only for τ → ∞ → r ( 0). Extending previous approaches 23,24,37,38 , we qualitatively describe the transformation of a wire medium into a grid metasurface by introducing, in parallel to the previous unit cell, an effective capacitance C periodically distributed with step p (Fig. 1(c)). As a result, the new unit cell impedance is In other terms, the introduction of a capacitance per unit length lowers the metal-dielectric transition frequency, modifying the asymptotic value of permittivity to ε ε ω ω ε . Since in both GMs the holes in the metal are square with side d, the capacitance can be assumed as composed by two armatures with area t·(p-d) and distance (p-d), so that ε = C t d and ε = = c pt C t / / d . Following 21 , one can empirically extend the geometrical dependence of the lumped inductance given for a wire medium to the case of a grid (two-dimensional) structure, using the expression π = l F ln(1/ )/ . In such a way, one can write the effective plasma frequency as where we assume that the metal behaves as a perfect conductor with ( τ = →∞ r 0 ( )). In other words, the transmitted radiation does not feel the metal dissipation because the conducting layer is thicker than the skin depth.
Following the schemes reported in Fig. 11(a,b), the filling factor F for each grid can be easily related to η = d/p