Investigation of Giant Nonlinearity in a Plasmonic Metasurface with Epsilon-Near-Zero Film

Abstract

Plasmonic metamaterials can exhibit a variety of physical optical properties that offer extraordinary nonlinear conversion efficiency for ultra-compact nanodevice applications. Furthermore, the optical-rectification effect from the plasmonic nonlinear metasurfaces (NLMSs) can be used as a compact source of deep-subwavelength thickness to radiate broadband terahertz (THz) signals. Meanwhile, a novel dual-mode metasurface consisting of a split-ring resonator (SRR) array and an epsilon-near-zero (ENZ) layer was presented to boost the THz conversion efficiency further. In this paper, to explore the mechanism of THz generation from plasmonic NLMSs, the Maxwell-hydrodynamic multiphysics model is adopted to investigate complex linear and intrinsic nonlinear dynamics in plasmonics. We solve the multiphysics model using the finite-difference time-domain (FDTD) method, and the numerical results demonstrate the physical mechanism of the THz generation processes which cannot be observed in our previous experiments directly. The proposed method reveals a new approach for developing new types of high-conversion-efficiency nonlinear nanodevices.

1. Introduction

Electromagnetic metamaterials with unique geometrical shapes and sizes provide some physical properties that cannot be produced by natural materials. The applications of metamaterials are numerous, such as in phone antennas and nonlinear optical devices. Nowadays, there are some typical metamaterials including left-handed materials and photonic crystals. The limitations of these two materials are their complex structures and restricted bandwidth. However, plasmonic materials can handle these limitations, and thus have attracted tremendous attention in regard to nonlinear optics and microwaves. Surface plasmon resonance (SPR) occurs when electrons on the metal surface are excited by light at a specific incident angle [1]. SPR is the key feature dominating the linear and nonlinear optical responses of the plasmonic metasurfaces based on the metallic nanostructures. Various photonic functions based on the nonlinear effects of plasmonic metasurfaces have flourished in recent years, such as optical sensing, ultrashort pulse generation, nanoantennas, and optical signal processing [2,3,4,5]. Because of the poor conversion efficiency of higher-order nonlinearities, nonlinear interactions in plasmonics are weaker than linear interactions and depend on the field amplitude of excitation. In experiments, high-intensity excitation is ideally suited for plasmonic metasurfaces to achieve strong field enhancement. Simultaneously, they can be improved by focusing light using resonant modes, and it is possible to modify and modulate the nonlinear nature of metamaterials by designing the geometry of the plasmonic structure. This leads to exciting and useful effects including second-harmonic generation, third-harmonic generation, and other high-harmonic generations [6].

Although plasmonic metasurfaces based on metallic nanostructures can be viewed as a general concept, the sources of nonlinearity of each plasmonic metasurface based on metallic nanostructures are quite diverse. Recently, ENZ films, including indium-tin oxide (ITO) and aluminum-doped zinc oxide (AZO) films, have emerged as promising options for enhancing nonlinear effects due to their ability to greatly enhance the electric field at the ENZ interface [7,8]. However, the overall photoelectric performance of AZO films is worse than that of ITO thin films. ITO films have nearly flat dispersion in ENZ mode, so they exhibit strong optical nonlinearities. No matter what kind of plasmonic metasurface based on a metallic nanostructure is considered, it is of great significance to study the optical response between the electromagnetic field and metasurface to reveal its working principle and design a working device.

Meanwhile, THz plays an increasingly critical role in a wide range of applications, including imaging, sensing, and communications. With the development of nanoscale fabrication techniques, devices for THz generation and operation now have more sophisticated structures. In the past decades, various generation methods have been tried, e.g., nonlinear organic crystals, topological materials, etc. [9]. However, these devices are spatially separated and the organic crystals are limited by quasi-phase-matching conditions. The advent of nonlinear metasurfaces will allow for the efficient generation and control of THz waves that are integrated into a single device. Researchers have demonstrated the possibility of using plasmonic metasurfaces as broadband THz sources. The combination of an ENZ material and metasurface can be guided to generate nonlinear waves through reasonable symmetry and geometry, creating a broadband THz signal powered by the difference-frequency generation (DFG) effect and freeing them from the quasi-phase-matching conditions [10]. This provides a promising solution for simultaneous THz emergence and manipulation, which will lead to lower costs and more efficient THz systems and drive the development of current THz fabrication.

In this paper, we investigate the use of an SRR array on an ultra-thin ITO film device to emit THz signals with high efficiency. In order to reveal the mechanism of THz generation from the metasurface, we utilize the Maxwell-hydrodynamic model to study the complex linear and nonlinear optical effects of electronic gases in metals and an ITO film [11,12,13]. Hybridization of the metasurface with an ITO film can help lower the threshold for nonlinear interactions by enhancing the strong near-field distribution of resonant meta-atoms at the desired frequency. The plasmonic resonant modes of the SRRs and the ENZ modes of the ITO film are spatially overlapped to achieve strong mode coupling, thus improving the nonlinear effect of THz generation in ITO. With the help of a parallel FDTD method, this paper obtains the numerical solution for a Maxwell-hydrodynamic multiphysics model for nanodevices with any geometry or parameters. Based on our simulation method, we numerically reveal the mechanism of significant THz enhancement resulting from the ENZ mode of the ITO film, which is essential for designing nonlinear effect devices.

2. Method and Dual-Mode Metasurface

2.1. Multiphysics Hydrodynamic Model

The interaction between the electromagnetic fields $\mathbf{E}$ and $\mathbf{H}$ and the nonmagnetic materials can by described by Maxwell’s equations:

$$\frac{1}{\mu }\nabla \times \mathbf{B}=\epsilon \frac{\partial \mathbf{E}}{\partial t}+\mathbf{J}$$

(1)

$$\frac{\partial \mathbf{B}}{\partial t}=-\nabla \times \mathbf{E}$$

(2)

where $\epsilon$ and $\mu$ are the material permittivity and magnetic permeability and $\mathbf{J}$ is the current density, i.e., the linear and nonlinear polarization currents generated by electromagnetic waves acting on free electron gases in metallic materials. The hydrodynamic equations describe the nonlinear response characteristics of the electromagnetic wave acting on the free electron gas within the metallic metamaterial [14,15,16,17].

$$\frac{\partial {\mathbf{v}}_e}{\partial t}+({\mathbf{v}}_e·\nabla ){\mathbf{v}}_e=-\frac{e}{m}(\mathbf{E}+{\mathbf{v}}_e\times \mathbf{B})-\gamma {\mathbf{v}}_e$$

(3)

$$\frac{\partial n_e}{\partial t}+\nabla ·(n_e{\mathbf{v}}_e)=0$$

(4)

where m and e represent the electron mass and charge, respectively. $\gamma$ is the phenomenological damping frequency constant (used to measure optical losses). $n_e$ and ${\mathbf{v}}_e$ describe the electron density and velocity. ($\mathbf{E}+{\mathbf{v}}_e\times \mathbf{B}$) is the Lorentz force. The basis for the physical origin of nonlinear and nonlocal effects is ${\mathbf{v}}_e·\nabla {\mathbf{v}}_e$, ${\mathbf{v}}_e\times \mathbf{B}$, and $n_e{\mathbf{v}}_e$. Equation (4) is the current continuity equation. In order to obtain the Lorentz force acting on the electron gas, it is necessary to couple the hydrodynamic equations with the macroscopic electromagnetic field determined via Maxwell’s equations, which can be linked by the macroscopic current density and charge density terms.

$$\mathbf{J}=-e{n}_e{\mathbf{v}}_e$$

(5)

$$\rho =e(n_e-n_0)$$

(6)

To ensure charge neutrality, the initial electron density $n_e$ is set to be equal to the positive ion density $n_0$ without excitation. When substituting Equations (5) and (6) into Equations (3) and (4), we can obtain:

$$\frac{\partial \mathbf{J}}{\partial t}={\epsilon }_{\infty }{\omega }_p^2\mathbf{E}-\gamma \mathbf{J}+\frac{e}{m}(\rho \mathbf{E}-\mathbf{J}\times \mathbf{B})+\nabla ·(\frac{1}{\rho +{\epsilon }_{\infty }m{\omega }_p^2/e}\mathbf{J}\mathbf{J})$$

(7)

$$\frac{\partial \rho }{\partial t}-\nabla ·\mathbf{J}=0$$

(8)

where ${\omega }_p=\sqrt{e^2{n}_e/{\epsilon }_{\infty }m}$ is the plasma frequency. The Equations (7) and (8) above equation can be used to simplify by substituting $\rho ={\epsilon }_{\infty }\nabla ·\mathbf{E}$ which is transformed into Equation (9).

$$\frac{\partial \mathbf{J}}{\partial t}={\epsilon }_{\infty }{\omega }_p^2\mathbf{E}-\gamma J+\frac{e}{m}({\epsilon }_{\infty }(\nabla ·\mathbf{E})\mathbf{E}-\mathbf{J}\times \mathbf{B})+\nabla ·(\frac{1}{{\epsilon }_{\infty }(\nabla ·\mathbf{E})+{\epsilon }_{\infty }m{\omega }_p^2/e}\mathbf{J}\mathbf{J})$$

(9)

Equations (1), (2) and (9) constitute a self-calculating set of equations for solving the response of free electrons to an applied electromagnetic field. The numerical solution of this set of multiphysics field equations using the FDTD method provides a description of the complex electron motion in metals and enables the study of the nonlinear behavior of the electron gases.

2.2. Numerical Approach

The nonlinear optical response of free conduction electrons can be described by the Maxwell-hydrodynamic model for analyzing the electromagnetic fields of the SRR array on an ITO film device. If the metamaterial is described by the effective medium theory, its constitutive parameters are continuous, and the material parameters will be very complicated. However, a metamaterial can be directly calculated using the FDTD method as a medium with electromagnetic properties. The FDTD scheme uses Yee cells to store the electric and magnetic fields in space and can directly simulate the scattering and time-domain propagation response of electromagnetic field signals [18]. In addition, the FDTD method solves various relatively complex electromagnetic problems with simple and high accuracy and can be combined with various parallel computing techniques to improve computational efficiency [19,20,21].

2.3. Computational Grid

To solve Equations (1), (2), and (9), we used the Yee scheme to arrange all components of $\mathbf{E}$ and $\mathbf{H}$ with a special interleaving, as shown in Figure 1. The electric field $\mathbf{E}$ is located at the center of the surface and defined at l + 1/2 time steps, and the magnetic field $\mathbf{H}$ and current density $\mathbf{J}$ are located at the edges of the grid and at the center of the surface, respectively, and defined at l time steps. A spatial and temporal sampling of electric and magnetic fields is performed alternately. This grid arrangement conforms to the Maxwell-hydrodynamics equations, which can capture the coupling effect between the electromagnetic field and the free electrons in the metal metasurface.

2.4. Final Discretized Formula

The auxiliary differential equation (ADE) technique combined with the FDTD method can extend the electric and magnetic fields, as well as charge and current densities, into first-order linear and second-order nonlinear responses [22]. This can be applied to complex dispersive dielectric material models. The ADE method solves linear effects explicitly and also obtains intermediate values that are used to update the nonlinear response expressly. That is, a two-step splitting method for solving Maxwell’s equations is used. Firstly, for the decomposition of Equation (9):

$$\frac{\partial \mathbf{J}}{\partial t}={\epsilon }_{\infty }{\omega }_p^2\mathbf{E}-\gamma \mathbf{J}$$

(10)

$$\frac{\partial \mathbf{J}}{\partial t}=\frac{e}{m}\left[{\epsilon }_{\infty }(\nabla ·\mathbf{E})\mathbf{E}-\mathbf{J}\times \mathbf{B}\right]+\nabla ·(\frac{1}{{\epsilon }_{\infty }(\nabla ·\mathbf{E})+{\epsilon }_{\infty }m{\omega }_p^2/e}\mathbf{J}\mathbf{J})$$

(11)

where the ${\epsilon }_{\infty }\nabla ·\mathbf{E}$ term represents the charge density. A linear intermediate current density term is obtained by central differencing Equation (10), while the nonlinear response is obtained by Equation (11). The intermediate current density term ${\mathbf{J}}^{(1)}$ at time step l + 1 is:

$${\mathbf{J}}^{(1),l+1}=\frac{1-0.5\Delta t\gamma }{1+0.5\Delta t\gamma }{\mathbf{J}}^l+\frac{{\epsilon }_0{\omega }_p^2\Delta t}{1+0.5\Delta t\gamma }{\mathbf{E}}^{l+1/2}$$

(12)

In view of the differences between $\mathbf{J}$, $\mathbf{E}$, and $\mathbf{B}$ in terms of their spatial and temporal locations, the interpolation method is used to approximate the values at the same points in space and time as $\mathbf{J}$. Take the example of the updated x-component of $\mathbf{J}$, which is located in the same center of the grid as $\mathbf{E}$ in space and is updated at the same time as $\mathbf{B}$ in time. Thus, the terms $\nabla ·\mathbf{E}$ and $\mathbf{B}$ are calculated by the following equation:

$$\begin{array}{l}\overline{\nabla ·{\mathbf{E}}_{i,j+1/2,k+1/2}^{t+1/2}}=\frac{1}{2}(\nabla ·{\mathbf{E}}_{i+1/2,j+1/2,k+1/2}^{t+1/2}+\nabla ·{\mathbf{E}}_{i-1/2,j+1/2,k+1/2}^{t+1/2})\\ \overline{B_{yi,j+1/2,k+1/2}^{t+1/2}}=\frac{1}{4}(B_{yi,j+1/2,k}^{t+1}+B_{yi,j+1/2,k+1}^{t+1}+B_{yi,j+1/2,k}^t+B_{yi,j+1/2,k+1}^t)\\ \overline{B_{zi,j+1/2,k+1/2}^{t+1/2}}=\frac{1}{4}(B_{zi,j,k+1/2}^{t+1}+B_{zi,j+1,k+1/2}^{t+1}+B_{zi,j,k+1/2}^t+B_{zi,j+1,k+1/2}^t)\end{array}$$

(13)

The results of the above steps allow for the explicit solution of Equation (11) to obtain the nonlinear intermediate current density term. We update the electric field in the next step based on the final current density calculated from the two intermediate terms obtained in the previous steps.

$${\overline{J}}_{xi,j+1/2,k+1/2}^{l+1/2}=\frac{1}{2}(J_{xi,j+1/2,k+1/2}^l+J_{xi,j+1/2,k+1/2}^{(1),l+1})$$

(14)

By using two-step splitting, the linear and nonlinear responses of the metallic metasurface are split, preserving their dominant nonlinear effects and time-domain expressions while greatly simplifying the calculation of the Maxwell’s hydrodynamic model and improving the computational efficiency by a significant amount compared to the implicit difference method.

2.5. GPU Programming Structure

With the development of nanoscale manufacturing technology, metal metamaterials have smaller and finer structures. The nonlinear optical behaviors and high-order harmonic signals can be generated by the high-energy incidence wave to stimulate the free electron gases in metamaterials. However, this requires large iteration steps and long iteration times to ensure the effective attenuation of pulses. Considering that this takes many computing resources, the FDTD method can be combined with parallel techniques, where parallel processing helps to speed up the computation [23,24,25].

Based on the implementation of the graphic processing unit (GPU) platform multiphysics model and the FDTD method, firstly, all data are initialized on the CPU (host) and then loaded into the GPU’s (device’s) global memory. Next, the EM field’s updated equations for each point are implemented on a set of parallel Compute Unified Device Architecture (CUDA) threads, and finally, the completed calculation is passed back to the CPU for post-processing. The whole computational process is shown in Figure 2. By implementing the FDTD code on the GPU, the computational efficiency can be significantly improved.

2.6. Dual-Mode Plasmonic Metasurface with ENZ Material

We explore the mechanism of the THz emission of the hybrid nonlinear metasurfaces by simulating a plasmonic with the proposed two-step split FDTD algorithm. The hybrid nonlinear metasurfaces consist of three parts: a glass substrate as the overall support at the bottom, an ITO film as the middle layer, and a periodically arranged SRR array as the top layer. Our simulations involve analytical analysis of individual cells to obtain the features and properties of the entire structure. The 23 nm thick ITO metasurface is fabricated between a glass substrate and a 40 nm thick gold SRR structure. The relative permittivity of glass is 2.25. The detailed geometries of the elements are also illustrated in Figure 3.

ENZ materials with small permittivity exhibit highly efficient nonlinear optical phenomena. ITO, characterized by excellent optoelectronic properties, is one of them. ITO can rapidly produce an enhanced nonlinear optical response in the permittivity near-zero spectral region. Simultaneously, modulating carrier concentrations in ITO films allows for tuning of the ENZ mode wavelength, opening up an even wider range of applications. The significant nonlinearity of ITO films is attributed to the strong field constraints associated with the ENZ mode, which can only be supported if the film is thin enough to prop up the ENZ mode, and the electric field distribution of the ENZ mode is perpendicular to the film surface; thus, the ENZ mode cannot be excited without any near-field coupling under normal illumination. THz is produced by making use of plasmonic resonant meta-atoms in this paper, enhancing the second-harmonic signal and thus increasing the nonlinear response of the ITO. The ENZ wavelength, as measured by spectroscopic ellipsometry, is approximately 1432 nm.

SRR metamaterials can generate broadband THz waves by optical rectification, while the magnetic dipole resonance of SRR excites THz emission [26,27,28]. The combination of plasmon resonances and ENZ modes in hybrid metasurfaces can amplify the field-enhanced effects and enhance the nonlinear response of the ITO material, resulting in efficient THz generation.

3. Results

3.1. Simulation Setup

The dual-mode metasurface, which consisted of a gold nano-SRR structure supported on a glass substrate coated with an ITO film, was utilized as an experimental test, while the SRR metasurface, which had the same configuration but without the ITO film, was used as a control array. The simulated metasurface unit is shown in Figure 4. The metallic material of SRR was made of gold, the periodic boundary conditions were in the x- and y-directions, and the total-field and scattered-field (TF/SF) technique was deployed to take the normal incident plane wave as the time-domain excitation signal in the z-direction of the cell structure, thus allowing the overall properties to be easily obtained with as few computational resources as possible by calculating the structure of the cell. Moreover, perfectly matched layers (PMLs) were applied in the z-direction to absorb broadband outgoing waves, truncating the computational domain from infinity to a finite size. In order to ensure convergence of the algorithm, the spatial step was $\Delta x=\Delta y=\Delta z=1\times {10}^{-9} \mathrm{m}$, and the time step was $\Delta t=1.5\times {10}^{-18} \mathrm{s}$.

It has been explained that the Drude model $\epsilon (\omega )={\epsilon }_{\infty }-{\omega }_p^2/({\omega }^2+i\omega \gamma )$ can be used to describe the dispersive property of metals, semiconductors, and plasmonics. ${\epsilon }_{\infty }$ is the infinite frequency permittivity, ${\omega }_p=\sqrt{n{e}^2/({\epsilon }_{\infty }m)}$ is the plasma frequency that determines the carrier density, and $\gamma$ is the phenomenological damping frequency constant. The ITO parameters were defined to be: ${\epsilon }_{\infty }=4.1$, ${\omega }_p=2.7297\times {10}^{15}\mathrm{rad}/\mathrm{s}$, and $\gamma =10.68\times {10}^{13}\mathrm{rad}/\mathrm{s}$. Thus, the ENZ wavelength of ITO referring to the zeros of the real permittivity ${\epsilon }^{\prime }(\omega )=0$ was approximately given by ${\lambda }_{ENZ}=2\pi c/\sqrt{{\omega }_p^2/{\epsilon }_{\infty }-{\gamma }^2}=1340\mathrm{nm}$. Thus, the parameters of the gold and ITO used in the hydrodynamic model are listed in

Table 1. Parameter list of gold and ITO used in the model.

	Background Permittivity ${\mathbf{\epsilon }}_{\mathbf{\infty }}$	Electron Density $\mathit{n}$	Collision Frequency $\mathbf{\gamma }$
Au	1	5.8613 × 1028	1.2 × 1014
ITO	4.1	2.3406 × 1027	2.9326 × 1014

.

The measured dual-mode and SRR metasurfaces with polarization-dependent transmission spectra are depicted in Figure 5. The resonances caused by the linear response of the dual-mode metasurface are exhibited at 1375 and 1750 nm. Resonance is surveyed in the transmission spectrum at 1350 nm, which is close to the ENZ wavelength of the ITO, that is, when the ITO film is added to the same excitation, this produces a transmission dip as the ENZ layer absorbs the scattered Ez field.

3.2. Revealing the Enhanced Nonlinearity from the Dual-Mode Metasurface

Practically, a sufficiently strong second harmonic is radiated to generate a THz signal when the polarization direction of the incident light satisfies the plasmonic excitation condition of the bare ITO film. Simultaneously, the resonant meta-atoms significantly reinforce the fundamental electric field component perpendicular to the ENZ film, which contributes to the THz signal generated from the ITO interlayer at a normal incidence. This behavior can be observed in Figure 6, where the evolution of the THz peak-to-peak amplitude is shown as the central wavelength of the pump laser changes on the SRR and dual-mode metasurfaces. Pumped at a wavelength of 1350 nm, which approaches the ENZ state of ITO, the SRR metasurface exhibits a THz generation amplitude that is ~90 times smaller, while at the same time, strongly enhanced THz generation occurs on the dual-mode metasurface. The THz enhancement has fundamental wavelength dependence, and when the experimental results are compared with the simulation results, they show a reasonable, good agreement with the experimental data.

We conducted simulations of THz emission spectra from the dual-mode metasurfaces without accounting for the ITO nonlinearity in order to illustrate the physical mechanism by which the ITO layer increases the THz radiation intensity of dual-mode metasurfaces. Considering only the nonlinearity of the Au SRR, the THz signal dropped by approximately two orders of magnitude, indicating the significant contribution of the ITO layer to the enhancement of THz emission.

The strong coupling between the magnetic dipole mode of the SRR, the fundamental mode of the gold antenna, and the bulk plasmon mode of the ITO layer plays a significant role in the THz generation by the dual-mode metasurface. In order to investigate the role of mode coupling, a silica spacer was sandwiched between the SRR and ITO films to tune the coupling. Figure 7 shows the linear transmission spectra of silica spacers of different thicknesses set in the dual-mode metasurface at a normal incidence. As the thickness of the silica spacer increases, the coupling between the SRR resonance and the ENZ mode weakens, resulting in a decrease in the fundamental electric field polarized along the z-direction and a sharp decrease in the THz signal, as depicted in Figure 7b.

4. Conclusions

We investigated the physical mechanism of electromagnetic radiation enhancement in ENZ materials by combining ITO substrates with metallic metamaterials. The Maxwell-hydrodynamic model was utilized to study the nonlinear response in metals and an ITO film under the excitation of incident electromagnetic fields. By comparing the THz electric field intensities of SRRs with and without ITO substrates, we demonstrated that the ENZ mode can be excited in the presence of the ITO substrate. The strong coupling between the ENZ mode and the plasmonic metasurface allows for a strong nonlinear response over a wide wavelength range, thereby validating our multiphysics model. In addition, ITO substrates can also be combined with other metallic nanostructures to capture the nonlinear response of metallic metamaterials and other quantum effects, which are expected to have applications in the design of nonlinear metamaterial devices.