Over the last couple of decades, particular attention has been given to the terahertz (THz) frequency band, ranging from 0.1 to15 THz, as it offers a great variety of important applications, including medical imaging [1,2], non-invasive security scanning [3], wireless communications with up to terabit-per-second data rates [4,5], and a plethora of applications in material sciences [6–8].

However, conventional THz generation systems suffer from some limitations that do not allow us to fully exploit this spectral regime. THz quantum cascade lasers operate at cryogenic temperatures [9,10], while Schottky diode multipliers and Gunn oscillators are limited by transit times and resistance-capacitive effects that limit their bandwidths and conversion efficiencies [11]. Difference frequency generation (DFG) and parametric oscillations in inorganic crystals, such as $\mathrm{ZnTe}/{\mathrm{LiNbO}}_3/\mathrm{GaAs}/\mathrm{GaP}$ pumped by femtosecond laser pulses, offer an alternative approach in the generation of THz emission [12]. Yet, these setups are usually bulky and require high-intensity pump lasers, while their produced spectra show gaps in the Reststrahlen region [13].

Recently, it was shown experimentally that THz radiation can be generated through optical rectification in nonlinear metasurfaces (NLMSs) constructed from split ring resonators (SRRs) [14]. These ultrathin plasmonic metasurfaces have been shown to provide significant second-order nonlinear susceptibilities [15–17] and THz emission amplitudes that are comparable to those achieved by thousand times thicker ZnTe crystals [14,18]. Moreover, it was shown that the ability to engineer the arrangement of nonlinear building blocks on the metasurface allows us to control the spatiotemporal properties of the emitted THz waves [18], or even enhance the generated THz emission [19].

In order to understand the generation process, and to be able to design novel metamaterial-based THz sources and components, there is a need to develop reliable numerical models of the THz generation by nonlinear metamaterials. Previous works address this need by treating the electron gas using a time-domain expansion of the Maxwell-Hydrodynamic model [20,21]. This enables the decrease in computational requirements that a full quantum approach would need. In this letter, we expand the hydrodynamic model perturbatively [22–25], in a way that allows us to describe the second-order nonlinearity that arises in metals in an equivalent surface current density.

The proposed Maxwell–Hydrodynamic model investigates the DFG and is simulated in the frequency domain with a commercial finite element method (FEM) software. The model assumes that the electromagnetic field ($\mathit{E}$, $\mathit{B}$) induced motion of free electrons inside a thin metal layer can be described by Euler’s equation [22–25]:

After expanding the fields’ expressions and applying the boundary conditions [23], DFG is formulated by

The effective nonlinear surface current density is defined as $K_{\mathrm{NL}}=i{\omega }_{\mathrm{gen}}{\int }_{\xi =0}^l{\mathit{P}}_{\mathrm{NL}}\mathrm{d}\mathit{r}$, where ${\omega }_{\mathrm{gen}}$ is the generated frequency. The integration is performed across the surface layer of the nanoparticles in a region $0<\xi \le l$, where the electron gas pressure is significant for the formation of quadratic nonlinear currents. Similar to previous works [27], the above integral can be calculated using the simplified formulas ${\int }_{\xi =0}^l(\nabla ·{\mathit{P}}_1){\mathit{P}}_2\mathrm{d}\mathit{r}=-\hat{\mathit{t}}{P}_1^{\perp }{P}_2^{\parallel }-\hat{\mathit{n}}\frac{1}{2}{P}_1^{\perp }{P}_2^{\perp }$ and ${\int }_{\xi =0}^l({\mathit{P}}_1·\nabla ){\mathit{P}}_2\mathrm{d}\mathit{r}=-\hat{\mathit{n}}\frac{1}{2}{P}_1^{\perp }{P}_2^{\perp }$, in which $\hat{\mathit{n}}$ is the normal vector on the metallic boundaries, and $\hat{\mathit{t}}$ points in the direction of $\hat{\mathit{n}}\times {\mathit{P}}_{\mathrm{NL}}$. Finally, the expression of the nonlinear surface currents is presented below:

Here we use this approach to investigate the role of the DFG process in SRRs, resulting in THz emission [14]. The building block of the NLMS is schematically presented in Fig. 1(a). The thickness of the SRRs is 40 nm, and the plasma parameters are taken to be $\gamma =10.68\times {10}^{13}\mathrm{rad}/\mathrm{s}$, $n_0=5.92\times {10}^{28}{\mathrm{m}}^{-3}$ [20], and the permittivity of gold is taken from Johnson and Christy [28]. Rounded corners of the SRRs with a radius of 5 nm are considered, in order to reduce numerical artifacts due to the field localization near the metallic corners. The metasurface has a square periodicity of 380 nm and is placed on top of a glass substrate with a refractive index of 1.5. Periodic boundary conditions are applied in $x$- and $y$- directions and absorbing boundary conditions in the $z$- direction. Figure 1(b) shows the linear response spectra when the metasurface is illuminated at normal incidence by an $x$- polarized plane wave.

 

Fig. 1. (a) Schematic illustration of the unit cell of the examined metasurface with a periodicity of 380 nm. The SRR’s thickness is 40 nm. (b) Calculated linear response spectrum when the metasurface is illuminated at normal incidence by an infrared laser pulse polarized along the x- direction. R, reflectance; T, transmittance; A, absorbance.

Download Full Size | PDF

For the THz emission studies, we consider illumination with a femtosecond near-infrared pump pulse. The incident pump pulse is described as $E(\omega )=E_o\exp (-{\tau }^2\frac{{(\omega -{\omega }_0)}^2}{8\log 2})$, where $\tau$ is the temporal width of the pulse; its amplitude is $E_o=2\times {10}^7\mathrm{V}/\mathrm{m}$, and the carrier angular frequency ${\omega }_0=1.257\times {10}^{15}\mathrm{rad}/\mathrm{s}$. The structural parameters of our meta-atoms are selected in a way that a magnetic dipole resonance is excited near the central carrier frequency of the pump (200 THz), enhancing the second-order nonlinearities that arise due to the symmetry breaking of the SRRs [27]. These generate nonlinear surface current distributions on the SRRs that follow the described Maxwell–Hydrodynamic model. The process of frequency generation is simulated using COMSOL Multiphysics, which is a FEM solver, and the problem is divided into three study steps. In the first two steps, the response to angular frequencies ${\omega }_1$ and ${\omega }_2$ is simulated [Fig. 2(a)]. These are used to calculate the surface currents that form on the resonators according to Eq. (3). These currents act as radiating sources, emitting at ${\omega }_1-{\omega }_2$.

 

Fig. 2. (a) Schematic of the studied nonlinear optical metasurface. The metasurface is illuminated at normal incidence by a femtosecond pump beam with frequency components ${\omega }_1$ and ${\omega }_2$, which are within the pump bandwidth. The generated DFG signal at ${\omega }_3={\omega }_1-{\omega }_2$ is studied in reflection. The spatial distribution of the (b) linear and (c) nonlinear surface currents.

Download Full Size | PDF

Figures 2(b) and 2(c) show the generated linear and nonlinear currents, respectively. It can be seen that the induced nonlinear currents flow along the arms of the resonators in phase and their radiated fields interfere constructively, resulting in a 90° polarization rotation of the nonlinear far-field emission with respect to the pump polarization.

The DFG process of all the frequency components within the bandwidth of the illuminating femtosecond pulse is expected to lead to broadband THz generation [21]. To validate this effect, we employ Eq. (3) of the hydrodynamic model. Even though any two combinations of the frequency components of the incident pulse interact with each other in DFG, the number of possible pairs that can be calculated by our method is limited due to the frequency-domain approach, bounding the simulated pairs to a finite set. However, this computational scheme enables us to obtain a quick and accurate prediction of the emission. In addition, it allows us to optimize the efficiency at the level of a single generated frequency (${\omega }_{\mathrm{gen}}$) by simulating a finite number of combinations $({\omega }_1,{\omega }_2)$ that comply with ${\omega }_1-{\omega }_2={\omega }_{\mathrm{gen}}$, thus avoiding a time-consuming full wave simulation of the complete DFG spectrum.

In Fig. 3(a) we present the results of the model for THz generation by a DFG process on the metasurface presented in Fig. 2. The results are obtained for near-infrared short pulses with different pulse lengths. It is evident that the generated emission has a broad and tunable bandwidth in the THz range that is strongly dependent on the bandwidth of the source [Fig. 3(b)], exhibiting a shift in the spectral power peak as well. The obtained results of our computationally efficient frequency-domain model are in perfect agreement with results obtained in [20,21] that used a time domain approach.

 

Fig. 3. (a) Normalized power spectral amplitudes of the THz emission measured in the far field, generated by Gaussian pump pulses of various durations. (b) Dependence of the FWHM and peak intensity of the DFG delivered power of the metasurface versus the pulse duration of the pump pulse. The width of the generated pulses shows a nonlinear behavior that is proportional to $1/\tau$, as the time–bandwidth product states.

Download Full Size | PDF

This model allows us to study different types of complex NLMS-based elements for THz generation and control. In the following section, we present a new waveguide platform for generation, control, and guiding of the THz emission, consequently allowing the use of THz radiation in photonic and imaging devices [29–31]. Unlike previous designs, our proposed parallel plate waveguide (PPWG) is free of impedance mismatch, which is a crucial parameter in the design of electronic and photonic devices [32]. The suggested structure, as shown in Fig. 4(a), consists of a nonlinear metamaterial based photonic crystal (NLMPC) [17] with a periodicity $\mathrm{\Lambda }$. A metallic plate is positioned in distance $s$ across the emission plane, and a two-output port PPWG is formed. The generated THz radiation by the NLMPC is bounded between the deep-subwavelength nano-resonators that act as a mirror and the bottom plate.

 

Fig. 4. (a) Geometrical illustration of the THz NLMS activated PPWG, consisting of $N$ periods of flipping SRRs. (b) Radiated field from two inverted SRRs. The two different orientations of the meta-atoms are illustrated in the insets of (c) and (d), leading to the excitation of diverse modes. Theoretical predictions (dashed lines) are also displayed and match perfectly with the simulations, for a constant periodicity $\mathrm{\Lambda }=157\mathrm{\mu m}$. $\mathcal{T}$ is defined as the ratio of the normal Poynting vector $\mathcal{P}$ in the output and the generated THz power in the input. Both excitations share the same colormap.

Download Full Size | PDF

Exploiting the nonlinear phase inversion of the THz field when the orientation of the SRRs is flipped, as presented in Fig. 4(b), we can position the NLMPC in a way that guarantees quasi-phase-matching conditions of specific THz waveguide modes. The direction of the THz emission follows the momentum-matched Raman–Nath diffraction, as it was previously studied for the SHG [17]. The momentum conservation condition in DFG is

The ability to manipulate the geometrical orientation of the meta-atoms, allows us to excite either the TE or the TM modes simply by orienting the SRRs with their arms parallel or perpendicular to the THz propagation direction, respectively [see the inset in Figs. 4(c) and 4(d)].

The simulated waveguides have a finite number of periods set to $N=20$ and infinite width along the $y$- direction. The generated THz signal is measured in one of the two possible PPWG ports, due to the symmetric nature of the structure. Figures 4(c) and 4(d) show the normalized collected THz signal at the output ports versus frequency and PPWG plate separation for the TE and TM configurations, respectively. The dashed lines mark the waveguide modes calculated by Eq. (6). It is noticeable that the collected THz signal fits perfectly to the supported waveguide modes of the structures. The lower peaks that appear on the power ratio plots are a consequence of the finite size of the metasurface and show a $\mathrm{sinc}$-like behavior. Nonetheless, they cannot be counted as guided modes, as they display low transmitted energy levels of around 10% at maximum. An increase in the number of periods $N$ will affect the power ratio diagrams by sharpening the resonance and decreasing the amplitude of the lower local maxima.

In terms of the guided THz energy, our structure inherits the tunability of the NLMS generated power shown in Fig. 2. The transmitted energy, as shown in Fig. 5, can be approximated as the product of the power ratio and the power spectral envelope of the NLMSs. Altogether, it allows us to widely tune the spectral properties of the generated THz. In fact, this may provide a new degree of freedom in possible spectroscopy applications.

 

Fig. 5. Transmitted energy of TM excited modes for two different incident pump pulses. By tuning the shape of the pump, we can continuously control the emission profile into the modes of interest.

Download Full Size | PDF

In summary, we suggest a perturbative hydrodynamic model that describes the quadratic nonlinear optical interactions supported by the electron gas of plasmonic meta-atoms. This model allows us to calculate the generated nonlinear surface currents in DFG processes. The numerical simulations demonstrate that by addressing the problem in the frequency domain we can achieve accurate and high-speed results that can be easily implemented by commercial EM simulation tools. We specifically used the model to study the recently reported broadband THz generation process of metallic SRRs and show that it can serve as a novel method in the generation of tunable THz sources. In addition, a new waveguide platform was studied that generates the THz radiation directly into the guided modes. The waveguide supports both TE and TM mode excitation. The specific modes that are excited, and their spectral properties, can be controlled by the geometrical design of the NLMS, plate separation, and exciting pump pulse, proving a versatile new platform for THz generation and manipulation. The suggested structure can be easily embedded in existing designs due to its absence of impedance matching requirements. This new excitation scheme can also be applied in rectangular waveguide structures, leading to a whole new set of active THz waveguide platforms. Considering its unique characteristics, we believe that it may play a major role in forthcoming THz applications.