Introduction

Advances in nonlinear optics empower a plethora of applications, such as attosecond light sources based on high harmonic generation (HHG) and photodetectors for sensitive terahertz (THz) detection at elevated temperatures1,2,3,4. Inherently, the nonlinear optical properties of materials are connected with their magnetic structures5,6, crystalline symmetries7,8, and electronic band topologies. In particular, nontrivial band topologies lead to exotic electronic dynamics and enhanced optical responses9,10,11,12,13,14. Notable examples include anomalous HHG in various classes of topological materials15,16,17,18. The observation of enhanced optical responses in topological materials have been found primarily in three-dimensional systems until now12,13,15,16,17,19. Designing two-dimensional (2D) platforms with strong optical responses is advantageous for optoelectronic applications at the nanoscale with easy controllability and scalability, but so far is limited to topologically trivial materials such as graphene20 and 2H-phase transition metal dichalcogenides (TMDs)21. A promising topologically nontrivial candidate are the monolayer Janus TMDs (JTMDs) in the distorted octahedral (1 T’) phase3. Similar to 1 T’ pristine TMDs22,23,24, 1 T’ JTMDs are topologically nontrivial with an inverted bandgap in the THz regime (tens of meV). Generally, a topologically protected band structure and small electronic bandgap result in larger Berry connections, larger electronic interband transition rate, and thus stronger optical response. In addition, by replacing the top layer chalcogen atoms (e.g., sulfur) in the monolayer 1 T’ TMDs with a different type of chalcogen (e.g., selenium), the resulting Janus structure has strong inversion asymmetry and electric polarization25,26, which can further improve the nonlinear optical response.

In this work, we report experimental observations of giant nonlinearities at THz frequencies in monolayer 1 T’ JTMDs, which are synthesized via a room-temperature atomic-layer substitution (RT-ALS) method27 under ambient conditions. It is revealed that, although the electromagnetic interaction occurs only in a single monolayer flake of 1 T’ MoSSe (~10–20 µm in transverse size), the generation of mid-infrared high harmonics, THz emission, and infrared second harmonic generation are all exceptionally efficient. Further comparison with topologically trivial TMDs and theoretical analyses indicate that the keys to such giant THz-frequency nonlinearities are strong inversion symmetry breaking and topological band mixing. Our results suggest that 1 T’ JTMDs is a promising material class that could lead to an era in THz/infrared sensing using atomically-thin materials. Our results also deepen the understanding of the fundamental mechanisms underlying strong nonlinear optical responses, which could have a profound influence in, for example, room-temperature THz detection and clean energy harvesting via the bulk photovoltaic effect2,28.

Results

Multimodal nonlinearity characterization of 1 T’ MoSSe

The schematic illustration of multimodal characterization methods is shown in Fig. 1a. Our experiments investigated the THz-frequency nonlinearities of monolayer 1 T’ MoSSe with three different techniques, i.e., high harmonic generation (HHG)29,30, THz emission spectroscopy (TES), and second harmonic generation (SHG). These techniques access nonlinear coefficients with different orders (2nd to 18th order) and spectral ranges (THz to infrared). As a comparison, we also studied the responses of monolayer 2H MoSSe, 1 T’ MoS2, and 2H MoS2 under the same measurement conditions. Such combined information unequivocally indicates giant THz-frequency nonlinearities for 1 T’ MoSSe. As shown in Fig. 1b, c, 1 T’ MoSSe and MoS2 have distorted octahedral structures, with band inversion between metal d-orbitals and chalcogen p-orbitals22. In contrast, the 2H phase is characterized by a trigonal prismatic structure and is topologically trivial. In this work, Janus 1 T’ MoSSe and 2H MoSSe (Fig. 1d) are respectively converted from 1 T’ MoS2 and 2H MoS2 by the room-temperature atomic-layer substitution method2. Highly reactive hydrogen radicals produced by a remote plasma were used to strip the top-layer sulfur atoms. Meanwhile, selenium vapor was supplied in the same low-pressure system to replace the missing sulfur, resulting in the asymmetric Janus MoSSe in 1 T’ phase and 2H phase. To confirm the fidelity of material conversion, Raman scattering measurements were performed due to their sensitivity to the crystal lattice structure (Fig. 1e). For Janus 2H MoSSe, the positions of the A1g mode (~288 cm−1) and E2g mode (~355 cm−1) are consistent with literature results2; Meanwhile, the multiple A’ modes of Janus 1 T’ MoSSe located at ~226.2 cm−1, ~298.4 cm−1, ~429.8 cm−1 agree well with the theoretical calculations as well, indicating the successful material substitution.

Fig. 1: Giant THz-frequency nonlinearities in 1 T’ MoSSe.
figure 1

a Schematic illustration of THz emission spectroscopy (TES), mid-infrared high-harmonic generation (HHG), and near-infrared second harmonic generation (SHG) in 1 T’ MoSSe. b The lattice structure of 1 T’ MoSSe, 1 T’ MoS2, 2H MoSSe, 2H MoS2. c Schematic illustration of the topological band inversion in the 1 T’ phase (left) compared with the topologically trivial band structure of the 2H phase. The colormap (red and blue) indicates the wavefunction contributions from contributed by different atomic electron orbitals (e.g., chalcogen p and metal d orbitals.), and the 1 T’ phase exhibits a hybridization between the original valence and the conduction bands. d Optical images of 1 T’ MoSSe (top) and 2H MoSSe (bottom). e Experimental and theoretical Raman spectrum of 1 T’ MoSSe and experimental Raman spectrum of 2H MoSSe.

Full size image

Efficient high-harmonic generation

We first show highly efficient HHG from a single monolayer flake of 1 T’ MoSSe. The excitation source for HHG is  mid-infrared (MIR) pulses with in-plane linear polarization at 5-µm wavelength, 1-kHz repetition rate, and ~ 20 MV/cm peak field strength (setup schematic shown in Supplementary Fig. 1). The HHG image acquired in 1 T’ MoSSe (Fig. 2a) contains at least up to 18th order response, limited by our detection scheme. The even-order HHG, which is absent in bulk TMDs30, is a direct consequence of the broken spatial symmetry of the monolayer Janus systems. We varied the incident MIR polarization and observed nearly perfect cancellation of HHG intensity at specific angles to one of the crystallographic axes, indicating the HHG signal originates from a single flake instead of an average over many flakes with random orientations (shown in Supplementary Fig. 4). This is consistent with the laser spot size (1/e2 size) ~100 µm and the sparse flake-flake spacing (shown in Supplementary Fig. 4). The HHG intensity of single flake 1 T’ MoSSe are further compared with that of millimeter-scale 2H MoS2 under the same condition. Despite the irradiated flake being generally ~10 times smaller than the laser spot, the HHG of 1 T’ MoSSe is over an order of magnitude stronger than that of the millimeter-scale 2H MoS2 with 100% coverage (Fig. 2b–d)31. The strong THz nonlinearity of 1 T’ MoSSe is further confirmed by comparing it with other reference samples (2H MoSSe and 1 T’ MoS2). Figure 2e shows the HHG spectrum of 2H MoSSe, which has much weaker even-order harmonics than those of 1 T’ MoSSe. Meanwhile, the HHG of 1 T’ MoS2, which is also topological nontrivial22, shows relatively strong odd-order harmonics but no even-order harmonics, due to the inversion symmetry (Fig. 2f). Further semi-quantitative HHG efficiency comparison with other literature16,30 is summarized in Table 1 showing clear advantages of 1 T’ MoSSe over most solid-state bulk or film samples.

Fig. 2: Efficient mid-infrared high harmonic generation (HHG) in 1 T’ MoSSe.
figure 2

a HHG images of 1 T’ MoSSe observed by CCD camera. HHG extends up to ~ 5 eV and is limited mainly by the cutoff of detection optics (e.g., aluminum mirrors and grating). b HHG spectrum of 1 T’ MoSSe shown as the blue curve is over an order of magnitude stronger than the HHG from macroscopic monolayer 2H MoS2 shown as the red dashed line. The inset shows 1 T’ MoSSe HHG intensity as a function of MIR incident polarization angles. The cancellation of a few orders at some polarization angles indicates the signal is generated from a single flake. c, d Quantitative comparison of HHG intensity between 1 T’ MoSSe and 2H MoS2. Although the coverage of the 1 T’ MoSSe sample is ~10 times lower than wafer-scale 2H MoS2, the HHG signal is enhanced by over an order of magnitude, and orders as high as 18 or more are observed in 1 T’ MoSSe while only up to 15th order can be observed in macroscopic monolayer 2H MoS2. e The HHG spectrum of 2H MoSSe taken under the same conditions. The asterisks are artifacts due to high-order spectrometer diffraction. f HHG spectrum of 1 T’ MoS2. Even orders of 2H MoSSe and 1 T’ MoS2 are weak and indetectable compared with 1 T’ MoSSe. The asterisks are artifacts due to high-order spectrometer diffraction.

Full size image
Table 1 High-harmonic generation efficiency comparison of different materials
Full size table

Enhanced terahertz emission and second-harmonic generation

Dramatic enhancements in the TES and SHG measurements further validate giant nonlinearities in 1 T’ MoSSe. Figure 3a shows the TES measurements under 800-nm laser excitation (details in Supplementary Fig. 2) on four kinds of chemical vapor deposition (CVD) grown samples (1 T’ MoSSe, 1 T’ MoS2, 2H MoSSe, and 2H MoS2), among which 1 T’ MoSSe shows distinctly higher THz emission efficiency. We do not observe a detectable signal in 1 T’ MoS2 with the same excitation fluence, consistent with its centrosymmetric structure, which forbids second-order nonlinear response. The weak TES signal in 2H MoS2 has been attributed to an inefficient surface photocurrent32,33. The augmented TES in 1 T’ MoSSe aligns with the theory that 1 T’ TMDs exhibit giant nonlinearities at THz frequencies3. The polarization analysis of the THz emission (Fig. 3b) reveals the emitted radiation is mainly polarized in the lab-frame x-direction and contains a slightly weaker y-direction component (axis definition shown in Fig. 3b). Based on our experimental configuration, the x-direction emission has contributions from both in-plane and out-of-plane photoresponses, while the y-direction emission originates only from in-plane photoresponses. Thus, the observation of emission in the y-direction indicates the existence of an in-plane current contributing to the TES signal, but the stronger emission in the x-direction indicates there are likely significant contributions as well from out-of-plane currents. Further experiments are needed to disentangle the in-plane and out-of-plane contributions. For fixed excitation fluence, the peak THz field as a function of the pump polarization exhibits a sinusoidal modulation with a periodicity of approximately \(\pi\) (Fig. 3c), reflecting the rank-two tensor nature of the photoresponses, which are second order in electric fields. Detailed analysis is included in Supplementary Figs. 5–10. Finally, the TES signal shows a linear dependence on the excitation fluence (Fig. 3d) at low fluences and continues to increase at higher fluences exceeding 60 µJ/cm2, albeit with a smaller slope. Such phenomena are likely due to the combined effects of photocurrent saturation due to carrier generation32 and nonlinearly increasing photocurrents (detailed discussion in Supplementary Note 5).

Fig. 3: THz emission spectroscopy and second harmonic generation of 1 T’ MoSSe, 1 T’ MoS2, 2H MoSSe and 2H MoS2.
figure 3

a Left plot shows a schematic of THz emission setup. 1 T’ MoSSe shows a dramatically enhanced THz emission signal compared with three other types of monolayer TMDs (right plot). b Left plot shows a schematic of polarization analysis of THz emission, which shows the emission contains a major in-plane component with a possible out-of-plane contribution (right plot). c Dependence of peak THz field on excitation polarization. No polarizers were used for the emitted THz signal. d Scaling of THz emission shows a linear dependence to incident power at low fluences and saturation at higher fluences. e Left plot shows a schematic of second harmonic generation with normally incident 800-nm excitation. The white arrow marks the incident 800-nm beam and the black arrow marks the SHG light. The right plot shows SHG intensity of five different flakes in each sample and shows SHG is enhanced in 1 T’ MoSSe and 2H MoSSe compared with 1 T’ MoS2 and 2H MoS2. f The left plot shows the schematic of the angle-resolved SHG setup that measures out-of-plane dipole. The white arrow marks the incident 800-nm beam and the black arrow marks the SHG light. The beam position at the objective back aperture is scanned perpendicular to the incident beam direction with a motorized stage, which tunes the incident angle θ accordingly. The right plot shows the angle-dependent SHG intensity ratio between p and s polarization (Ip and Is) in 1 T’ MoSSe, 2H MoSSe, and 2H MoS2. In the 1 T’ MoSSe and 2H MoSSe, the Ip/Is ratio increases at non-normal incidence angles, indicative of out-of-plane dipoles. In 2H MoS2, almost no change is observed as the incident angle varies.

Full size image

Figure 3e shows the SHG measurements on the four kinds of CVD-grown samples excited with 800-nm pulses (details in Supplementary Fig. 3). The SHG of several different flakes from each sample was measured to estimate the average intensity and flake-to-flake deviation. 2H MoS2, 1 T’ MoSSe, and 2H MoSSe show high SHG efficiency, and no SHG signal is detected in 1 T’ MoS2. In 1 T’ MoSSe and 2H MoSSe, SHG is further enhanced by a factor of 4 and 3 compared to monolayer 2H MoS2 respectively, for which high SHG efficiency has been extensively reported21,34,35,36. This highlights the importance of augmented inversion symmetry breaking in Janus structures, which improves even-order nonlinearities. The SHG efficiency in Janus-type samples is further amplified in an angle-resolved SHG measurement (Fig. 3f) that is particularly sensitive to out-of-plane dipoles25. In this experiment, the incident angle of the 800-nm fundamental beam deviates from the normal incidence so that the tilted incident beam provides a vertical electric field and interacts with the out-of-plane dipoles effectively. To exclude other geometric factors, an s-polarized SHG Is induced by an in-plane dipole with the same collection efficiency is measured and used to normalize p-polarized SHG Ip that contains out-of-plane dipole contribution at non-normal incidence. For 1 T’ MoSSe and 2H MoSSe, Ip/Is symmetrically increases as a function of the incident angle, while 2H MoS2 shows much smaller angle-dependent changes. This confirms the presence of out-of-plane dipoles in Janus-type samples.

Theoretical origin of giant terahertz-frequency nonlinearity

The experimental results above indicate that the optical nonlinearity of 1 T’ MoSSe can be orders-of-magnitude (e.g., >50 times higher for 18th order HHG; >20 times higher for TES) stronger than those of 2H MoSSe. To understand this effect, we examine the microscopic mechanism underlying the strong THz-frequency nonlinear responses in 1 T’ MoSSe. The band structures of 1 T’ MoSSe is shown in Fig. 4a. The band inversion of 1 T’ MoSSe happens around the \(\varGamma\)-point. Due to spin-orbit interaction, there is a band reopening at the \(\pm \varLambda\)-points (marked in Fig. 4a). When the Fermi level is inside the bandgap, the interband transition dipole (Berry connection) \({r}_{{mn}}(k)\equiv \langle {mk}\left|r\right|{nk}\rangle\) plays an essential role in optical processes37, because it determines the strength of the dipole interaction between electrons and the electric fields. Here \(r\) is the position operator, while \({|mk}\rangle\) is the electron wavefunction at band \(m\) and wavevector \(k\). In Fig. 4b, c, we plot \(\left|{r}_{{vc}}\left(k\right)\right|\) of 2H and 1 T’ MoSSe, where \(v\) (\(c\)) denotes the highest valence (lowest conduction) band. For 2H MoSSe, the maximum value of \(\left|{r}_{{vc}}\left(k\right)\right|\) is around \(\sim 2\,\mathring{\rm A}\) near the band-edge (\(\pm {K}\) points), while for 1 T’ MoSSe, \(\left|{r}_{{vc}}\left(k\right)\right|\) can reach \(\sim 50\,\mathring{\rm A}\) near the band-edge (\(\pm \varLambda\) points). Consequently, electrons in 1 T’ MoSSe would have stronger dipole interaction and hence faster interband transitions under light illumination. This is attributed to the topological enhancement, that is, band inversions in topological 1 T’ MoSSe lead to wavefunction hybridization and hence larger wavefunction overlap between valence and conduction bands near the band edge, which accelerates the interband transitions3,38,39. The calculated first, second, and third-order nonlinear susceptibilities of 2H and 1 T’ MoSSe are shown in Fig. 4d–f. For \(\omega \, \lesssim \, 0.5\,{{{{{\rm{{eV}}}}}}}\), the responses of 1 T’ MoSSe are significantly stronger than those of 2H MoSSe. For \(\omega \, \gtrsim \, 1\,{{{{{\rm{{eV}}}}}}}\), the responses of 1 T’ and 2H MoSSe are relatively close, consistent with experimental HHG, TES, and SHG measurements at different wavelengths. In the insets of Fig. 4d–f, we plot the \(k\)-resolved contributions \({I}^{(i)}(k)\) to the optical susceptibility at \(\omega=1\,{{{{{\rm{{eV}}}}}}}\) (see “Methods” for details). Notably, the maximum value of \({I}^{(i)}(k)\) of 1 T’ MoSSe, located around the \(\pm \varLambda\) points, is larger by orders-of-magnitude than that of 2H MoSSe. This indicates that \(k\)-points near the \(\pm \varLambda\) points, which are influenced by topological enhancement, make major contributions to the total susceptibility even at \(\omega=1\,{{{{{\rm{{eV}}}}}}}\). Note that the bandgap of 1 T’ MoSSe is on the order of \(10\,{{{{{\rm{{meV}}}}}}}\), and thus interband transitions of electrons near the \(\pm \varLambda\) points are far off-resonance with \(\omega=1\,{{{{{\rm{{eV}}}}}}}\) photons. However, the contributions near the \(\pm \varLambda\) points still dominate those at other \(k\)-points where resonant interband transitions could happen. This again suggests the importance of the topological enhancement and the large interband transition dipoles near the \(\pm \varLambda\) points. The topological enhancement gradually decays at large \(\omega\). Consequently, the optical responses could be stronger in 2H MoSSe with \(\omega \, \gtrsim \, 1\,{{{{{\rm{{eV}}}}}}}\). Other in-plane elements of the SHG tensor are shown in Supplementary Fig. 14. We also note that the theoretical calculations here should only be considered as qualitative estimations, and several issues can lead to inaccuracies. For example, the density functional theory calculations suffer from intrinsic errors regarding some electronic properties, including bandgaps. Some many-body interactions are also ignored in the calculations here. Besides, the theoretical calculations deal with ideal materials, which should be distinguished from the real samples used in experimental that are influenced by doping levels, etc. Future theoretical and experimental developments could yield more accurate information for quantitative theory-experiment comparison.

Fig. 4: Theoretical calculations of the nonlinear optical response of 2H and 1 T’ MoSSe.
figure 4

a Band structure of 1 T’ MoSSe. The energy is offset to the valence band maximum. The band edge of 1 T’ MoSSe is located at the \(\varLambda\) point, which is marked in (a). b, c The magnitude of the interband transition dipole \(\left|{r}_{{vc}}\left(k\right)\right|\) for (b) 2H and (c) 1 T’ MoSSe in the first Brillouin zone. The colormap is in the unit of \(\mathring{\rm A}\). d–f First, second, and third-order optical responses of 1 T’ and 2H MoSSe as a function of incident light frequency ω. Insets of (d–f) show the \(k\)-resolved contributions (with arbitrary units) to the total response function at \(\omega=1\) eV for (left) 2H and (right) 1 T’ MoSSe. Red (blue) arrows in (b, c) and insets of (d–f) denote the \(\pm K\) \((\pm \varLambda )\) points, which are the bandedge of 2H (1 T’) MoSSe. In insets of (d–f), the Brillouin zone is zoomed in around the \(\pm \varLambda\) points for 1 T’ MoSSe to give better visibility. g The NEP (left y axis) and photo-responsivity (right y axis) of 1 T’ JTMD THz detector as a function of the light frequency \(\omega\).

Full size image

Discussion

The giant nonlinearities of 1 T’ JTMD, corroborated by both experimental and theoretical results above, support the giant THz frequency photocurrent responses of 1 T’ JTMDs predicted by theory3 and prelude that 1 T’ JTMD could serve as efficient dark-current-free THz detectors via the nonlinear bulk photovoltaic effect10. Our calculations indicate that the intrinsic photo-responsivity and noise equivalent power of the 1 T’ JTMD THz detector can outperform many current room-temperature THz sensors based on Schottky diodes or silicon field-effect transistors40,41, albeit lower than the best pyroelectric detectors and bolometers40 (Fig. 4g, see also Supplementary Note 1 and Supplementary Fig. 11). We foresee stacking multiple monolayer 1 T’ JTMDs and using field-enhancement structures42 can further enhance the responsivity43 and enable a facile usage of this detector for THz sensing purposes.

In conclusion, we demonstrate giant nonlinear responses in monolayer 1 T’ MoSSe, a prototype Janus topological semiconductor. Comparative experiments with different crystal phases (2H vs. 1 T’) and symmetry types (Janus vs. non-Janus) indicate that 1 T’ MoSSe possesses orders-of-magnitude enhancement in HHG and second-order THz emission efficiency, and a few times enhancement in infrared SHG. Supported by theoretical calculations, our results elucidate that the remarkable enhancements originate from augmented structural asymmetry in Janus-type structures and topological band-mixing in 1 T’ phases. The boosted HHG efficiency and the high fabrication versatility27 of 1 T’ JTMDs prelude a plethora of applications in light-wave electronics44,45 in the monolayer limit. Meanwhile, the giant THz-frequency nonlinearities observed in this work could enable THz detection46,47 with a large photo-responsivity at sub-A/W level and noise equivalent power down to the \({{{\mathrm{pW}}}}/\sqrt{{{\rm{Hz}}}}\) level.

Methods

Growth of 1 T′ MoS2 monolayer flakes

The precursor K2MoS4 was prepared according to the previously reported synthesis procedures48. The growth of 1 T’ MoS2 monolayer flakes was carried out in a standard CVD furnace with a 1-inch quartz tube under atmospheric pressure. A fresh-cleaved fluorophlogopite mica substrate with K2MoS4 powders were placed in the center of the furnace. After the system was purged with Ar for 10 min, the furnace was heated up to 750 °C in 40 min with 100 sccm Ar. When the temperature of the furnace reached 750 °C, 10 sccm H2 was introduced and the flow rate of Ar was decreased to 90 sccm. After 5 min, the mica substrate was rapidly pulled out of the furnace heating zone. After cooling down to room temperature, the 1 T’ MoS2 monolayer flakes were obtained.

Synthesis of 1 T’ MoSSe monolayer flakes

The synthesis of monolayer 1 T’ MoSSe is realized by a room-temperature atomic layer substitution method27 from 1 T’ MoS249. A remote commercial inductively coupled plasma (ICP) system was used to substitute the top-layer sulfur atoms of monolayer 1 T’ MoS2 with selenium. The potassium-assisted CVD-grown monolayer 1 T’ MoS2 was placed in the middle of a quartz tube. The plasma coil placed at the upstream of CVD furnace. The distance between the sample and the plasma coil is around 10 cm, with the selenium powder placed on the other side. At the beginning of the process, the whole system was pumped down to a low mTorr to remove air in the chamber. Then, hydrogen was introduced into the system with 10 sccm and the plasma generator was ignited for 20 min. The hydrogen atoms assist the removal of the sulfur atoms on the top layer of MoS2, at the same time, the vaporized selenium filled in the vacancy of the sulfur atoms, resulting in the asymmetric Janus topological structure of MoSSe. The whole process was performed at room temperature. After the reaction, Ar gas was purged into the system with 100 sccm to remove the residual reaction gas, and the pressure was recovered to atmospheric.

HHG, THz emission, and SHG setups

For HHG, the fundamental laser beam has a wavelength of 5.0 μm with a repetition rate of 1 kHz and a pulse duration of ~70 fs. It is generated by the difference frequency of the signal and idler beam from an optical parametric amplifier pumped with a Ti: Sapphire chirped-pulse amplifier (6 mJ, 1 kHz). The fundamental MIR beam with a 10 µJ pulse energy was focused on the sample with a ZnSe lens with a focal length of 15 cm. The measurements are performed in a transmission geometry at normal incidence. Generated HHG signals transmitted through the sample are collected by a CaF2 lens and directed and dispersed in a spectrometer (Princeton Instruments HRS-300) and detected by a charge-coupled device (CCD) camera (Princeton Instruments PIXIS 400B). We note that the distorted shapes of high-order harmonics are due to chromatic aberration when focusing and imaging the HHG from the entrance slit of the spectrometer to the CCD camera.

For THz emission measurement, the ultrafast laser excitation was provided by a mode-locked Ti:Sapphire laser with pulses of 40-fs (FWHM) duration and 5.12-MHz repetition rate. After focusing the laser beam on the sample, the refocused THz radiation from the sample was detected using the electrooptic (EO) effect in a non-centrosymmetric crystal (1-mm-thick ZnTe or 258-µm-thick GaP). The induced birefringence in the EO crystal was recorded at different delay times by a laser probe pulse passing through a polarizing beamsplitter (Wollaston prism) and impinging on a balanced photodetector. The power imbalance was fed into a lock-in amplifier synchronized with modulation of the excitation beam at 320 kHz by an acousto-optic modulator. By scanning the time delay between the excitation and probe pulses, the temporal profile of the transient THz electric field could be mapped. A pair of wire-grid polarizers were used to determine the THz emission polarization.

For SHG measurement, the fundamental pulses are provided by a mode-locked Ti:Sapphire oscillator at 800-nm wavelength, 5.12-MHz repetition rate, and 40-fs pulse duration. They are focused on the sample with a ×20 objective, and the generated SHG light from a single flake is filtered and detected by a photomultiplier tube. For out-of-plane measurement, a collimated p-polarized pump beam with a 1 mm spot size is guided to the objective back aperture (D = 7.6 mm). The beam was focused at the sample with a tilted angle and generated an oscillating vertical electrical field to drive the out-of-plane dipole for SHG. The SHG (green) is collected with the same objective and analyzed by a polarizer. The beam position at the objective back aperture can be scanned along the x-direction with a motorized stage, which tunes the incident angle accordingly.

Ab initio calculations

The ab initio density functional theory (DFT)50,51 calculations are performed using the Vienna ab initio simulation package (VASP)52,53. The exchange-correlation interactions are included using generalized gradient approximation (GGA) in the form of Perdew–Burke–Ernzerhof (PBE)54. Core and valence electrons are respectively treated by projector augmented wave (PAW) method55 and plane-wave basis functions. The first Brillouin zone is sampled by a \(13\times 17\times 1\) and \(17\times 17\times 1\) \(k\)-mesh for 1 T’ and 2H structures, respectively.

Nonlinear optical susceptibility calculations

After the DFT results are obtained, a tight-binding (TB) Hamiltonian in the Wannier basis is built using the Wannier90 package56. The TB Hamiltonian is then used to interpolate the relevant properties on a denser \(k\)-mesh.

The first-order susceptibility is calculated within the velocity gauge

$${\chi }_{{ij}}^{\left(1\right)}\left(\omega {{{{{\rm{;}}}}}}\omega \right)=-\frac{{e}^{2}}{{\varepsilon }_{0}\omega }\int \frac{{d}^{3}{{{{{\boldsymbol{k}}}}}}}{{\left(2\pi \right)}^{3}}{\sum }_{{mn}}\frac{{f}_{{mn}}}{{E}_{{mn}}}\frac{{v}_{{nm}}^{i}{v}_{{mn}}^{j}}{{E}_{{mn}}-{{\hslash }}\omega }$$
(1)

Here \(m,n\) label the electron states, while \({v}_{{mn}}^{i}\equiv \left\langle m\left|{v}^{i}\right|n\right\rangle\) is the velocity operator. \({E}_{{mn}}\) and \({f}_{{mn}}\) are respectively the difference in energy and occupation number between \(\left|m\right\rangle\) and \(\left|n\right\rangle\). \({\varepsilon }_{0}\) is the vacuum permittivity. The second-order susceptibility is calculated within the length gauge57

$${\chi }_{{ijk}}^{\left(2\right)}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)={\zeta }_{{ijk}}^{{II}}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)+{\eta }_{{ijk}}^{{II}}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)+{\sigma }_{{ijk}}^{{II}}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)$$
(2)

where

$${\zeta }_{{ijk}}^{{II}}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)=\frac{{e}^{3}}{{\varepsilon }_{0}}\int \frac{{d}^{3}{{{{{\boldsymbol{k}}}}}}}{{\left(2\pi \right)}^{3}}{\sum }_{{mnl}}\frac{{r}_{{nm}}^{i}{r}_{{ml}}^{j}{r}_{{ln}}^{k}}{{E}_{{{ln}}}-{E}_{{ml}}\,}\left(\frac{{f}_{{ml}}}{{E}_{{ml}}-{{\hslash }}\omega }+\frac{{f}_{{{ln}}}}{{E}_{{{ln}}}-{{\hslash }}\omega }+\frac{2{f}_{{nm}}}{{E}_{{mn}}-{{\hslash }}\omega }\right)$$
(3)
$${\eta }_{{ijk}}^{{II}}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)=\frac{{e}^{3}}{{\varepsilon }_{0}}\int \frac{{d}^{3}{{{{{\boldsymbol{k}}}}}}}{{\left(2\pi \right)}^{3}}\left\{{\sum }_{{mnl}}{E}_{{mn}}{r}_{{nm}}^{i}\left\{{r}_{{ml}}^{j}{r}_{{{ln}}}^{k}\right\}\left[\frac{{f}_{{nl}}}{{E}_{{{ln}}}^{2}\left({E}_{{{ln}}}-{{\hslash }}\omega \right)}+\frac{{f}_{{lm}}}{{E}_{{ml}}^{2}\left({E}_{{ml}}-{{\hslash }}\omega \right)}\right]\right. \\ -\left. 8i{{\hslash }}{\sum }_{{nm}}\frac{{f}_{{nm}}{r}_{{nm}}^{i}}{{E}_{{mn}}^{2}\left({E}_{{mn}}-2{{\hslash }}\omega \right)}\left\{{\Delta }_{{mn}}^{j}{r}_{{mn}}^{k}\right\}-2{\sum }_{{nml}}\,{f}_{{nm}}{r}_{{nm}}^{i}\frac{\left\{{r}_{{ml}}^{j}{r}_{{{ln}}}^{k}\right\}\left({\omega }_{{{ln}}}-{\omega }_{{ml}}\right)}{{E}_{{mn}}^{2}({E}_{{mn}}-2{{\hslash }}\omega )}\right\}$$
(4)
$${\sigma }_{{ijk}}^{{II}}\left(2\omega {{{{{\rm{;}}}}}}\omega,\omega \right)= \frac{{e}^{3}}{{2\varepsilon }_{0}}\int \frac{{d}^{3}{{{{{\boldsymbol{k}}}}}}}{{\left(2\pi \right)}^{3}}\left\{{\sum }_{{nml}}\frac{{f}_{{nm}}}{{E}_{{mn}}^{2}({E}_{{mn}}-{{\hslash }}\omega )}\left[{E}_{{nl}}{r}_{{lm}}^{i}\left\{{r}_{{mn}}^{j}{r}_{{nl}}^{k}\right\}\right. \right. \\ -\left. \left.{E}_{{lm}}{r}_{{nl}}^{i}\left\{{r}_{{lm}}^{j}{r}_{{mn}}^{k}\right\}\right]+i{{\hslash }}{\sum }_{{nm}}\frac{{f}_{{nm}}{r}_{{nm}}^{i}}{{E}_{{mn}}^{2}\left({E}_{{mn}}-2{{\hslash }}\omega \right)}\{{\Delta }_{{mn}}^{j}{r}_{{mn}}^{k}\}\right\}$$
(5)

Here the inter-band position matrix is \({r}_{{mn}}^{i}=\frac{{{\hslash }}{v}_{{mn}}^{i}}{{{iE}}_{{mn}}}\) for \(m\, \ne \, n\) and \({r}_{{mn}}^{i}=0\) for \(m=n\). \({\Delta }_{{mn}}^{i}\equiv {v}_{{mm}}^{i}-{v}_{{nn}}^{i}\) is the difference in band velocities. Meanwhile, for two numbers \(A\) and \(B\) one has \(\left\{{AB}\right\}\equiv \frac{1}{2}({AB}+{BA})\).

The ab initio theory for calculating the third-order susceptibility is not well-developed. Here we use the velocity gauge formula58

$${\chi }_{{ijkl}}^{\left(3\right)}\left(3\omega {{{{{\rm{;}}}}}}\omega,\omega,\omega \right)= \frac{i{e}^{4}}{4{\varepsilon }_{0}{\omega }^{3}}\int \frac{{d}^{3}{{{{{\boldsymbol{k}}}}}}}{{\left(2\pi \right)}^{3}}{\sum }_{{nmpq}}\frac{{v}\,_{{nm}}^{i}}{{E}_{{mn}}-3{{\hslash }}\omega }\\ \left[\frac{{v}\, _{{mp}}^{j}\left({g}_{{qn}}^{l}{v}_{{pq}}^{k}-{g}_{{pq}}^{l}{v}_{{qn}}^{k}\right)}{{E}_{{pn}}-2{{\hslash }}\omega }+\frac{{v}\, _{{qn}}^{j}\left({g}_{{mp}}^{l}{v}_{{pq}}^{k}-{g}_{{pq}}^{l}{v}_{{mp}}^{k}\right)}{{E}_{{mq}}-2{{\hslash }}\omega }\right]$$
(6)

Here \({g}_{{mn}}^{i}\equiv \frac{{f}_{{nm}}{v}_{{mn}}^{i}}{{E}_{{nm}}-\hslash \omega }\). Equation (6) experiences a spurious divergence in the real part of \({\chi }_{{ijkl}}^{\left(3\right)}\). Therefore, we first calculate the imaginary part of \({\chi }_{{ijkl}}^{\left(3\right)}\), and then obtain the real part from the Kramers–Kronig relations59.

The \({{{{{\boldsymbol{k}}}}}}\)-resolved contributions to the total susceptibility \({I}^{\left(i\right)}\left({{{{{\boldsymbol{k}}}}}}\right)\), which are shown in the inset of Fig. 4b–d, are defined as the integrand of the Brillouin zone integration in Eqs. (1, 2, 6). The Brillouin zone integrations is performed by a \({{{{{\boldsymbol{k}}}}}}\)-mesh sampling, \(\int \frac{{d}^{3}{{{{{\boldsymbol{k}}}}}}}{{\left(2\pi \right)}^{3}}=\frac{1}{V}{\sum }_{{{{{{\boldsymbol{k}}}}}}}{w}_{{{{{{\boldsymbol{k}}}}}}}\), where \(V\) is the volume of the unit cell and \({w}_{{{{{{\boldsymbol{k}}}}}}}\) is the weight factor. Since 2D materials do not have well-defined volume, we use \(V=S{l}_{{{{{{\rm{eff}}}}}}}\), where \(S\) is the area of the 2D unit cell, while \({l}_{{{{{{\rm{eff}}}}}}}\) is taken as \(6\,\mathring{\rm A}\) for all materials. The convergence in the \({{{{{\boldsymbol{k}}}}}}\)-mesh is tested.