## Abstract

The THz response of slit structures and split-ring resonators (SRRs) featuring extremely small gaps on the micro- or nanoscale is investigated numerically. Both structures exhibit strong field enhancement in the gap region due to light-induced current flows and capacitive charging across the gap. Whereas nanoslits allow for broadband enhancement the resonant behavior of the SRRs leads to narrowband amplification and results in significantly higher field enhancement factors reaching several 10,000. This property is particularly beneficial for the realization of nonlinear THz experiments which is exemplarily demonstrated by a second harmonic generation process in a nonlinear substrate material. Positioning nanostructures on top of the substrate is found to result in a significant increase of the generation efficiency for the frequency doubled component.

©2011 Optical Society of America

## 1. Introduction

The last decades have seen significant progress in the development of terahertz (THz) technologies in a wide variety of fields such as chemical recognition, material inspection, or security control. One of the major and yet unmet limitations is that, compared to the optical regime, the pulse energies supplied by current THz sources are still rather limited. This is particularly problematic for the realization of nonlinear THz experiments which have been a topic of considerable interest in recent years. The highest average THz power is currently available from large-scale electron particle accelerators [1]. Using table-top sources generation of high energy THz pulses has been demonstrated using for example large area photoconductive switches [2], frequency mixing in a laser generated plasma [3], or optical rectification in nonlinear crystals [4].

In order to extend THz experiments into the nonlinear regime the quantity to be optimized is essentially the electric field strength *E* since a nonlinear process of order *n* scales with *E ^{n}*. The field strength is linked to the pulse energy

*Q*through the relation $E\propto \sqrt{Q/\Delta \tau \xb7A}$ where Δ

*τ*is the temporal width of the pulses and

*A*denotes the cross-sectional beam area. The available pulse energy is naturally limited by the THz system at hand and the pulses are typically already single-cycle so that one cannot increase the field strength by adjusting the parameters

*Q*and Δ

*τ*. The cross-sectional beam area can be minimized through tight focusing which is typically realized using e.g. paraboloidal mirrors. However, the diffraction limit imposes a lower boundary on

*A*thereby restricting the obtainable field strength. This can be overcome by using metallic nanostructures that collect the incident radiation and focus it in a subwave-length volume leading to strong field enhancement. This concept was introduced by Seo and coworkers [5] who have shown that a single nanoslit in a thin gold film (Fig. 1(a)) may act as a nanocapacitor. The incident radiation induces a current flow on the metal surface that leads to an accumulation of charge carriers in the gap region. This capacitive charging in turn results in an in-gap enhancement of the electric field by orders of magnitude. For a 70 nm wide slit and a frequency of 0.1 THz field enhancement factors on the order of 1,000 were reported by the authors. A slightly different concept using arrays of slits with different lengths has recently also been applied to obtain broadband enhanced THz radiation [6].

Here, we propose a different geometry in order to further increase the obtainable field enhancement. The structure consists of split-ring resonators (SRRs) [7] that are resonant in the THz regime and which, compared to their side length, feature extremely small gaps on the micro- or nanoscale (Fig. 1(b)). At the fundamental SRR resonance, which is also denoted LC resonance due to its analogy with an LC resonator, strong circulating currents are excited in the ring that lead to a build up of charge across the gap [8]. Due to the resonant behavior the resulting field strengths in the gap region are significantly higher than for the non-resonant nanoslits. For SRRs with gap widths on the nanoscale giant enhancement factors of several 10,000 are obtained. In order to maximize the integrated nonlinear response the overall size of the volume comprising the high field strengths also plays an important role. We therefore investigate how the nonlinear response can be maximized by adjusting the structural parameters of both nanoslits and microgap SRRs. Finally, we exemplarily use a second harmonic generation (SHG) process in a nonlinear substrate material to compare the strength of the nonlinearities induced by the two structures.

## 2. Simulations

The numerical modeling is based on the finite element method (FEM) [9] using a commercial software package [10]. The simulations were carried out in both the frequency- and the time-domain. The frequency dependent simulations were used in order to determine the obtainable field enhancement in single slits and SRRs. This approach has the advantage that the spectral dependence of the metal may be included using a complex conductivity as obtained from the Drude model where we used the values of bulk gold [11] for the plasma frequency and damping constant. Even though thin metallic films have been shown to exhibit a smaller conductivity than the bulk materials [12] we observe that the effect of a tenfold increased damping constant on the field enhancement is negligible. Modeling of the slit structures could be carried out in two dimensions assuming infinite slit lengths whereas for the SRRs a three dimensional geometry was employed. The slit (SRR) structures were positioned in the center of a rectangular (box shaped) simulation domain and the refractive index of the surrounding dielectric was set to unity. A plane harmonic wave having normal incidence onto the structure plane was launched from one of the boundaries and the frequency *ν* was scanned parametrically. For all other boundaries scattering boundary conditions were chosen. In order to include nonlinear effects such as SHG, however, we analyze the structure’s response in the time-domain. For this purpose a multi-cycle THz transient was employed. Also, for the examination of the SHG process we chose to investigate the interaction of the THz pulse with slit and SRR arrays rather than single structures. For this purpose the corresponding boundaries were set to periodic conditions and the size of the simulation domain was chosen according to the periodicity.

## 3. Results

#### 3.1. Obtainable field enhancement

### 3.1.1. Nanoslit structures

We begin by investigating the response of single nanoslits as schematically shown in Fig. 1(a). The thickness of the gold film was set to *h* = 60 nm, the slit width *D* was varied from 50 μm down to 40 nm, and the incident wave was polarized perpendicular to the slit’s long axis. From the frequency dependent simulations we determine the field enhancement factor

*E*

_{gap}(

*ν*) denotes the absolute amplitude of the electric field in mid-gap and

*E*

_{inc}(

*ν*) the one of the incident electric field. We find that this is an adequate measure for the field enhancement allowing for convenient comparison between different geometries even though small inhomogeneities across the gap are observed. Seo and coworkers [5] pointed out that the field enhancement obtainable in such structures depends on both the wavelength

*λ*and the slit width

*D*as shown in Fig. 2(a). We find that for slits with

*D*≪

*λ*the field enhancement scales approximately linear with the ratio

*λ /D*. This is shown in Fig. 2(b) where we plot

*F*as a function of this ratio for four different frequency components. The different curves overlap which shows that the funneling is a purely geometric effect and that the smaller the slit on the scale of the wavelength the larger is the obtainable field enhancement. For the smallest considered slit width and frequency (

*D*= 40 nm,

*ν*= 100 GHz) a maximum enhancement factor approaching 3000 is obtained.

### 3.1.2. Nanogap split-ring resonators

Figure 1(b) gives a schematic representation of the SRR structures under study. The spectral position of the SRR resonances can be controlled by adjusting the side length *L* [13, 14]. Here, we chose *L* = 300 μm to obtain resonances at the low frequency side of the THz range. The wire width was set to *w* = 5 μm, the structure height was *h* = 1 μm, and the gap width *d* was varied from 100 μm down to 100 nm. The LC resonance can be excited through either an electric field component across the SRR gap, a magnetic field component perpendicular to the SRR plane, or by a combination of the previous two [15, 16]. We start by considering *E*-field excitation where the *E*, *H*, and *k* triad of the incoming field is oriented along the *x*, *y*, and *z* axes. In Fig. 3(a)–3(c) we exemplarily show the field enhancement *F*(*ν*) obtained for three SRRs with decreasing gap width.Within the considered frequency range the curves exhibit a single peak corresponding to the fundamental SRR resonance. With decreasing gap width the resonance is redshifted which can be qualitatively understood by the analogy between a SRR and a simple LC resonator [8]. A smaller gap width translates into a larger capacitance which in turn results in a reduced resonance frequency. Also, for decreasing *d* the in-gap field strength is found to be significantly enhanced. In analogy to the nanoslits we find that the field enhancement scales linearly with the ratio *λ _{r}*/

*d*where

*λ*denotes the wavelength at resonance. This is shown in Fig. 3(d) where the red line is a linear fit to the data points. For the smallest considered gap width of

_{r}*d*= 100 nm giant enhancement factors approaching 40,000 are obtained which is more than an order of magnitude larger than the maximum field enhancement obtained for the nanoslits. This is explained by the resonant behavior where the oscillating currents excited in the SRR are much stronger than the non-resonant current flow induced on the metallic surface of the slit structures.

Note, that with the fixed values of *h* and *w* the aspect ratio of the gap volume becomes more and more extreme for decreasing gap width requiring an increasing amount of mesh points which in turn results in a significant increase of the computational demands. Simulating gap sizes smaller than 100 nm could therefore not be realized. We expect, however, that the inverse scaling of the field enhancement with *d* extends to smaller gap widths so that even larger field strengths could be realized using e.g. SRRs with 40 nm gap width. Here, the threshold for breakdown fields (on the order of 10^{8} V/cm) imposes an upper limit on the obtainable in-gap field strengths. The fabrication of nanostructures featuring such extreme aspect ratios, on the other hand, has already been demonstrated [17].

It has been shown previously [18] that compared to only *E*-field excitation the LC resonance may be excited more strongly using a combination of *E*- and *H*-field excitation. Using this excitation scheme, where the *E*, *H*, and *k* triad is oriented along the *x*, *-z*, and *y* axes, we find that the in-gap field enhancement can be increased even further. This is exemplarily shown in Fig. 3(e) for a SRR with 1 μm gap width. By simply rotating the structure the peak values of the in-gap fields can approximately be increased by another factor of 1.3.

#### 3.2. Optimizing the geometry for maximum nonlinear response

### 3.2.1. Slit structures arranged in one dimensional arrays

The integrated nonlinear response depends not only on the obtainable field strength but also on the size of the volume comprising those high field strengths. Since for a single nanoslit the total volume and, thus, the overall transmission is rather small the obvious approach is to arrange the slits into one dimensional arrays with periodicity *p* as shown in Fig. 1(c). Since decreasing the periodicity increases the slit density the transmission is expected to scale inversely with *p*. For small periodicities, however, neighboring slits start sharing the incident field energy so that, compared to single slits, the obtainable field enhancement is reduced [19]. The use of nanoslit arrays is consequently accompanied by a tradeoff between an increased volume and a decreased field enhancement. For the realization of nonlinear *χ*
^{(2)} and *χ*
^{(3)} effects the properties to be optimized are the ratios *F*
^{2}(*p*)/*p* and *F*
^{3}(*p*)/*p*, respectively. Their dependence on the periodicity is shown in Fig. 4(a) as obtained for 40 nm wide slits at a frequency of 500 GHz. For the two curves maxima are obtained at different values of *p* indicating that the ideal periodicity depends on the order of the nonlinearity to be optimized. In addition, the behavior also exhibits a strong frequency dependence as exemplarily shown for a *χ*
^{(2)} process in Fig. 4(b).

### 3.2.2. Split-ring resonators featuring extended capacitive faces

Due to the ring geometry of SRRs, at least in planar arrays, the possibility of increasing the nonlinear response by decreasing the periodicity is limited. We therefore follow another approach to increase the size of the high field strength volume that is based on adding capacitive faces of length *s* to the gap [20] as shown in Fig. 5(a). The electromagnetic energy stored in the SRR and, thus, the amount of charge carriers accumulating in the gap region is finite. For increasing split length the charge accumulations are therefore distributed over a larger distance so that the field enhancement is found to decrease. This is shown in Fig. 5(a) for single SRRs with *d* = 1 μm gap width. As in the case of nanoslits this implies a tradeoff between an increased volume and a decreased field enhancement. For the realization of *χ*
^{(2)} and *χ*
^{(3)} effects the integrals

*s*is shown in Fig. 5(b). For

*U*

^{(2)}and

*U*

^{(3)}maxima at different values of

*s*are obtained which again shows that the ideal split length depends on the order of the nonlinearity.

#### 3.3. Nanostructure induced second harmonic generation in lithium tantalate

We exemplarily consider a second harmonic generation (SHG) process in lithium tantalate (LiTaO_{3}). The purpose of this study is not to optimize the conversion efficiency of the SHG process but to demonstrate that the structures under study can be used to induce THz nonlinearities in adjacent materials. We note, that similar concepts have recently been used to control the THz transmission through a VO_{2} film [21] or to produce extreme ultraviolet light through high harmonic generation [22]. LiTaO_{3} is a trigonal 3*m* crystal having one extraordinary axis. The incident field is polarized along this axis and, for simplicity, we neglect the birefringence in the simulations. According to [23] the refractive index of LiTaO_{3} at THz frequencies can then be set to n = 6.15. The nonlinearity was included by introducing a nonlinear polarization *P*
^{(2)} = *ε*
_{0}
*χ*
^{(2)}
*E*
^{2} into the governing time dependent equation. For *χ*
^{(2)} at THz frequencies a constant numerical value could be derived from the potential energy surface as shown in Appendix A. For simplicity the imaginary part has been neglected. The central frequency of the exciting multi-cycle pulse was 138 GHz and the peak field strength was set to *E*
_{inc} = 10 kV/cm. Note, that this value is about an order of magnitude smaller than what has already been demonstrated using table-top THz sources [4].

### 3.3.1. Nanoslit arrays

We consider a nanoslit array deposited on top of a flat LiTaO_{3} substrate. The slit width was 40 nm and a periodicity of 30 μm was used optimized for maximum SHG efficiency at the center frequency. Note, that the introduction of the LiTaO_{3} substrate leads to a modification of the dielectric environment which results in an effectively smaller wavelength *λ*
_{eff} = *λ /n*
_{eff} where *n*
_{eff} denotes the effective refractive index. According to our previous findings this leads to a decrease of the field enhancement as compared to a freestanding structure. The black curve in Fig. 6(a) shows a reference spectrum of the incident pulse transmitted through an unstructured LiTaO_{3} substrate. The nonlinearity of LiTaO_{3} together with the given field strength of the incident pulse is small so that hardly any SHG signal can be observed. The red curve, on the other hand, shows the corresponding spectrum of a pulse transmitted through the nanoslit array on the substrate. The induced field enhancement leads to a much stronger SHG signal even though the nonlinear substrate material was only present below and not in the gap. Figure 6(b) shows the field enhancement in a *xz*-slice cutting through the gap. The high field strength region extends a few tens of nanometers into the LiTaO_{3} substrate. Even though, the volume comprising the high field strengths within LiTaO_{3} is rather small, the field enhancement within this volume is large enough to cause a significant increase of the SHG efficiency. The generation efficiency can be scaled up even further if the nonlinear material is directly inserted into the slit, for example by using liquids [24]. For demonstration we have also analyzed the case of a nanoslit array on substrate where the slit volume was additionally filled with LiTaO_{3}. The obtained transmission spectrum is included in Fig. 6(a) as the green curve. Compared to SRRs with air-filled gaps (red curve) the SHG signal is increased by another factor of 2.3.

### 3.3.2. Microgap split-ring resonator arrays

For the SRRs the presence of the LiTaO_{3} substrate results in a redshift of the LC resonance [13] due to the SRR responding to the effectively reduced wavelength *λ*
_{eff}. This redshift was compensated by choosing a correspondingly smaller side length of *L* = 70 μm for which the SRR resonance coincides with the center frequency of the incident spectrum. The SRRs were arranged in a two dimensional array and the periodicity was set to 100 μm. Since time-dependent simulations in three dimensions are particularly time consuming an intermediate gap width of *d* =1 μm was chosen as a compromise between obtainable field enhancement and a manageable amount of mesh points.

Figure 6(c) shows the spectrum of the incident pulse transmitted through the sole substrate (black curve, reference), a microgap SRR array deposited on LiTaO_{3} (red curve), and the SRR array on substrate where the gap volume was additionally filled with LiTaO_{3} (green curve). The exhibited behavior qualitatively resembles the one of the nanoslits except that the amplitudes of the generated second harmonic spectra are much higher. The slight discrepancy between the frequency doubled components of the black curves in Fig. 6(a) and 6(c) is assigned to different numerical noise levels in the two and three dimensional simulations. In Fig. 7(a) we plot the peak amplitudes of the SHG spectra obtained with the unfilled structures as a function of the incident field strength. Here, the data points were referenced to the peak amplitude of the incident spectrum at the highest considered field strength. From the log-log plot we obtain a slope of 2 for both curves indicating a quadratic dependence on the field strength as expected for a *χ*
^{(2)} process. Much higher SHG efficiencies are obtained for the SRRs even though an optimized nanoslit geometry was employed whereas the SRR conversion efficiency can still be increased by (i) using narrower gaps, (ii) adding extended capacitive faces, and (iii) using a combined *E*- and *H*-field excitation. Note, however, that the strong field concentration may cause photorefractive damage so that the natural damage threshold imposes an upper limit on the obtainable nonlinearites.

In Fig. 7(b) we plot the normalized SHG spectra of the slit and SRR structures on substrate (corresponding to the red curves in Fig. 6(a) and 6(c)). In addition, we include the ideal frequency doubled spectrum as calculated from the incident waveform (dashed black curve). The spectrum obtained with the SRR samples is somewhat narrower than the calculated spectrum which is explained by the SRRs providing field enhancement only at their resonances. The SHG spectrum of the nanoslits, on the other hand, has approximately the same width as the calculated spectrum except that it is slightly redshifted which is attributed to the field enhancement scaling with *λ*. The difference in the behavior exhibited by the two structures can consequently be used to realize either broadband or narrowband conversion. The SRRs are thereby of particular interest because the narrowband conversion property allows for tunability of the nonlinearity by adjusting the dielectric environment or the structural parameters. This behavior can be seen in Fig. 6(c) where the red and green curves exhibit a dent in the incident spectrum corresponding to the fundamental SRR resonance. Filling the gap with LiTaO_{3} results in a small redshift which is reproduced in the second harmonic spectrum. For better visibility of this effect the corresponding SHG spectra are shown again in Fig. 7(c) where the curves were normalized. The redshift induced by the modification of the dielectric environment can clearly be seen. The frequency position and linewidth of SRR resonances can be tailored by adjusting the spilt-ring’s geometry [13] or by introducing a coupling between neighboring rings [14, 25]. This property can in turn be used to realize structural tunability of the induced nonlinearities as demonstrated in Fig. 7(c). The dashed black line represents the SHG spectrum obtained with smaller SRRs positioned on LiTaO_{3}. The SRR side length *L* was reduced by 10% which results in a blueshift of the resonance that is consequently reproduced in the SHG spectrum.

## 4. Conclusion

We have analyzed the response of slit structures as well as THz split-ring resonators featuring extremely small gaps on the micro- or nanoscale. Both structures can be used to enhance the electric field by orders of magnitude. Whereas broadband enhancement is exhibited by the nanoslits the resonant nature of the SRRs leads to a narrowband behavior. Due to strong oscillatory currents induced in the ring at resonance giant enhancement factors approaching 40,000 were observed. These structures are valuable tools for the realization of nonlinear THz experiments and, in order to maximize the nonlinear response, we have investigated the interplay between the size of the high field strength volume and the obtainable field enhancement. Application of the structures was exemplarily demonstrated by a second harmonic generation process in a nonlinear substrate material. For microgap SRRs deposited onto the substrate significantly higher conversion efficiencies were obtained than for nanoslit arrays and suggestions were made on how the SHG efficiency can be increased even further. SRRs additionally allow for tailoring of the induced nonlinearity due to the narrowband conversion efficiency together with the possibility of structural resonance tuning. Also, whereas for higher harmonic generation in bulk crystals a correct crystal orientation is crucial in order to obtain phase matching, this condition is rendered irrelevant if nanostructures are used due to highly localized high field strength regions. Based on the enormous field enhancement, the investigated structures are expected to find wide ranging application in the field of nonlinear spectroscopy or functional THz elements based on various nonlinearities such as the Kerr or Pockel’s effect.

## Appendix A. Nonlinearity of LiTaO_{3} at THz frequencies

## A.1. Potential energy surface

The potential energy surface of LiTaO_{3} may well be represented by [26, 27]

*Q*is the net ionic displacement,

*M*is the oscillator mass, −

*is the frequency of the transverse optical phonon mode, and*

_{T}*V*

_{0}is the value of the potential at its minima (well depth).

## A.2. Expressions for the nonlinear susceptibilities

In nonlinear optics the relation between the polarization density *P* and the electric field *E* is usually expanded in terms of a power series [28]

*ε*

_{0}is the vacuum permittivity and

*χ*

^{(}

^{n}^{)}with

*n*= 1,2,3, … denotes the susceptibility of different order. In order to include nonlinear effects of second and third order in our simulations we need to obtain numerical values for

*χ*

^{(2)}and

*χ*

^{(3)}of LiTaO

_{3}at THz frequencies. We start by assuming that the crystal consists of independent harmonic oscillators that couple to the electric field through their dipole moment [23]. The equation of motion is then given by

*Q*is the ionic displacement, Γ is a phenomenological damping constant and

*q*denotes the effective charge. Using a Taylor expansion of

*V*(

*Q*) where we consider only terms up to third order we obtain where We now follow the approach in [28] and use a procedure analogous to Rayleigh-Schrödinger perturbation theory. The electric field

*E*(

*t*) is replaced by

*λ E*(

*t*) and

*Q*is expanded in a power series

*Q*=

*λQ*

^{(1)}+

*λ*

^{2}

*Q*

^{(2)}+

*λ*

^{3}

*Q*

^{(3)}. The terms proportional to the different orders of

*λ*

^{(}

^{n}^{)}must satisfy Eq. (6) separately so that we obtain a set of three differential equations

*n*is linked to the displacement

*Q*

^{(}

^{n}^{)}and the electric field

*E*through Assuming harmonic time dependence of the form we can obtain solutions for Eq. (8)–(10) which can be combined with Eq. (11) to obtain expressions for the susceptibilities of different order and

## A.3. Numerical values of constants used for the calculation of the nonlinear susceptibilities

In order to calculate *χ*
^{(2)} and *χ*
^{(3)} we need numerical values for *ω _{T}*, Γ,

*V*

_{0},

*N*,

*M*, and

*q*. The first three are found in the literature and are given in Table 1 where we assume a polarization along the extraordinary axis. The oscillator density

*N*is calculated by considering the inverse of a 10-atom unit cell with hexagonal lattice in a rhombohedral setting [29]

*a*and

*c*can be taken from x-ray and neutron scattering data [30, 31]. The oscillator mass

*M*is the reduced mass given by [29] where

*m*

*refers to the Li, Nb, and O atoms within a unit cell. The values calculated for*

_{j}*M*and

*N*are included in Table 1. In order to also obtain a value for

*q*we use the permittivity

*ε*and

_{s}*ε*

_{∞}are found in the literature. The corresponding values along the extraordinary axis as well as the calculated value for

*q*are included in Table 1. The nonlinear susceptibilities are then calculated to

*χ*

^{(2)}= (1.75 –

*ι ·*0.08) · 10

^{−}

^{8}m/V and

*χ*

^{(3)}= (2.1 –

*ι ·*0.2) · 10

^{−18}m

^{2}/V

^{2}. Note, that the imaginary part of

*χ*

^{(2)}causes an absolute phase shift of the second harmonic field whereas the imaginary part of

*χ*

^{(3)}leads to two-photon absorption, i.e. loss.

## Acknowledgment

This work was supported by the Swiss National Research Foundation (SNSF) through the NCCR-MUST and through project 200020-119934. A.B. acknowledges financial support by Limat.

## References and links

**1. **G. P. Williams, “Filling the THz gap-high power sources and applications,” Rep. Prog. Phys. **69**, 301 (2006). [CrossRef]

**2. **E. Budiarto, J. Margolies, S. Jeong, and J. Song, “High-intensity terahertz pulses at 1-kHz repetition rate,” IEEE J. Quantum Electron. **32**, 1839 (1996).

**3. **T. Bartel, P. Gaal, K. Reimann, M. Woerner, and T. Elsaesser, “Generation of single-cycle THz transients with high electric-field amplitudes,” Opt. Lett. **30**, 2805–2807 (2005). [CrossRef] [PubMed]

**4. **K. Yeh, M. C. Hoffmann, J. Hebling, and K. A. Nelson, “Generation of 10 μJ ultrashort terahertz pulses by optical rectification,” Appl. Phys. Lett. **90**, 171121 (2007). [CrossRef]

**5. **M. A. Seo, H. R. Park, S. M. Koo, D. J. Park, J. H. Kang, O. K. Suwal, S. S. Choi, P. C. M. Planken, G. S. Park, N. K. Park, Q. H. Park, and D. S. Kim, “Terahertz field enhancement by a metallic nano slit operating beyond the skin-depth limit,” Nat. Photonics **3**, 152–156 (2009). [CrossRef]

**6. **H. R. Park, Y. M. Park, H. S. Kim, J. S. Kyoung, M. A. Seo, D. J. Park, Y. H. Ahn, K. J. Ahn, and D. S. Kim, “Terahertz nanoresonators: Giant field enhancement and ultrabroadband performance,” Appl. Phys. Lett. **96**, 121106 (2010). [CrossRef]

**7. **J. B. Pendry, A. J. Holden, D. J. Robbins, and W. J. Stewart, “Magnetism from Conductors and Enhanced Nonlinear Phenomena,” IEEE Trans. Microwave Theory Tech. **47**, 2075–2084 (1999). [CrossRef]

**8. **S. Linden, C. Enkrich, M. Wegener, J. Zhou, T. Koschny, and C. M. Soukoulis, “Magnetic Response of Metamaterials at 100 Terahertz,” Science **306**, 1351–1353 (2004). [CrossRef] [PubMed]

**9. **J. Jin, *The Finite Element Method in Electromagnetics*, 2nd ed. (Wiley-IEEE Press, 2002).

**10. **COMSOL Multiphysics 3.5.

**11. **M. A. Ordal, R. J. Bell, R. W. Alexander, L. L. Long, and M. R. Querry, “Optical properties of fourteen metals in the infrared and far infrared: Al, Co, Cu, Au, Fe, Pb, Mo, Ni, Pd, Pt, Ag, Ti, V, and W,” Appl. Opt. **24**, 4493–4499 (1985). [CrossRef] [PubMed]

**12. **N. Laman and D. Grischkowsky, “Terahertz conductivity of thin metal films,” Appl. Phys. Lett. **93**, 051105 (2008). [CrossRef]

**13. **M. Kafesaki, T. Koschny, R. S. Penciu, T. F. Gundogdu, E. N. Economou, and C. M. Soukoulis, “Left-handed metamaterials: detailed numerical studies of the transmission properties,” J. Opt. A, Pure Appl. Opt. **7**, S12–S22 (2005). [CrossRef]

**14. **A. Bitzer, J. Wallauer, H. Helm, H. Merbold, T. Feurer, and M. Walther, “Lattice modes mediate radiative coupling in metamaterial arrays,” Opt. Express **17**, 22108–22113 (2009). [CrossRef] [PubMed]

**15. **P. Gay-Balmaz and O. J. F. Martin, “Electromagnetic resonances in individual and coupled split-ring resonators,” J. Appl. Phys. **92**, 2929–2936 (2002). [CrossRef]

**16. **N. Katsarakis, T. Koschny, M. Kafesaki, E. N. Economou, and C. M. Soukoulis, “Electric coupling to the magnetic resonance of split ring resonators,” Appl. Phys. Lett. **84**, 2943–2945 (2004). [CrossRef]

**17. **S. Gorelick, V. A. Guzenko, J. Vila-Comamala, and C. David, “Direct e-beam writing of dense and high aspect ratio nanostructures in thick layers of PMMA for electroplating,” Nanotechnology **21**, 295303 (2010). [CrossRef] [PubMed]

**18. **J. García-García, F. Martín, J. D. Baena, R. Marqués, and L. Jelinek, “On the resonances and polarizabilities of split ring resonators,” J. Appl. Phys. **98**, 033103 (2005). [CrossRef]

**19. **M. Shalaby, H. Merbold, M. Peccianti, L. Razzari, G. Sharma, R. Morandotti, T. Ozaki, T. Feurer, A. Weber, L. Heyderman, H. Sigg, and B. Patterson, “Concurrent field enhancement and high transmission of THz radiation in nanoslit arrays,” in preparation (2011).

**20. **D. Schurig, J. J. Mock, B. J. Justice, S. A. Cummer, J. B. Pendry, A. F. Starr, and D. R. Smith, “Metamaterial Electromagnetic Cloak at Microwave Frequencies,” Science **314**, 977–980 (2006). [CrossRef] [PubMed]

**21. **J. Kyoung, M. Seo, H. Park, S. Koo, H. sun Kim, Y. Park, B.-J. Kim, K. Ahn, N. Park, H.-T. Kim, and D.-S. Kim, “Giant nonlinear response of terahertz nanoresonators on VO_{2} thin film,” Opt. Express **18**, 16452–16459 (2010). [CrossRef] [PubMed]

**22. **S. Kim, J. Jin, Y.-J. Kim, I.-Y. Park, Y. Kim, and S.-W. Kim, “High-harmonic generation by resonant plasmon field enhancement,” Nature **453**, 757–760 (2008). [CrossRef] [PubMed]

**23. **T. Feurer, N. S. Stoyanov, D. W. Ward, J. C. Vaughan, E. R. Satz, and K. A. Nelson, “Terahertz Polaritonics,” Annu. Rev. Mater. Res. **37**, 317–350 (2007). [CrossRef]

**24. **M. C. Hoffmann, N. C. Brandt, H. Y. Hwang, K. Yeh, and K. A. Nelson, “Terahertz Kerr effect,” Appl. Phys. Lett. **95**, 231105 (2009). [CrossRef]

**25. **H. Liu, D. A. Genov, D. M. Wu, Y. M. Liu, Z. W. Liu, C. Sun, S. N. Zhu, and X. Zhang, “Magnetic plasmon hybridization and optical activity at optical frequencies in metallic nanostructures,” Phys. Rev. B **76**, 073101 (2007). [CrossRef]

**26. **I. Inbar and R. E. Cohen, “Comparison of the electronic structures and energetics of ferroelectric LiNbO_{3} and LiTaO_{3},” Phys. Rev. B **53**, 1193 (1996). [CrossRef]

**27. **V. Romero-Rochin, R. M. Koehl, C. J. Brennan, and K. A. Nelson, “Anharmonic phonon-polariton excitation through impulsive stimulated Raman scattering and detection through wave vector overtone spectroscopy: theory and comparison to experiments on lithium tantalate,” J. Chem. Phys. **11**, 3559 (1999). [CrossRef]

**28. **R. W. Boyd, *Nonlinear Optics*, 1st ed. (Academic Press, 1992).

**29. **D. W. Ward, “Polaritonics: An Intermediate Regime Between Electronics and Photonics,” Ph.D. thesis, Department of Chemistry – Massachusetts Institute of Technology (2005).

**30. **H. Boysen and F. Altorfer, “A neutron powder investigation of the high-temperature structure and phase transition in LiNbO_{3},” Acta Crystallogr., Sect. B **50**, 405 (1994). [CrossRef]

**31. **R. Hsu, E. N. Maslen, D. du Boulay, and N. Ishizawa, “Synchrotron X-ray Studies of LiNbO_{3} and LiTaO_{3},” Acta Crystallogr., Sect. B **53**, 420 (1997). [CrossRef]