Colossal infrared nonlinear optical anisotropy in a 2D charge-transfer Mott insulator

TL;DR

A 2D material called VOCl shows extremely high infrared nonlinear optical anisotropy due to its unique electronic structure and symmetry.

Contribution

The study reveals a synergistic effect of charge-transfer Mott insulator behavior and symmetry breaking in VOCl, leading to a record THG anisotropy.

Findings

VOCl exhibits a THG anisotropy ratio of 187 at 1280 nm, the highest in van der Waals materials.

THG anisotropy increases 72-fold across a broadband IR range from 2028 to 1280 nm.

Layer-independent χ(3) and band structures in VOCl are due to weak interlayer coupling.

Abstract

Mott insulators are a unique class of materials whose insulating state originates from strong electron-electron correlations: the interactions localize charge carriers, and the resulting on-site Coulomb repulsion opens a charge gap, fundamentally different from conventional insulators, making these systems an exceptional platform for exploring exotic physical phenomena. Significantly, the interplay between strong correlations and charge transfer not only stabilizes the antiferromagnetic ground state but also endows the material with enriched properties, particularly in optics. Herein, we demonstrate a 2D antiferromagnetic charge-transfer Mott insulator, Vanadium Oxychloride (VOCl), which shows giant third-harmonic generation (THG) anisotropy (ρTHG = Ix/Iy, where Ix and Iy represent the THG intensities corresponding to the excitation polarization parallel to crystal’s x- and y-axes), with ρTHG reaching up to 187 at 1280 nm excitation wavelength. Notably, it is the highest THG anisotropic ratio within the van der Waals materials family. The nonlinear anisotropy is further modulated across a broadband infrared (IR) excitation wavelength range from 2028 to 1280 nm, during which ρTHG rises from 2.6 to 187, corresponding to a 72-fold enhancement relative to its value at 2028 nm. Additionally, VOCl demonstrates layer-independent third-order susceptibilities (χ(3) ~ 10-19 m2/V2) and band structures attributed to its extremely weak interlayer electronic coupling. Moreover, the colossal THG anisotropic ratio in 2D VOCl can be ascribed to the synergistic effect of the correlated charge-transfer Mott insulator behavior and intrinsic C3 symmetry breaking, as supported by theoretical calculations. The colossal nonlinear optical anisotropy in 2D VOCl positions it as an excellent candidate for nanophotonic and optoelectronic applications, enabling next-generation nanodevices based on 2D correlated Mott insulators.

Giant THG anisotropy in 2D VOCl emerges from a synergistic coupling of charge-transfer Mott insulator behavior with intrinsic C3 symmetry breaking.

Figures

Crystal structure and synthesis of the 2D semiconductor VOCl.a Schematic of the charge-transfer Mott insulator state using the Hubbard U treatment, where the parameters Δ and U indicate the charge-transfer gap and on-site Coulomb interaction, respectively. b Structure of 2D VOCl along the *a*-axis. c Atomic-resolution STEM image of exfoliated VOCl flakes viewed along the [001], showing interplanar distances of the (200) and (020) planes as 0.169 nm and 0.192 nm, respectively. d–f Orbital-projected band structures for V 3 d, O 2p, and Cl 3p orbitals

THG process in VOCl.a Schematic representation of the THG response in VOCl flakes. b THG process-dependent on excitation intensity at 1558 nm. c THG intensity vs. the excitation intensity. The slope of 3.0 implies the occurrence of a THG process. d Variation of THG intensity with excitation wavelength. e Variation of χ^(3)^ with thickness. f Band gap vs. the number of layers in VOCl flakes. g Calculated electronic structures for monolayer (left) and bulk (right) VOCl

Anisotropic THG process in VOCl flakes.a THG mapping for VOCl flakes collected under 1558 nm with parallel (top) and perpendicular (bottom) geometry. b, c Angle-dependent polarized-THG spectra at 1558 nm. Red (blue) spots represent the *x* (*y*)-component of the THG process, respectively. 0° indicates the direction along the *x*-axis. The solid lines correspond to theoretical fitting curves, indicating the giant THG anisotropy. d Polarized *x*(*y*)-component THG spectra (*I*_x_ and *I*_y_) collected under different excitation wavelength (from left to right: 1280, 1470, 1583, 2082 nm, respectiv

Interpretation of the large in-plane anisotropy of THG in layered VOCl.a $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document}$ for 1 ~ 4 layers and bulk VOCl. b, c The THG magnitude susceptibilities of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mat

Full text

Introduction

Nonlinear optical processes and their anisotropy are critical components in advanced optical technologies, serving as the foundation for the efficiency and capabilities of modern optics.^1–7^ In particular, two-dimensional (2D) materials’ nonlinear optics and anisotropy open a new platform for researchers to modulate and manipulate light, unlocking a broad array of applications in nanoscale photonics, such as quantum photonics, on-chip photonics, ultrafast photonic devices, data security chips, signal processors, all-optical switching, and multiplexers.^8–15^ However, as the photonic industry rapidly evolves, 2D materials with strong nonlinear responses and large optical anisotropy are urgently required to meet the increasing performance demands of nanophotonic devices.

Notably, Mott insulators exhibit superior nonlinear optical properties compared to traditional materials, making them promising candidates in the field of nonlinear optics.^16–19^ These materials feature strong electron correlations that lead to electron localization and a charge gap arising from intense on-site Coulomb interactions—distinct from conventional band insulators.^20–22^ This unique electronic structure results in an interplay across charge and spin manifolds, making the optical excited states distinctly different from those in traditional insulators.^22^ As a result, Mott insulators display a range of novel phenomena in optical and other related fields.^23^ Their enhanced nonlinear optical response primarily stems from the combined effects of electron correlation and charge transfer.^24^ Neither factor alone is sufficient to enhance nonlinear optical effects, however, when both coexist, the material’s ground state can enter an antiferromagnetic state, while the excited state undergoes charge transfer, leading to a significant enhancement in nonlinear optical effects.^24^

This study reveals that the 2D charge-transfer Mott insulator (CTMI), VOCl, shows strong nonlinear optical process, and giant linear and nonlinear optical anisotropy in the infrared regime. Because VOCl possesses inversion symmetry, electric-dipole second-harmonic generation (SHG) is forbidden,^25,26^ we therefore focus on its THG response. VOCl exhibits strong third-order susceptibility (χ^(3)^ ~ 10^-19^ m^2^/V^2^) that is independent of layer number, with its band gap (~2.0 eV) displaying a similar characteristic, which suggests that VOCl has extremely weak interlayer electronic coupling. Notably, VOCl nanoflakes exhibit colossal linear and nonlinear optical anisotropy. The anisotropic ratio for photoluminescence (PL) and THG intensities reaches up to 6.8 and 187, respectively. Moreover, as the excitation wavelength reduces from 2028 to 1280 nm, the THG anisotropic ratio increases by a factor of 72. As far as we are aware, the THG anisotropy ratio represents the highest reported to date, surpassing those observed in other van der Waals materials. Based on theoretical calculations, the enormous nonlinear optical anisotropy can be ascribed to the synergistic effect of the correlated charge transfer in the Mott insulator state and intrinsic C3 symmetry breaking. Our research underscores the promise of the 2D antiferromagnetic CTMI VOCl, with its strong nonlinear optical anisotropy, for on-chip nanophotonic devices, while simultaneously advancing the development of 2D Mott insulators in nonlinear optics.

Results

Crystal characterization of VOCl

As previously mentioned, Mott insulators are a class of materials with exceptional nonlinear optical properties.^17,24^ Here, density functional theory (DFT)-based simulations were conducted to investigate the 2D antiferromagnetic semiconductor VOCl. Its crystal symmetry is orthorhombic Pmmn (No.59) which lacks C3 rotational symmetry. As shown in Fig. 1b, VOCl displays a layered structure viewed along the c-axis. In the asymmetric unit, a V atom binds with two Cl and four O atoms to form deformed octahedra VO_4_Cl_2_, named A-type and B-type octahedron (Figs. 1b and S1). Octahedral chains are built by the A-type (B-type) VO_4_Cl_2_ octahedrons by sharing Cl atoms along the b-axis, and the octahedral chains connect to each other via Cl and O atoms (shared edges in the octahedron) to form A-type (B-type) octahedral plane along the a-axis (Figs. 1b and S1). Finally, the monolayer VOCl is constructed by A-type and B-type octahedral planes via sharing three O atoms (a plane in the octahedron) along the c-axis (Figs. 1b and S1), resulting in the puckered configuration of V atoms along the (001) direction (see Fig. S1). In the DFT calculation process, the lattice of VOCl parameters were optimized to a = 3.35 Å, b = 3.77 Å, and c = 15.79 Å. They match the experimental results (see Figs. S2 and 1c), in which the X-ray diffraction (XRD) pattern displays d (001) is ~0.791 nm and the scanning transmission electron microscopy (STEM) characterization displays the interplanar distances of (200) and (020) are ~0.169 and 0.192 nm, respectively. Based on theoretical calculation, VOCl exhibits metallic behavior with the treat of semilocal functional (Fig. S3), however, experimentally, it intrinsically shows a semiconductor band gap of ~2.1 eV (Figure S4). Therefore, VOCl should be a typical correlated 2D insulator. To fully consider the strong electron-electron interaction in d-orbitals, an effective Hubbard on-site energy was adopted to correct the electronic structure, realizing a Mott transition, as shown in Fig. 1a. The Hubbard U term is applied to the V ions, opening a Mott-Hubbard band gap in the V 3 d band, with the O 2p and Cl 3p band positioned between the V 3 d orbitals. As shown in Fig. 1d-f, the lowest conduction band is primarily contributed by V_dyz, while V_dxy, O_px, and Cl_py jointly contribute to the near EF valence band. For the topmost valence band, V_d orbital plays the dominant role in the whole of Brillouin zone. However, Cl_py orbitals localize around X and U points and O_px orbitals mainly distribute from S to Γ point and R to Z point. Then, the charge transfer should mainly occur between V 3 d, O 2p, and Cl 3p orbitals, with a theoretical band gap of ~ 2.0 eV, matching the experimental observations and previous report.^27^ Thus, electronic-structure analysis identifies layered VOCl as a typical correlated CTMI (see Fig. 1a).Fig. 1. Crystal structure and synthesis of the 2D semiconductor VOCl.a Schematic of the charge-transfer Mott insulator state using the Hubbard U treatment, where the parameters Δ and U indicate the charge-transfer gap and on-site Coulomb interaction, respectively. b Structure of 2D VOCl along the a-axis. c Atomic-resolution STEM image of exfoliated VOCl flakes viewed along the [001], showing interplanar distances of the (200) and (020) planes as 0.169 nm and 0.192 nm, respectively. d–f Orbital-projected band structures for V 3 d, O 2p, and Cl 3p orbitals

Inspired by the calculations, VOCl crystals with a lateral size of up to ~15 mm are successfully grown by the chemical vapor transport (CVT) approach (Fig. S2 inset). XRD pattern (Fig. S2) shows its high quality and pure phase. The X-ray photoelectron spectroscopy (XPS, Fig. S5) was adopted to quantify elemental ratios of V, O, and Cl, yielding values of 1.00:0.92:1.05. In addition, XPS spectra exhibit that V 2p1/2 and V 2p3/2 are located at 524 and 516 eV, O 1 s is positioned at 532 eV, and Cl 2p3/2 and Cl 2p1/2 are found at 199 and 200 eV, respectively. As shown in Figs. 1c and S6, the STEM images and fast Fourier-transform (FFT) pattern reveal the orthogonal symmetry of VOCl crystals and their high quality. The simulated STEM image is well fitted to the selected atomic-resolution STEM image and calculated crystal model, and all manifest the orthogonal arrangement of V atoms in the plane of (001). Moreover, the STEM EDX mapping for the exfoliated nanoflake presents the homogeneous distribution of V, O, and Cl. The magnetic characterizations for VOCl crystals are plotted in Fig. S7, which indicate the in-plane antiferromagnetic behavior for VOCl with TN ~ 79 K.^28^

THG process in VOCl

Guided by the Mott insulator nature of VOCl, the THG properties were investigated using a home-built optical system under transmission configurations (see Fig. 2a and Method).^2,29^ Fig. 2b shows the THG spectra vs. pump power under 1558 nm irradiation, giving the THG emission peak at ~519 nm (λ/3). The log-log scaling in Fig. 2c gives a slope of ~3, confirming third-order behavior. The irradiation wavelength significantly influences the THG intensity in 2D materials. Figure 2d compares spectra of 2D VOCl acquired at identical average power for different pump wavelengths, revealing a pronounced variation with pump wavelength and a peak THG response at 500 nm, corresponding to 1500 nm excitation. It agrees with the absorption peak of VOCl and the theoretically calculated band gap, suggesting the strong resonance enhancement driven by interband transitions (see Figs. 2g, S4, and S8). Furthermore, monolayer VOCl exhibits a moderate third-order nonlinear susceptibility (χ^(3)^) of ~ 1.9 × 10^-19^ m^2^/V^2^ (Fig. S9), on par with BP^30^, graphene^31^, MoS_2_, MoSe_2_, WS_2_ and MoSe_2_^32^ ( 10^-19^ m^2^/V^2^), smaller than PdPSe (6.2 × 10^-19^ m^2^/V^2^)^2^. Notably, the χ^(3)^ of VOCl remains nearly unchanged as the thickness increases, indicating the negligible interlayer coupling in VOCl (Fig. 2e). To further illustrate the suppressed interlayer interactions in VOCl, the layer-resolved electronic band structures were calculated via first-principle methods. As displayed in Figs. 2f, g, and S8, the band gap remains nearly constant ( 2.0 eV) from the monolayer to the bulk^33^, in opposition to the behavior of other 2D materials e.g., PdSe_2_^34^, PtSe_2_^35^, PtS_2_^29^, BP^36^, where strong interlayer coupling leads to significant changes in their electronic structures and band gaps. Additionally, the layer-dependent Raman spectra show almost no shift in Raman peaks with increasing thickness (see Fig. S10), further confirming the weak interlayer electronic coupling in VOCl. This behavior is distinctly different from other 2D materials and can be attributed to VOCl’s unique structural characteristics. In VOCl’s layered structure (see Fig. 1b), the intralayer connections are primarily covalent bonds formed by V-O atoms, creating a puckered 2D framework. The Cl atoms, which encase the V-O corrugations, become inert after extracting electrons from V atoms, forming an inert outer shell that separates the layers and weakens interlayer interactions. This aligns well with the orbital-projected band structures, where the 3p orbitals of Cl atoms have negligible contribution to the valence band maximum (Fig. 1d-f). Furthermore, it is noteworthy that the near EF valence bands of VOCl, spanning from the Г- to Z-points, are primarily contributed by V 3 d, O 2p, and a minor contribution from Cl 3p orbitals. These regions exhibit minimal dispersion, namely flat band, caused by highly localized in-plane electronic states and weak interlayer coupling in 2D VOCl. This can be attributed to its unique 2D puckered V-O configuration and the inherent characteristics of a correlated CTMI.^21,37–39^Fig. 2THG process in VOCl.a Schematic representation of the THG response in VOCl flakes. b THG process-dependent on excitation intensity at 1558 nm. c THG intensity vs. the excitation intensity. The slope of 3.0 implies the occurrence of a THG process. d Variation of THG intensity with excitation wavelength. e Variation of χ^(3)^ with thickness. f Band gap vs. the number of layers in VOCl flakes. g Calculated electronic structures for monolayer (left) and bulk (right) VOCl

Colossal THG anisotropy in 2D VOCl

Motivated by electron localization and flat bands arising from strong electron correlations,^21,22^ along with C3 symmetry breaking,^40^ 2D CTMI VOCl is expected to exhibit exceptional optical anisotropy. As shown in Fig. S11, VOCl displays strong anisotropic linear optical properties. Its PL peak is located at ~2.15 eV, consistent with the experimental and theoretical band gap (Figs. S4 and 2f). The polarization-resolved PL was collected using a parallel configuration, namely, the excitation and detection polarizations are parallel, while the rotation angle θ varies from 0° to 360° (in the work, all θ = 0° means the excitation polarization direction is aligned with VOCl’s b-axis, see Fig. S12). As shown in Fig. S11b, angle-dependent PL intensities for exfoliated VOCl flakes exhibit two-fold symmetry, with the PL intensity reaching a maximum at θ = 0°. Notably, the anisotropic ratio of PL intensities is calculated using the equation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho =\,\frac{{I}_{\max }}{{I}_{\min }}$$\end{document} , and reaches up to 6.8, surpassing that of many other 2D materials, such as GaTe (5.0)^41^, ReS_2_ (4.3)^13,42^, TiS_3_ (3.2)^43^, and GeAs (1.3)^13,44^. Figure S11c displays the angle-dependent reflectance spectra spanning 450-950 nm, revealing that VOCl flakes have two absorption peaks located at ~ 565 nm and 818 nm, and two reflectance peaks located at ~680 and 850 nm. The absorption peak position of 565 nm corresponds to 2.19 eV, slightly larger than the PL peak of ~2.15 eV. Meanwhile, over 450-950 nm, the reflectance for b-axis-polarized light is lower than that for a-axis-polarized light, which suggests that more photons are absorbed under b-axis-aligned excitation polarization. This is consistent with the polarization-dependent PL spectra and absorption spectra (Fig. S4a). Additionally, the polarization-dependent reflectance intensity at a wavelength of 850 nm is extracted from Fig. S11c, which shows the typical two-fold symmetry and suggests strong optical anisotropy with a ratio of ~1.4 (Fig. S11d). Furthermore, angle-resolved polarized Raman spectra (ARPRS) for 2D VOCl flakes were also collected under parallel and perpendicular geometry, respectively (Fig. S13). As shown in Fig. S13, all peaks under the parallel configuration demonstrate two-fold symmetry. These polarized linear optical properties indicate the presence of highly localized valence-band states and the breaking of C3 symmetry.

Beyond linear optical anisotropy, the anisotropy of THG was also examined, driven by the highly localized electronic CTMI states and C3 symmetry breaking in 2D VOCl.^40^ It has a Pmmn (NO.59) space group and mmm point group, thus THG susceptibility tensor of VOCl can be expressed as^25,45^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }^{(3)}=\left[\begin{array}{cccccccccc}{\chi }_{11}^{(3)} & 0 & 0 & 0 & 0 & {\chi }_{16}^{(3)} & 0 & {\chi }_{18}^{(3)} & 0 & 0\\ 0 & {\chi }_{22}^{(3)} & 0 & {\chi }_{24}^{(3)} & 0 & 0 & 0 & 0 & {\chi }_{29}^{(3)} & 0\\ 0 & 0 & {\chi }_{33}^{(3)} & 0 & {\chi }_{35}^{(3)} & 0 & {\chi }_{37}^{(3)} & 0 & 0 & 0\end{array}\right]$$\end{document}

where the first subscript (1-3) stands for the x (namely, b-axis in VOCl crystal structure), y (a-axis), and z (c-axis), respectively. The second index specifies the term defined below:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{cccccccccc}xxx & yyy & zzz & yzz & yyz & xzz & xxz & xyy & xxy & xyz\\ 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 0\end{array}$$\end{document}

For the incident light, its electric field is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\vec{\rm{E}}}}={E}_{{x}}{\vec{\rm{x}}}+{E}_{y}{\vec{\rm{y}}}+{E}_{z}{\vec{\rm{z}}}$$\end{document} . Since the sample is pumped in the x-y plane in our geometry, the z component is neglected, thus the incident light can be further given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\vec{\rm{E}}}={{E}}_{0}({\vec{\rm{x}}}\,\cos (\theta )+{\vec{\rm{y}}}\,\sin (\theta ))$$\end{document} , where θ denotes the polarization angle (0° along x).

The χ^(3)^-induced polarization component in VOCl is expressed as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{P}}^{(3\omega )}=\left[\begin{array}{c}{{P}}_{x}^{(3\omega )}\\ {{P}}_{y}^{(3\omega )}\\ {{P}}_{z}^{(3\omega )}\end{array}\right]={\varepsilon }_{0}{E}_{0}^{3}\left[\begin{array}{c}{\chi }_{11}^{(3)}{\cos }^{3}(\theta )+3{\chi }_{18}^{(3)}\cos (\theta ){\sin }^{2}(\theta )\\ {\chi }_{22}^{(3)}{\sin }^{3}(\theta )+3{\chi }_{29}^{(3)}\sin (\theta ){\cos }^{2}(\theta )\\ 0\end{array}\right]$$\end{document}

Thus, the THG intensity resolved into x- and y-polarized components is given by:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{x}^{(3\omega )}\propto ({\chi }_{11}^{(3)}{\cos }^{3}(\theta )+3{\chi }_{18}^{(3)}\cos (\theta ){\sin }^{2}{(\theta ))}^{2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{y}^{(3\omega )}\propto ({\chi }_{22}^{(3)}{\sin }^{3}(\theta )+3{\chi }_{29}^{(3)}\sin (\theta ){\cos }^{2}{(\theta ))}^{2}$$\end{document}

The detected total THG intensity should be expressed as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I={I}_{x}+{I}_{y}$$\end{document}

As shown in Fig. S14, the total THG intensity vs. polarization angle was explored, revealing giant nonlinear optical anisotropy with two-fold symmetry. Subsequently, the polarized THG images were collected without a polarizer in the detection path. As plotted in Fig. 3a, under b-axis-aligned pump polarization, the THG signal reaches its maximum. In contrast, under b-axis-perpendicular polarization, i.e., parallel to the a-axis, the THG signal is at its lowest intensity. This behavior is consistent with the polarized PL and reflection polarization measurements. Figure 3b, c plot the x- and y-polarized THG responses on the incident polarization angle for the VOCl nanoflake under 1558 nm excitation. Equations (3) and (4) were used to model the measured results theoretically, producing good agreement with the data. Thus, the relative amplitudes of χ_11_, χ_18_, χ_22_, and χ_29_ were extracted with the ratio of χ_11_ : χ_18_ : χ_22_ : χ_29_ = 1:0.05:0.02:0.026 based on 2D VOCl with a thickness of ~35 nm. Furthermore, the THG ellipsometry consistently shows a zero value for the pump polarization angles spanning 0°-90°, indicating a linearly polarized THG response (Fig. S15a). Additionally, the orientation of the THG polarization aligns well with the incident polarization parallel to x- and y-axes, respectively (Fig. S15b).Fig. 3. Anisotropic THG process in VOCl flakes.a THG mapping for VOCl flakes collected under 1558 nm with parallel (top) and perpendicular (bottom) geometry. b, c Angle-dependent polarized-THG spectra at 1558 nm. Red (blue) spots represent the x (y)-component of the THG process, respectively. 0° indicates the direction along the x-axis. The solid lines correspond to theoretical fitting curves, indicating the giant THG anisotropy. d Polarized x(y)-component THG spectra (Ix and Iy) collected under different excitation wavelength (from left to right: 1280, 1470, 1583, 2082 nm, respectively). e THG anisotropic ratio as a function of excitation wavelength for a few-layer VOCl flake. The maximum THG anisotropic ratio is ~187 under 1280 nm, the highest reported value. f Comparison of THG anisotropic ratios for VOCl and other materials^2,15,53–65^. The shaded region indicates that the THG anisotropic ratio is below 10

The THG anisotropic ratio ρTHG is defined by the equation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rho }_{{THG}}=\,\frac{{I}_{x}}{{I}_{y}}$$\end{document} , which is proportional to the squared ratio between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{11}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{22}$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rho }_{{THG}}\propto {\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document} ). In this equation, Ix and Iy represent the THG intensities corresponding to the excitation polarized parallel to x- and y-axes (namely, VOCl’s b- and a-axes, respectively). As shown in Fig. 3b, c, Ix is much larger than Iy, suggesting a giant THG anisotropic ratio ρTHG ~ 43 for VOCl (~ 35 nm) at a 1558 nm excitation wavelength. According to the anisotropic band dispersion and optical selection rules (see Fig. 2g), the excitation wavelength is an effective strategy to modulate the THG anisotropic ratio ρTHG. Figure 3d demonstrates the THG intensities of the x- and y-components under selected excitation wavelengths (1280, 1470, 1583, 2082 nm), implying effective modulation of ρTHG by varying the excitation wavelength. Figure 3e illustrates the dependence of ρTHG on the excitation wavelength. It is noteworthy that ρTHG reaches a maximum value of ~187 under 1280 nm excitation, which is, to date, the largest reported value compared to other materials (see Fig. 3f) observed by far. This underscores the substantial leverage that the excitation wavelength exerts on ρTHG, enabling up to a 72-fold modulation, providing strong support for the application of 2D VOCl across a broad wavelength range. Moreover, the ρTHG behavior can be divided into two distinct regimes. In phase I (2082 nm down to ~1550 nm), ρTHG gradually increases, whereas in phase II (~1550 to 1280 nm), ρTHG rises sharply as the wavelength decreases. The striking difference in wavelength-dependent anisotropy between the two phases will be elucidated in the following theoretical analyses. Additionally, the thickness dependence of ρTHG was also explored under a 1558 nm excitation wavelength (see Fig. S16). As the thickness increases, ρTHG declines, which may be attributed to two factors. First, quantum confinement in ultrathin VOCl gives rise to more localized interband transitions concentrated along one in-plane direction, producing an anisotropic band structure and high THG anisotropy. As the thickness grows (i.e., through interlayer stacking), this confinement diminishes, broadening the electronic states that shape the optical response and thereby reducing anisotropy. Although VOCl exhibits interlayer electronic coupling (its bandgap and χ^(3)^ remain nearly unchanged with the thickness), the stacking nonetheless introduces slight mixing of layer-specific local fields and transitions, gradually enhancing isotropy in thicker samples. In essence, thicker flakes adopt a more three-dimensional arrangement, creating an “averaged” in-plane optical response that raises isotropy and thus lowers ρTHG.

To unveil the underlying mechanism of the high THG anisotropy observed in the 2D VOCl, we calculated its THG response using a tight-binding Hamiltonian derived from DFT wavefunctions. Based on the relationship of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\rho }_{{THG}}\propto {\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document} , Fig. 4a shows \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document} for 1 L to 4 L and bulk VOCl, where the tendency of the ratio with the excitation wavelengths shows a good agreement with our experimental observations (see Fig. 3e). As excitation wavelength shortens (but remains in phase I), both χ_11_ and χ_22_ all increase (see Fig. 4b and c), yet their ratio \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document} remains modest (only a few to tens), accounting for the weak anisotropy in phase I of Fig. 3e. Here, the third-harmonic photon energy in phase I (~ 1.69-2.30 eV) lies below the interband transition, thus unable to trigger charge transfer from O_p and Cl_p to V_d orbitals (Figs. 1a and d–f), leading to a small ρTHG. By contrast, once the excitation wavelength decreases further (phase II, ~ 2.3-3.0 eV), χ_11_ increases rapidly while χ_22_ remains unchanged or even diminishes, causing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document} to rise sharply and peak around 1200 nm (see Fig. 4a-c). To reveal the origin of the high THG anisotropy, the transition dipole moment (TDM) was calculated along a high symmetric k-point, as shown in Fig. 4d. The optical transitions demonstrate a large anisotropy between the x and y components around the X, Y, U, and T points, marked with light blue shaded regions. The large and expanded TDM in the x direction corresponds to the band transitions indicated by red arrows (see Fig. 4e), explaining the swift increase in ρTHG during phase II and its peak near 1200 nm. From the orbital-projected band structure, the transitions primarily occur from V_dxy, O_px, Cl_py orbitals to V_dyz orbitals (see Fig. 1d–f) ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|{d}_{xy},{p}_{y}\rangle +|{d}_{xy},{p}_{x}\rangle \to |{d}_{yz}\rangle$$\end{document} ), accompanied by a ~ 3.0 eV charge-transfer gap matching the THG response at ~1200 nm. Based on the composition of the orbital-resolved band structure, we define layered VOCl as a typical CTMI (Fig. 1a), which is regarded as a potential candidate for generating strong nonlinear optical responses and significant nonlinear optical anisotropy.^24^ To clarify the role of charge transfer in the THG process and its anisotropy, the charge populations of the ground state and excited state on each atom were extracted, as shown in Fig. 4f. An evident decrease in charge population is observed in the excited state relative to the ground state. Moreover, the decrease in charge population of the Cl and O atoms is larger than that of the V atom, signifying a substantial p → d transfer at resonance energies around 3.0 eV, consistent with the slight p-type conductivity in VOCl revealed by electrical transport measurements (Fig. S17). Thus, combining the analyses of orbital-resolved band structures and optical transitions (see Figs. 1d–f and 4d–f), it can be concluded that when the photon energy is tuned to the scale of the charge-transfer gap, a significant charge transfer occurs in the excited state between p-orbitals and d-orbitals, driving a rapid rise and exceptionally high THG anisotropy near 1200 nm. Excitation-induced charge transfer preserves the ground-state antiferromagnetic ordering, therefore, the giant THG ratio and THG response present a strong correlation with such separate phases.Fig. 4. Interpretation of the large in-plane anisotropy of THG in layered VOCl.a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\left(\frac{{\chi }_{11}}{{\chi }_{22}}\right)}^{2}$$\end{document} for 1 ~ 4 layers and bulk VOCl. b, c The THG magnitude susceptibilities of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{11}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{22}$$\end{document} components for 1 ~ 4 layers and bulk VOCl. d The x- and y-component transition dipole moments. The light blue shaded regions indicate the optical transitions that contribute to the THG around 1200 nm. e The optical transitions contributing to the THG anisotropic ratio peak located around the X, Y, U and T high symmetry k-point, as indicated by red arrows, and corresponding to the transition vector labeled in (d). f The ground- and excited-state charge populations on V, O, and Cl atoms

Discussion

In brief, we report a 2D CTMI VOCl that exhibits substantial linear and nonlinear optical anisotropy. First-principle calculations were applied to demonstrate the correlated Mott insulator nature for 2D VOCl. The magnetic properties were measured to emphasize its antiferromagnetic features. Comprehensive characterizations were utilized to demonstrate the linear and nonlinear optical properties, illustrating giant anisotropy ratios in PL and THG, with values of 6.7 and 187, respectively. It is worth noting that this is the largest THG anisotropy ratio reported to date. Furthermore, the THG anisotropic ratio can be significantly enhanced from 2.6 to 187 as the excitation wavelength shifts from 2028 to 1280 nm, resulting in 72 times improvement. The giant, tunable THG anisotropy of VOCl makes them promising building blocks for integrated polarization beam splitters, polarized long-wave upconversion photodetector, and polarized ultrafast lasers, opening new avenues for next-generation polarization-engineered nanophotonic devices. Furthermore, theoretical calculations indicate that the giant THG anisotropy can be explained by the synergistic effect of the CTMI behavior and C3 symmetry breaking. Besides, VOCl exhibits extremely weak interlayer coupling, as demonstrated by its thickness-independent band gap and χ^(3)^. Taking a broader view, our result provides a pathway for 2D correlated Mott insulators in nanophotonic and optoelectronic applications.

Materials and methods

Crystal preparation

Bulk crystals VOCl were synthesized via a CVT approach. A quartz ampoule was loaded with 0.5 g of V, V_2_O_5_, and VCl_3_ powder (molar ratio of 4:3:5), evacuated to < 10^-2 ^Pa, flame-sealed, and subsequently positioned in the center of a horizontal double-zone furnace. Then, the two ends of the furnace were heated to 850 °C and 950 °C over 30 hours, respectively. After two weeks, bulk crystals with a brown color were obtained in the cold zone. The obtained single crystals were mechanically exfoliated onto a SiO_2_/Si substrate (285 nm SiO_2_ film), leaving VOCl nanoflakes with different thicknesses and sizes on the substrate.

Materials characterization

Powder XRD measurement was performed on a Bruker D8 Advance powder diffractometer (λ = 1.5406 Å, Cu−Kα radiation). Elemental analysis of the bulk crystals was conducted via energy-dispersive X-ray spectroscopy (EDX, Oxford INCA) integrated into a field emission scanning electron microscope (FESEM, JEOL JSM-7600F). XPS measurements were carried out using a Kratos AXIS Supra equipped with a dual-anode Al-Kα (1486.6 eV) monochromator, and energies were referenced to C1s = 284.8 eV. AFM imaging was carried out on an Asylum Research Cypher S and a Veeco Dimension 3100. The atomic structure of exfoliated VOCl flakes was recorded using a modified JEOL 2100 F scanning transmission electron microscope (STEM). PL and Raman spectra were obtained with aid of a WITec CRM200 confocal system (λex = 532 nm). Absorption was determined by comparing the transmission through VOCl flakes with that of a fused silica substrate. Reflection spectrum was obtained with a Fourier transform spectrometer featuring a cooled mercury- cadmium-telluride (MCT) detector.

THG process measurements

THG process was evaluated by a custom-designed optical setup capable of operating in both reflective and transmissive configurations.^2,29^ The nonlinear optical responses were probed with a femtosecond laser system tunable over 650-2600 nm, delivering 250 fs pulses at 100 kHz repetition rate. The excitation laser was directed through a neutral density (ND) filter, a polarizer, a half-wave plate, and a NA 0.45, 50× objective lens, resulting in a ~ 3 μm focused spot size. The transmitted laser was obtained by a NA = 0.45, 20× objective lens, while a short-pass filter was employed to suppress the excitation light. The nonlinear optical process was subsequently collected by a spectrometer equipped with a cooled detector. All THG process measurements are collected at room temperature.

First-principles calculations

The simulations of the third harmonic generation (THG) in layered VOCl were performed using density functional theory (DFT) coupled to Wannier-based interpolation. In DFT calculations, a Hubbard correction of U = 4.5 eV was adopted to obtain the nature of the semiconductor for VOCl within the framework of PBE-GGA functional^46^ implemented in VASP package^47,48^. Structural relaxations were deemed converged when the absolute force value on any atom was < 0.01 eV/ Å. optB86b-vdW functional^49,50^ scheme was utilized to assess van der Waals interlayer interactions. Γ-centered Brillouin-zone samplings of 15 × 15 × 3 (bulk) and 15 × 15 × 1 (layered VOCl) were adopted. The cutoff energy was set to 500 eV and spin-orbit coupling was included for all DFT calculations. Using maximally localized Wannier functions, the tight-binding Hamiltonian was constructed using Wannier90 package^51^ to reproduce the band structure of DFT.

Based on the constructed tight-binding Hamiltonian, the calculations of third-order optical response were implemented in ERETCAD package. Nonlinear susceptibility ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }^{(3)}$$\end{document} ) is composed of interband contribution ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{\chi }^{(3)}$$\end{document} ) and intraband contribution ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }_{\sigma }^{(3)}$$\end{document} )^52^:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\chi }^{(3)}={\chi }_{\chi }^{(3)}+{\chi }_{\sigma }^{(3)}$$\end{document}

Alternatively, the interband contribution is described as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{rcl}\frac{{\chi }_{\chi }^{(3)}}{C} & = & \mathop{\sum }\limits_{l,m,n,p,{\bf{k}}}\frac{{r}_{mn}^{d}}{{\omega }_{nm}-{\omega }_{3}}\left[\frac{{r}_{nl}^{c}}{{\omega }_{lm}-{\omega }_{2}}\left(\frac{{r}_{lp}^{b}{r}_{pm}^{a}{f}_{mp}}{{\omega }_{pm}-{\omega }_{1}}-\frac{{r}_{lp}^{a}{r}_{pm}^{b}{f}_{pl}}{{\omega }_{lp}-{\omega }_{1}}\right)-\frac{{r}_{pm}^{c}}{{\omega }_{np}-{\omega }_{2}}\left(\frac{{r}_{nl}^{b}{r}_{lp}^{a}{f}_{pl}}{{\omega }_{lp}-{\omega }_{1}}-\frac{{r}_{nl}^{a}{r}_{lp}^{b}{f}_{\mathrm{ln}}}{{\omega }_{nl}-{\omega }_{1}}\right)\right]\\ & & +i{\mathop{\sum }\limits_{l,m,n,{\bf{k}}}\frac{{r}_{mn}^{d}}{{\omega }_{nm}-{\omega }_{3}}\left[\frac{1}{{\omega }_{nm}-{\omega }_{2}}\left(\frac{{r}_{nl}^{b}{r}_{lm}^{a}{f}_{ml}}{{\omega }_{lm}-{\omega }_{1}}-\frac{{r}_{nl}^{a}{r}_{lm}^{b}{f}_{ln}}{{\omega }_{nl}-{\omega }_{1}}\right)\right]}_{;{k}^{c}}\\ & & +i\mathop{\sum }\limits_{l,m,n,{\bf{k}}}\frac{{r}_{mn}^{d}}{{\omega }_{nm}-{\omega }_{3}}\left[\frac{{r}_{nl}^{c}}{{\omega }_{lm}-{\omega }_{2}}{\left(\frac{{r}_{lm}^{a}{f}_{ml}}{{\omega }_{lm}-{\omega }_{1}}\right)}_{;{k}^{b}}-\frac{{r}_{lm}^{c}}{{\omega }_{nl}-{\omega }_{1}}{\left(\frac{{r}_{nl}^{a}{f}_{ln}}{{\omega }_{nl}-{\omega }_{1}}\right)}_{;{k}^{b}}\right]\\ & & -\mathop{\sum }\limits_{m,n,{\bf{k}}}\frac{{r}_{mn}^{f}}{{\omega }_{nm}-{\omega }_{3}}{\left[\frac{1}{{\omega }_{nm}-{\omega }_{2}}{\left(\frac{{r}_{nm}^{a}{f}_{nm}}{{\omega }_{nm}-{\omega }_{1}}\right)}_{;{k}^{b}}\right]}_{;{k}^{c}}\end{array}$$\end{document}

The intraband contribution is described as:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{array}{l}\frac{{\chi }_{\sigma }^{(3)}}{C}=\frac{1}{i{\omega }_{3}^{2}}\mathop{\sum }\limits_{l,m,n,{\bf{k}}}{\varDelta }_{nl}^{d}\frac{{r}_{nl}^{c}}{{\omega }_{ln}-{\omega }_{2}}\left(\frac{{r}_{lm}^{b}{r}_{mn}^{a}{f}_{nm}}{{\omega }_{mn}-{\omega }_{1}}-\frac{{r}_{lm}^{a}{r}_{mn}^{b}{f}_{ml}}{{\omega }_{lm}-{\omega }_{1}}\right)-\frac{1}{{\omega }_{3}^{2}{\omega }_{2}}{\mathop{\sum }\limits_{m,n,{\bf{k}}}{\varDelta }_{nm}^{d}\left(\frac{{r}_{nm}^{b}{r}_{mn}^{a}{f}_{nm}}{{\omega }_{mn}-{\omega }_{1}}\right)}_{;{k}^{c}}\\ \begin{array}{rcl} & & +\end{array}\frac{1}{{\omega }_{3}^{2}}{\mathop{\sum }\limits_{m,n,{\bf{k}}}{\varDelta }_{nm}^{d}\frac{{r}_{nm}^{c}}{{\omega }_{mn}-{\omega }_{2}}\left(\frac{{r}_{mn}^{a}{f}_{nm}}{{\omega }_{mn}-{\omega }_{1}}\right)}_{;{k}^{b}}-\frac{i}{{\omega }_{3}}\mathop{\sum }\limits_{m,n,l,{\bf{k}}}\frac{{({r}_{nl}^{c})}_{;{k}^{d}}}{{\omega }_{mn}-{\omega }_{2}}\left(\frac{{r}_{lm}^{b}{r}_{mn}^{a}{f}_{nm}}{{\omega }_{mn}-{\omega }_{1}}-\frac{{r}_{lm}^{a}{r}_{mn}^{b}{f}_{ml}}{{\omega }_{lm}-{\omega }_{1}}\right)\\ \begin{array}{rcl} & & +\frac{1}{{\omega }_{3}}\mathop{\sum }\limits_{m,n,{\bf{k}}}\frac{{({r}_{nm}^{c})}_{;{k}^{d}}}{{\omega }_{mn}-{\omega }_{2}}{\left(\frac{{r}_{mn}^{a}{f}_{nm}}{{\omega }_{mn}-{\omega }_{1}}\right)}_{;{k}^{b}}-\end{array}\frac{1}{i{\omega }_{3}{\omega }_{2}}\mathop{\sum }\limits_{m,n,{\bf{k}}}\frac{{r}_{nm}^{b}{r}_{mn}^{a}{f}_{nm}}{{\omega }_{mn}-{\omega }_{1}}\left(\frac{\partial {\varepsilon }_{nm}^{c}}{\partial {k}^{d}}-\frac{\partial {\varepsilon }_{nm}^{d}}{\partial {k}^{c}}\right)\end{array}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C={e}^{3}K/{\hslash }^{2}m$$\end{document} with K accounting for the usual factors, r is the position operator, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{mn}$$\end{document} is the frequency offset between bands m and n, and f denotes the Fermi-Dirac distribution. Furthermore, the optical transition dipole moment is evaluated from the constructed wavefunctions,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{\rm{TDM}}}_{a\to b}=\frac{i\hslash }{({E}_{b}-{E}_{a})m}\langle {\psi }_{a}|{\bf{p}}|{\psi }_{b}\rangle =e\int {\psi }_{b}^{\ast }({\bf{r}}){\bf{r}}{\psi }_{a}({\bf{r}}){{\rm{d}}}^{3}{\bf{r}}$$\end{document}

where ψa and ψb are the wavefunctions of eigenstates Ea and Eb, respectively.

Supplementary information

Supplementary Information for Colossal infrared nonlinear optical anisotropy in a 2D charge-transfer Mott insulator