Terahertz (THz) technologies and devices have drawn great interest for their extensive applications in communications, signal processing, drug sensing, and nondestructive testing^[1–6]. However, the lack of nonlinear THz technologies and devices has been severely limiting the further development of THz-based fingerprint sensing, nonlinear neural networks, and imaging technologies^[7–13]. Some investigations on THz frequency based on topological insulators or Weyl semimetals were performed, where the nonlinear susceptibility may be enhanced by the carriers in the topological bands^[14,15]. However, in topological photonics, most THz-frequency research does not involve nonlinearity since, for a long time now, THz nonlinearity was considered very weak. After efficient high-harmonic generation (odd-order only) of THz waves was achieved on graphene in 2018^[16], researchers realized four-wave mixing (third-order) in topological graphene plasmon systems two years later^[17]. In 2021, we found that the THz wave could achieve giant nonlinearity on a LiNbO3 crystal based on phonon polaritons (second-order nonlinearity)^[18], which makes it possible to further achieve giant THz nonlinearity on crystal solids. In addition, topological photonics, featuring phenomenal surface states immune to disorder and impurities^[19–25], has led to new applications in the microwave, visible light, and THz frequency bands. Introduction of topological photonic systems in the THz regime has achieved THz signal transmission, polarization splitting, polariton edge transmission, and tuning of THz confinement^[26–29]. However, due to the low efficiency of the THz nonlinearity, the nonlinear harmonic generation of THz waves, especially in topological systems, is unexplored.

In this work, we have introduced a scheme for the THz harmonic generation in a 2D topological valley polaritonic microcavity. The obtained results suggest that the topological polaritonic microcavity is a good platform for generating and controlling THz harmonics. By employing the finite-difference time-domain (FDTD) method, based on the quantum valley Hall (QVH) system, we have designed a triangular microcavity using two topological valley photonic crystals (VPCs) based on polar crystals. The interface between two VPCs results in a triangular microcavity, as shown schematically in Fig. 1(a). Both centrosymmetric and noncentrosymmetric nonlinear materials, with the nonlinear susceptibility evaluated by nonlinear Huang equations, are selected to study the harmonic generation process. The harmonic waves are further tuned by changing the amplitude of the THz source power.

In our design, the honeycomb lattice is adopted to deform two different VPCs. For simplicity, a 2D model within a dielectric material with linear susceptibility, χ(1)=25, is considered, where absorption, dispersion, and nonlinearity are absent. Here, the prototype model is adopted from a LiNbO3 photonic slab, on which air-hole arrays can be fabricated^[28]. A triangular microcavity is constructed by the interface between two types of deformed honeycomb-like lattices, as shown in Fig. 1(a), where the lattice constant is a=150μm and the deformed lattice site diameters are D=80μm and d=22.6μm. The undeformed lattice is C6 symmetric, where D=d=51.4μm, which shows a Dirac cone at 0.31 THz in the band structure, and K point in the momentum space, as shown in Fig. 1(b). Here, the energy band is calculated along the boundary of the first Brillouin zone using the plane wave expansion method, where the hot-color mapping represents the calculated band E˜(ω,k→). Upon deforming the honeycomb lattices so that D≠d, the lattice changes to C3 symmetric with its inversion symmetry broken. Thus, as the degeneracy of the Dirac cones at K/K′ points is broken, a bandgap originates in the distorted lattice, which is indicated in Fig. 1(b). Here, we optimized the value of D−d=0.38a to reach a broad bandgap for further nonlinear applications^[30,31].

In our settings, the bandgap of both deformed lattices ranges from 0.30 to 0.36 THz, which is well located in the 6G communication band. In general, the inequivalent lattice deformation results in different topological properties. In Fig. 1(c), we have shown the magnetic field (Hz) distribution in the first band at the K valley for VPC1 and VPC2, where two different modes are depicted. The phase vortex profile of the Hz field at the K valley in the first band is also calculated and given in Fig. 1(d), which shows a clockwise/anticlockwise rotation by a 2π phase around the center of unit cells. This phase vortex indicates a valley pseudospin in the deformed lattice, which is also known as a photonic valley degree of freedom, and is related to a topologically nontrivial phase^[28,32]. The different distribution of the field amplitude and phase vortex of Hz suggests a valley-chirality locking property, thus VPC1 and VPC2 are assigned with different Chern numbers at the K/K′ valley^[28,32].

In the FDTD calculation, a 2D topological valley polaritonic microcavity with χ(1)=25 is constructed, and a standard perfect-match layer boundary is employed to avoid the influence of the reflected THz waves. In order to clarify the optical properties of the valley edge modes and the microcavity confinement, a THz source (dipole source, y-polarized with a center frequency of 0.32 THz) is placed at the left-bottom corner of the triangular microcavity, as depicted in Fig. 1(a). This setting is optimized for the experimental reproduction since the THz waves are generated in the slab as a transverse electric mode^[28].

After calculation, a time-resolved field Ey(x,y,t) was obtained. Thus, the spectral information can be obtained through a Fourier transform of E˜y(x,y,ω)=∫Ey(x,y,t)e−iωtdt and E˜y(ω)=∫E˜y(x,y,ω)dxdy. By calculating the spectrum in the microcavity, an excellent single mode can be seen in the result, which is presented in Fig. 2(a) in logscale, and the linear scale one is also shown in the inset. Further, at 0.32 THz, the Ey field distribution is calculated to obtain the triangular topological valley microcavity mode, as shown in Fig. 2(b). This result is in good agreement with that in the previous research.

In order to reveal the nonlinear property of the topological edge microcavity and provide a scheme for efficient nonlinear harmonic generation, now we have changed the base materials from linear to nonlinear. Phonon polaritons can assist in giant enhancement of THz nonlinearity. Therefore, a noncentrosymmetric phonon-polariton system is studied with χ(1)=25, χ(2)=10−6m/V, where the second-order nonlinear susceptibility is reproduced from our previous results^[18].

Employing the similar parameters and the same method as in the linear case, the results in the noncentrosymmetric nonlinear regime can be obtained. As the amplitude of the THz dipole source increases, the nonlinear effects gradually emerge. Figure 3(a) shows the spectra in the triangular microcavity, in which one can notice the gradual growth of nonlinear harmonics as the source power increases. For the source amplitudes smaller than 1×104V/m, no nonlinear signals are observed, as shown by the blue curve in Fig. 3(a). For the source amplitudes 5×104V/m and higher, second-harmonic waves begin to appear, as shown by the red curve. Additionally, third- harmonic waves are generated as the source amplitude increases to 1×106V/m, as indicated by the green curve. Both the order numbers and the amplitude of the harmonic waves increase further, as the source power gets far stronger. For instance, for the source amplitude 5×106V/m, all the harmonics from the linear order to the 8th-order can be observed in the spectrum. By choosing the source amplitude of 5×106V/m, the mode distributions of the Ey field of all the harmonic modes can be calculated, and the obtained Ey for the linear order to the 8th-order harmonic waves are presented in Figs. 3(b)–3(i). The geometry of the topological microcavity is designed for linear localization, as shown in Fig. 2. When nonlinearity involves, although the linear mode at 0.32 THz is well localized, the high-order harmonics is not located in the topological bandgap [see Fig. 1(b)]. Thus, the high-order harmonics would travel through the whole microcavity after it was generated at the interface.

In nonlinear optics and photonics, many nonlinear materials only show odd-order nonlinearity because of their crystal lattice and potential energy being centrosymmetric. In these materials, the most significant contribution to the nonlinearity is the third-order nonlinear susceptibility, χ(3). In such materials, if the phonon-polariton system is considered, the pertinent nonlinear susceptibility associated with it can also be evaluated using nonlinear Huang equations as follows^[18,33]:$$\begin{gathered} x¨+γx˙+ω02x=(Nm)−1/2bE−βx3+(q/m)Eex, \\ P=(Nm)1/2bx+ϵ0(ε∞−1)E, \end{gathered}$$(1)where x is the amplitude of motion, ω0 represents the eigenfrequency of the ions, and γ and β stand for the decay rate and the coefficient of the nonlinear restoring force, respectively. An external driving field Eex is exerted upon the ions whose effective mass and electrical charge are given by m and q, respectively. N shows the oscillator density, whereas E and P indicate the macroscopic electric field and polarization, respectively, in the crystal. The phonon polariton’s contribution is characterized by the constant b=ω0ϵ0(ε0−ε∞), with ε0 and ε∞ representing the low- and high-frequency relative permittivity of the material, respectively. The permittivity of the vacuum is ϵ0=8.85×10−12F/m.

Compared to the nonlinear Huang equations in the noncentrosymmetric materials^[18], here the nonlinear restoring force is described by the third-order enharmonic restoring force βx3 rather than the second order one ax2. The parameters in Eq. (1) are similar to the previous results, and coefficient β=(ω0/d)2, where d is the average lattice constant^[18,34]. Here, we take parameters from the LiNbO3 crystal, and all the calculation details are similar to the second-order case that can be found in Ref. [18]. Ignoring the evaluation process, the obtained third-order nonlinear susceptibility is about 10−13m2/V2. Although the LiNbO3 crystal is noncentrosymmetric, the χ(3)=10−13m2/V2 shows a good evaluation of a phonon-polariton polar crystal.

In the FDTD calculation, the nonlinear material is set to be χ(1)=25 and χ(3)=10−13m2/V2, where all other parameters remain identical to the above calculations. In Fig. 4, the results of the nonlinear harmonics in the triangular topological valley microcavity are displayed. For the source amplitudes smaller than 5×105V/m, no nonlinear harmonics are observed. The third-harmonic generation (THG) are produced upon raising the source amplitude to 1×106V/m, as shown in Fig. 4(a).

From the Ey field distribution, obtained at the source amplitude of 3×106V/m, a good edge mode is observed around the triangular microcavity for the linear modes, as shown in Fig. 4(b). However, the mode distribution of the THG waves is shown as a bulk mode rather than as an edge mode around the microcavity, which is displayed in Fig. 4(c). Generally, the field distribution of the topological edge states induced by third harmonics should match with the topological localization pattern^[35]. It is worth noting that while the frequency of the topological edge states varied from the middle to the edge in the bandgap, the field distribution of the third harmonics that it induced evolved from the edge states to the bulk states. In this case, field distribution at 0.32 THz presented typical topological edge modes, while the third harmonics that it induced distributed as bulk modes, which is quite beyond our expectation. Here, we speculate that this is caused by some complicated mechanisms between topological photonics and nonlinear effects. This bulk cavity mode is found to be tunable, which can also be observed in the near-infrared wavelengths^[35].

In order to verify the topological robust property, we remove a unit cell from the triangular microcavity in both centrosymmetric and noncentrosymmetric nonlinear materials, as shown in Fig. 5(a). Figure 5(b) shows the topological robust property of the linear mode against this defect, while the THG signal does not maintain a good topological bulk state due to the leaking energy, as seen in Fig. 5(c). For noncentrosymmetric materials, the topological robust was well protected for linear, second-order, and third-order modes against the defect, as indicated in Figs. 5(d)–5(f).

In conslusion, we introduced a scheme for the THz harmonic generation in a triangular topological valley polaritonic microcavity, which suggests it to be a good platform for generating and controlling the THz harmonics. Using the FDTD method, a triangular microcavity is designed using two topological VPCs. Based on the phonon-polariton enhanced nonlinearity using nonlinear Huang equations, we have studied the nonlinear harmonic generation in both noncentrosymmetric and centrosymmetric nonlinear materials. The results have indicated the prominent topological microcavity edge mode of the high-harmonic generation in the noncentrosymmetric nonlinear materials and the bulk topological microcavity mode of third-harmonic waves in the centrosymmetric nonlinear materials. The topological robust properties against defects are also verified. Outcomes of our study can assist as a valuable reference for nonlinear devices, harmonic control, and 6G signal nonlinear processing in the THz regime.