The platform of LNOI is opening new avenues to optical frequency conversion for its strong second-order nonlinear susceptibility tensor. This platform is formed from nanometer-thick LiNbO₃ films bonded on SiO₂ grown on an insulator (silica) substrate using a smart-cut technique, which results in devices with much increased index contrast (Δn > 0.6) and reduced modal size. The main obstacles are the challenges of the growth of thin crystalline films on low refractive index substrates and the subsequent nanostructuring. Figure 1A shows the schematic representation of the LiNbO₃ metasurface. In order to utilize the largest second-order nonlinear susceptibility d₃₃ maximizing THz generation, an x-cut LNOI wafer of 300 nm thickness bonded on a 4-µm thick SiO₂ is fabricated on a 500-µm thick Si buffer layer (NanoLN, Jinan Jingzheng Electronics Co., Ltd.). We refer to the x-axis of the LiNbO₃ principal crystallographic coordinate and the polarization of laser as the x-axis and z-axis of the experimental coordinate system, respectively. The optical axis of LiNbO₃ lies in the y-z plane and can rotate around the x-axis. Figure 1B shows a scanning electron microscopy image of the realized metasurface with a period constant of 950 nm. The magnified image shows the well-fabricated nanoholes of the LiNbO₃ metasurface with a diameter of 550 nm. The pores are directly drilled by electron-beam lithography.

Periodic pores in the LiNbO₃ film modulate the refractive index of the metasurface, resulting in a guided-mode resonance [45–47] which can couple the normally incident light into the metasurface. The guided-mode resonance is very sensitive to the geometric parameters and the refractive indexes of structured metasurfaces. To design optimized geometry, numerical simulations of the structure are carried out by commercial software (COMSOL Multiphysics) by means of the finite element method. When the collimated light incidents normally onto the LiNbO₃ metasurface, sharp peaks in the simulated reflection spectra are observed, as shown in Figure 1C. When the polarization of the excitation is parallel to the optical axis of LiNbO₃ (extraordinary electric field, E_e), the resonant wavelength is 1,375 nm with a bandwidth of 5 nm. The resonant wavelength shifts to 1,395 nm when the polarization of excitation is perpendicular to the optical axis (ordinary electric field, E_o). The anisotropic refractive indexes of x-cut LiNbO₃ result in the resonant wavelength varying with excitation polarization. The extraordinary refractive index of 2.14 for LiNbO₃ near the resonant wavelength is less than the ordinary refractive index of 2.21, which results in a shorter resonant wavelength [48].

By inducing strong electromagnetic resonances, the electric field intensity is locally enhanced, as compared to unstructured materials. THz waves are mainly contributed by the electric field oriented along the optic axis in LiNbO₃ due to the overwhelmingly dominant d₃₃ coefficient. Figure 1D shows electric fields in the y-z (left) and x-z (right) sections, respectively. The amplitude of electric fields in the bottom of Figure 1D shows a dramatic enhancement in the metasurface compared to the unstructured film. The arrows at the top of Figure 1D show the direction of electric fields and indicate weak components of E_y and E_x . In contrast to the plasmonic enhancement base on the metal metasurface [33, 34], the enhanced electric field extends over a larger region around the pores in the LiNbO₃ metasurface, which will, in general, boost the efficiency of THz generation.

Figure 2 schematically illustrates the experimental setup and key elements. The experimental setup is based on a Ti: sapphire amplifier (Coherent, Libra) operating at 800 nm with 4 mJ pulse energy, 100 fs pulse duration, and 1 kHz repetition rate. The laser pulses are split into two beams. The main portion (80%) of the pulse energy is used to pump an optical parametric amplifier (OPA) to generate near-infrared (NIR) light with a pulse duration of approximately 140 fs and a central wavelength tunable from 1,100 nm to 1,610 nm. NIR pulses from the OPA incident normally onto the 300-nm thick LN metasurface with a beam diameter of 5 mm. The generated THz waves are gathered using three 90° off-axis parabolic mirrors and focused on a 1-mm thick ZnTe wafer. The other portion (20%) of pump energy serves as the probe beam for electro-optical (EO) detection. A motorized translation stage is used to change the time delay between THz and probe pulses. A wire grid polarizer is inserted in the collimated THz path to analyze the polarization of the generated THz waves. The experimental coordinate system is defined by the pump pulses. The reverse direction of pump pulses is the x-axis. The polarization direction of pump pulses is fixed to be vertical and regarded as the z-axis. The optical axis of the x-cut LiNbO₃ metasurface lies in the y-z plane, which is rotated to vary the angle φ between the optical axis and y-axis in the experiment. When the optical axis is vertical and parallel to the polarization of the excitation, φ = 90°. Although the optical axis is horizontal and perpendicular to the polarization, φ = 0°.

Figure 3 shows the enhanced THz generation from the LiNbO₃ metasurface. In simulation, the Fano-type resonant wavelength of the LiNbO₃ metasurface is located at 1,375 nm and 1,395 nm, when the optical axis is parallel and perpendicular to the polarization of the pump pulses, respectively. In the experiment, the THz electric field depending on pump wavelength indicates the same result. The THz polarizer ensures that THz polarization is parallel to the pump polarization.

Figure 3A shows the anisotropy of the resonant enhancement of the THz electric field from the LiNbO₃ metasurface, and the pump intensity is 6 GW/cm². When the optical axis of the LiNbO₃ metasurface is rotated in the vertical direction (φ = 90° in Figure 2) and the wavelength is tuned to 1,375 nm, the extraordinary THz waves (E_e) are recorded as the red curve. The detected THz electric field is near 0.4 kV/cm. The conversion efficiency of the THz electric field is estimated to be 5 × 10^–9. When the optical axis is rotated to the horizontal direction (φ = 0° in Figure 2) and the wavelength is tuned to 1,395 nm, the ordinary (E_o) THz waves are recorded as the blue curve. Figure 3A shows that the extraordinary THz electric field is eight times larger than the ordinary counterpart because d₃₃ dominates the second-order nonlinear process [49].

Figure 3B shows the wavelength selectivity of the resonant enhancement for THz generation from the LiNbO₃ metasurface. In the experiment, the wavelength of pump pulses is tuned to search for the resonance wavelength when the optical axis of the LiNbO₃ metasurface is the z- and y-axis in Figure 2. The THz peak-to-peak electric field for extraordinary (φ = 90°) and ordinary (φ = 0°) waves are recorded as red and blue balls in Figure 3B by tuning the wavelength. The central resonant wavelengths of the LiNbO₃ metasurface are 1,375 nm and 1395 nm for E_e and E_o, respectively. This exactly agrees well with our simulation of the designed LiNbO₃ metasurface. At zero detuning, the largest fraction of the excitation pulse spectrum—which is much wider than the width of the resonance—is resonantly coupled to the metasurface, and the THz electric field is strongest. When detuning is slight, the THz electric field decreases monotonically with increasing detuning of the excitation pulse. This is induced by decreasing spectral power overlap between the laser spectrum and the metasurface resonance. This confirms that the THz enhancement is indeed contributed by the resonance-induced electric field enhancement on a metasurface. The full width at half maximum of the THz electric field is about 50 nm in Figure 3B, which is much wider than 5 nm in the simulated reflection spectra, as shown in Figure 1C. It is speculated that the fabrication error and the roughness of the LiNbO₃ metasurface cause the discrepancy.

Figure 3C demonstrates the resonant enhancement of THz generation in the LiNbO₃ metasurface compared with the unstructured LiNbO₃ nanofilm of the same thickness. THz waves are generated at an incident excitation intensity of 6 GW/cm². When the optical axis of the LiNbO ₃ metasurface is along the pump polarization and the excitation wavelength is tuned at the resonant wavelength of 1,375 nm, the extraordinary THz electric field (red) is enhanced by a factor of 10 compared to the unstructured LiNbO₃ nanofilm. The ordinary THz electric field (blue) from the resonant metasurface is also ten times larger. The normalized Fourier contents of the THz electric field are the same as shown in the inset. The metasurface only enhances the THz electric field and does not modify the spectral distribution. This demonstrates the THz bandwidth of the metasurface up to 3 THz that is limited mostly by the excitation pulse of duration of ∼140 fs (spectral width ∼13 meV or 3.2 THz) and by the Reststrahlen region of the ZnTe wafer. The enhancement of the THz electric field is one order of magnitude lower than that in the simulation. The spectral width of the pump pulse spectrum (12 nm) is much wider than the width of the resonance (5 nm). The poor overlap of spectral power between the laser spectrum and the metasurface resonance partially explains the obvious discrepancy.

Figure 3D shows the THz electric field as a function of the incident pump power. The pump polarization is parallel to the optical axis of the metasurface. As illustrated in Figure 3D, the THz electric field is linearly proportional to the incident pump power. THz waves from 0.5-mm thick LiNbO₃ or 2-mm thick ZnTe crystal wafers are believed to originate from optical rectification. The enhanced THz electric field originates from the resonant field enhancement in the LiNbO₃ metasurface compared with that in the unstructured LiNbO₃ nanofilm. In the experiment, we used the femtosecond pulses with plane wavefront to excite resonantly the metasurface, which limited the available pump intensity. The THz generation from the metasurface as a function of the excitation intensity is still perturbative and not saturated in the experiment. The measured THz electric field is the same order of magnitude as the ZnTe(110) emitter pumped with the same wavelength and intensity and much higher than plasmon-based metallic nanostructures [33, 34].

To further investigate the contribution of different nonlinear coefficients to THz generation, the THz peak-to-peak field versus the azimuthal angle φ is measured by rotating the optical axis of different LiNbO₃ emitters. In the experiment, the THz polarizer can analyze polarization-resolved THz waves, which are excited by polarization along the z-axis. Figure 4 shows the measured (dots) and fitted (lines) THz peak-to-peak fields for a 0.5-mm thick LiNbO₃ crystal wafer and LiNbO₃ metasurface excited with the corresponding resonance wavelength. The first row of Figure 4 is the E z T H z   component of the THz electric field where the THz polarization is parallel to the pump polarization. The middle row of Figure 4 is the E y T H z component of the THz electric field where the polarization of THz waves is perpendicular to that of excitation. In the polar plots, the azimuthal angles φ = 90 ° correspond to the optical axis along the z-axis (parallel to pump polarization), and φ = 0 ° correspond to the optical axis along the y-axis (perpendicular to pump polarization).

For the 0.5-mm thick LiNbO₃ wafer, Figure 4A shows the E z T H z component of the THz electric field when THz polarization is parallel to pump polarization. It clearly shows that E z T H z starting from a maximum at φ = 90 ° crosses zero at φ = 180 ° and returns back to a maximum at φ = 270 ° , with a π phase shift from the origin. This eight-like shape is due to the dominating d₃₃ coefficient, which couples z-polarized excitation fields with THz waves. This has been reported in second-harmonic generation from non-resonant LiNbO₃ thin films [49]. The influence of other coefficients can be seen in the E y T H z component of the THz electric field shown in Figure 4B. Here, the maximal E y T H z component is about one half of the maximal E z T H z and is achieved in an azimuthal angle which is approximately φ = 45 °   and φ = 225 ° in the y-z plane. This indicates that E y T H z is generated mostly by the d₃₁ coefficient, which combines z-polarized and y-polarized excitation fields. However, there is an asymmetry in the pattern, which is not exactly along the 45° direction. This can be attributed to the action of the d₂₂ coefficient coupling y-polarized excitation with y-polarized THz waves. This contribution is interfering constructively or destructively with the lobes originating from d₃₁. Comparing the THz electric field with orthogonal polarization in Figures 4A,B, the linearly polarized THz waves are validated even if the polarization plane rotates relative to the orientation of the optical axis. The phase-resolved THz waves allow us to unambiguously determine the second-order nonlinear susceptibility of the LiNbO₃ crystal by fitting the experimental data in Figures 4A,B. If we normalized the largest component d₃₃ of effective second-order nonlinear susceptibility as ‒1, the other two small components d₃₁ and d₂₂ are -0.1879 and 0.0969, respectively [49]. These data agree well with those obtained from second-harmonic generation. Figure 4C shows the total THz electric field. The azimuthal dependence is similar with the eight-shape of E z T H z . However, the THz electric field is not zero at φ = 0 °   due to the considerable contribution from E y T H z .

The LiNbO₃ metasurface excited in the resonant wavelength of 1,375 nm shows a similar nonlinear behavior as for the LiNbO₃ wafer, as discussed before. However, Figure 4E shows the different azimuthal dependence on the orientation of the optical axis for the E y T H z component of the THz electric field. The direction of the maximal E y T H z rotates about 90 ° , and there is an isotropic nonzero component compared with Figure 4B. The aforementioned two differences are attributed to the resonant enhancement E_x of the excitation field in the LiNbO₃ metasurface. As shown in Figure 1D, there exist considerable E_x and E_y components when the LiNbO₃ metasurface is resonantly excited with the polarization along the z-axis. Although the THz electric field in Figure 4E is generated mostly by d₃₁, d₂₂ is interfering constructively or destructively with the lobes. The relative values of E_x and E_y of the excitation field in the LiNbO₃ metasurface determine the amplitude and phase of the generated THz electric field E y T H z . There is no E_x of the excitation field in the LiNbO₃ wafer, while E_x is larger than E_y in the resonantly excited LiNbO₃ metasurface. This introduces a π phase shift of E y T H z in the LiNbO₃ metasurface compared with the wafer. The interference with E y T H z originated from the d₃₁ coefficient rotates in the maximum direction of E y T H z to about φ = 135 ° . The E_x of the excitation field in the LiNbO₃ metasurface is independent on the orientation of the optical axis, which contributes isotropic E y T H z , as shown in Figure 4E. The eight-like shape in Figure 4F is more stretched in the direction φ = 90 ° and φ = 270 ° , indicating a stronger nonlinear response for these polarization angles than in the case of the wafer. This is also due to the resonant electric field enhancement, which increases the THz electric field in this direction but is absent in the orthogonal polarization direction owing to the birefringence.

The LiNbO₃ metasurface excited in the resonant wavelength of 1,395 nm shows two eight-like patterns for E z T H z and E y T H z of the THz electric field. As discussed before, E z T H z is dominantly determined by the d₃₃ coefficient. Due to the resonant enhancement of the E_y excitation field in the LiNbO₃ metasurface, E y T H z is dominantly determined by the d₂₂ coefficient, which also produces the eight-like pattern, as shown in Figure 4G. The d₃₃ coefficient is ten times larger than the d₂₂ coefficient. By coincidence, the THz electric field is also ten times enhanced when the excitation wavelength is 1,395 nm, and its polarization is perpendicular to the optical axis. These two facts produce the nearly isotropic pattern of THz waves by coincidence, as shown in Figure 4I.

Metasurfaces composed of ensembles of dielectric nanoresonators with subwavelength dimensions have shown great flexibility in achieving versatile functions in THz technologies. The coexistence of resonant nonlinear and local field enhancement in properly designed metasurfaces opens up an alternative route for high-efficiency nonlinear devices. The geometry of the individual meta-atoms along with their crystalline orientation and lattice structure enables one to control the polarization state of the emitted light based on that of the excitation light. In order to explore the polarization of THz waves generated from the LiNbO₃ metasurface, Figure 5 shows the dependence of THz polarization on the orientation of the optical axis. By means of polarization-resolved THz detection, Figure 5A plots THz waves in three-dimension which is resonantly driven at the wavelength of 1,395 nm. It is found that THz waves are linearly polarized. The polarization direction is in accordance with the optical axis. Figure 5B confirms this result at φ = 10 °   and φ = 50 ° . Combined with the polarization-independent pattern shown in Figure 4I, it provides a simple method to control the polarization direction of THz waves by rotating the optical axis of the LiNbO₃ metasurface. Figure 5C further demonstrates that the polarization of THz waves is determined by the orientation of the optical axis of the LiNbO₃ metasurface due to the prominent contribution of d₃₃. This analysis demonstrates that for high efficient enhancement of THz waves, the excitation field within each nanoresonator should be predominantly polarized along the optical axis.