Large sliding regulation in van der waals layered nonlinear optical ternary chalcogenides

Abstract

Interlayer sliding modulation in van der Waals (vdW) layered materials opens opportunities to seek tunable applications for photonics and optoelectronics at the nanoscale. In this work, focusing on ternary chalcogenides, dozens of typical vdW layered materials with the formula of A₂MZ₄ (A, M = metal cation; Z = S/Se) were systematically classified and studied to illustrate the interlayer sliding modulation performance on nonlinear optical (NLO) effect and underlying mechanism. Based on first-principles analysis, the screened Cd₂GeSe₄ exhibits large second harmonic generation (SHG) deviation (>36 times between high- and low-SHG states) under sliding among all studied materials. This phenomenon arises from the polarization modulation induced by the vdW interlayer charge redistribution of two tetrahedral primitives in the ternary edge-sharing tetrahedral structures. Our findings present a strategy to substantially modulate NLO properties by altering the interlayer charge coupling which is commonly thought to be very weak in vdW layered materials.

Introduction

Van der Waals (vdW) layered materials are receiving significant attention in nanoscale optoelectronic applications thanks to their two-dimensional (2D) functional arrangement, easy integration, and high adjustability^1,2,3,4. Recently, interlayer sliding has drawn increasing research interest as a practical modulation method for vdW layered materials. A typical example is sliding ferroelectricity^5,6,7, which has been experimentally confirmed in a series of vdW materials^{8,9,10,11,12,13,14} after being proposed by first-principles calculations in 2017¹⁵. By changing the layer stacking mode via sliding, the interlayer charge transfer in a vdW material would be manipulated and the polarization direction is altered, thus the material can be switched between ferroelectric and antiferroelectric states, promoting the emerging concept of ‘slidetronics’^13,14,15,16. More importantly, the interlayer sliding in vdW materials usually entails a small energy barrier of several to tens meV per atom in 2D phases, which can bring great convenience to device applications⁵.

Second order nonlinear optical (NLO) and ferroelectric response in a material are closely related since both effects originate from non-centrosymmetric structures and involve the response of electrons to external field^17,18,19, meaning that ferroelectric tuning approaches like sliding can also be regarded as potential ways to regulate the NLO properties in vdW layered systems. Recently, studies on sliding-modulation of NLO properties relating to electrical polarization (e.g. bulk photovoltaic effects) are beginning to emerge²⁰. However, the sliding-modulation of second harmonic generation (SHG) is rarely reported in layered systems. Some recently discovered vdW layered SHG systems do have different stacked phases introduced by interlayer sliding, but these phases show nearly identical SHG susceptibilities. For example, layered 18MR-B₂O₃ reported by Li et al. has two NLO-active stacking types distinguished by a (1/3, 2/3) layer translation with only a minor SHG difference between them led by the difference in layer distances²¹. Therefore, it arouses our curiosity to discover vdW layered materials with sliding-sensitive SHG response and to identify their modulation mechanism.

The NLO regulating properties in vdW layered materials require large changes in NLO properties triggered by external stimulators, such as stacking or twisting of layers, and external electric field or strain. It has attractive optoelectronic and photonic application prospects including signal processing, sensors, and broad-spectral optical frequency converters, etc^{22,23,24,25,26,27,28,29}. Traditional vdW layered NLO semiconductors are concentrated in binary (e.g. MoS₂ and h-BN with second order NLO properties) or single-element systems (e.g. graphene and black phosphorus with third order NLO properties), and their atomic orbitals are coupled to form a single functional group. Differing from their counterparts with single or binary elements, ternary materials have become an emerging class of vdW layered semiconductors due to more abundant adjustability³⁰. With the addition of a third element, two or multiple different functional groups coexist in a ternary material. Their synergistic effect brings an additional degree of freedom to regulate the material performances by changing the elemental compositions, which may further produce more appealing applications. Several layered ternary materials have been reported with tunable NLO properties, such as ASb₅S₈ (A = K, Rb, and Cs) system synthesized by Chen et al. in 2021²⁵ and Cu₂MoS₄ theoretically predicted by our group in 2022³¹. To achieve the NLO regulation performances, NLO regulation ratio and reversibility are two important concerns²⁵. While the NLO switch in ASb₅S₈ (A = K, Rb, and Cs) is thermal driven, Cu₂MoS₄, as a typical transition metal ternary chalcogenide (TMTC)^30,32, exhibits a vdW layered structure with two different stacking patterns (namely AA and AB type, corresponding to P- and I-phases, respectively) distinguished by a (1/2, 1/2, 0) interlayer shift^32,33,34. Both phases are SHG-active with non-centrosymmetric structures, and the second harmonic generation (SHG) strength χ⁽²⁾ can be modulated by the AA-AB stacking change via the (1/2, 1/2, 0) shift, with ${\chi }^{\left(2\right)}(\text{AB})/{\chi }^{(2)}(\text{AA})$ being ~1.7 for bilayer systems and ~2 for bulk phases, as well as a small energy barrier for the interlayer sliding (~10 meV per atom)³¹. This means that Cu₂MoS₄ has an eligible reversibility of sliding driven NLO modulation. Nevertheless, its SHG deviation induced by interlayer sliding is not sufficiently large for tunable NLO applications, limiting its NLO regulation ratio. However, Cu₂MoS₄ provides a platform for SHG modulation with the synergistic effect of two types of NLO-active tetrahedral motifs. It then becomes intriguing whether larger SHG deviation can be achieved with the same structure as Cu₂MoS₄ by changing elemental composition while keeping the benefits of easy adjustment.

In this research article, we proposed a series of P- and I-phase vdW layered A₂MZ₄ materials derived from Cu₂MoS₄ by element substitution, where A (M) represents the metal cation taking up the Cu (Mo) site in the original Cu₂MoS₄ and Z is the chalcogen element (S or Se). The SHG effects in a total of 34 proposed structures are investigated by the first-principles calculations, which demonstrate that the SHG deviation between P- and I-phases differs significantly with different elemental compositions. Large SHG deviation is found in several substituted layered A₂MZ₄ materials, topped by Cd₂GeSe₄ with a ${\chi }^{(2)}$ deviation of more than 30 times between P/I phases, showing promising potential as a tunable/switchable layered NLO material. Underlying analysis shows that the interlayer charge redistribution, originally thought to be very weak in vdW layered materials, turns out to be a vital source of the large SHG difference, revealing a mechanism via interlayer charge coupling in the modulation of NLO properties in ternary vdW layered materials.

Results and discussion

Construction of A₂MZ₄ materials

Based on the parent Cu₂MoS₄, all designed A₂MZ₄ derivatives share a porous vdW layered structure, as illustrated in Fig. 1. They all have P- and I-phases with $\text{P}\bar{4}2\text{m}$ and $\text{I}\bar{4}2\text{m}$ space groups, corresponding to two different stacking patterns (namely AA type for the P-phase and AB type for the I-phase), respectively^32,33,34. Both phases are composed of almost identical A₂MZ₄ layers formed with anti-parallel [AZ₄] and [MZ₄] tetrahedrons separated by sharing the edge Z atoms as shown in Fig. 1. In the P-phase, all atoms from the same sites in adjacent layers are piled directly above each other, while in the I-phase a (1/2, 1/2, 0) layer shift occurs, and the AB stack of layers is displayed. In total, three types of isomorphic A₂MZ₄ materials (a total of 34) are designed from the Cu₂MoS₄ structures as follows:

Fig. 1: Evolution of A₂MZ₄-type stacked structures.

Each A₂MZ₄ layer is composed of anti-parallel [AZ₄] and [MZ₄] tetrahedrons. P- and I-phases are distinguished by a (1/2, 1/2, 0) interlayer shift.

Type-I: The A, M, and Z sites are substituted by the elements from the same groups as Cu, Mo, and S in Cu₂MoS₄, respectively, i.e., A = Cu, Ag; M = Mo, W; Z = S, Se.

Type-II: The A site is occupied by bi-valence d¹⁰ cations; the M site is occupied by tetra-valence cations; and Z site is occupied by chalcogen elements, i.e., A = Cd, Hg; M = Si, Ge; Z = S, Se.

Type-III: The A site is occupied by tri-valence IIIA group cations, M site is occupied by bi-valence d¹⁰ cations; and Z site is occupied by chalcogen elements, i.e., A = Al, Ga, In; M = Cd, Hg; Z = S, Se.

In all the materials, Cu₂MoS₄, Cu₂WS₄, and Cu₂WSe₄ have been experimentally synthesized in both P- and I-phases, while Cu₂MoSe₄ is considered less stable and may not have the I-phase³⁴. Ag₂WS₄ has been synthetized in I-phase³⁵. The SHG effects in all these materials have not been reported. Other structures are obtained based on the first-principles structural optimizations (see the computational details in the Supporting Information).

Evolution of NLO properties

The calculated band gaps E_g (with HSE06 functionals) and static second order NLO susceptibility ${\chi }_{123}^{\left(2\right)}$ for all types of A₂MZ₄ derivatives are listed in Table 1. All layered A₂MZ₄ compounds have only one non-vanishing component ${\chi }_{123}^{\left(2\right)}$. It can be found that type-I derivatives have positive ${\chi }_{123}^{\left(2\right)}$, while most type-II and III structures have negative ${\chi }_{123}^{\left(2\right)}$. Generally, type-I derivatives show the similar SHG divergence in P- and I-phases as the case of original Cu₂MoS₄ that the ${\chi }_{123}^{\left(2\right)}$ for I-phase is larger than that for P-phase in type-I structures. In type–II and III structures, however, P-phase often has larger absolute values for ${\chi }_{123}^{\left(2\right)}$. Within each type of compounds, heavier A-site element often leads to decreased E_g and increased absolute values for ${\chi }_{123}^{\left(2\right)}$ in both phases, and the opposite is true for M-site except several type-III compounds.

Table 1 Calculated results for all types of layered A₂MZ₄ derivatives.

Since our core concern in the layered A₂MZ₄ system is their discrepant SHG performance distinguished by stacking order, a statistical analysis is performed to assess the ${\chi }_{123}^{\left(2\right)}$ deviation between P/I phases for all designed derivatives, by defining a Deviation Factor D as:

$$D=\left|{\chi }_{123}^{\left(2\right)}\left({\rm{P}}\right)-{\chi }_{123}^{\left(2\right)}\left({\rm{I}}\right)\right|/\min \left(\left|{\chi }_{123}^{\left(2\right)}\left({\rm{P}}\right)\right|,\left|{\chi }_{123}^{\left(2\right)}\left({\rm{I}}\right)\right|\right)$$

(1)

Note that, although this index becomes infinite for the cases where phase transitions induce inversion symmetry break, all layered A₂MZ₄ materials studied in the present study are intrinsically non-centrosymmetric, and thus the definition of D is valid for the current purposes. Basically, a larger D value stands for a greater NLO regulation ratio in the material system. For a certain A₂MZ₄ derivative, ${\chi }_{123}^{\left(2\right)}\left({\rm{P}}\right)$ and ${\chi }_{123}^{\left(2\right)}\left({\rm{I}}\right)$ are the SHG susceptibilities of its P-phase and I-phase, respectively. When ${\chi }_{123}^{\left(2\right)}\left({\rm{P}}\right)={\chi }_{123}^{\left(2\right)}\left({\rm{I}}\right)$, D takes the minimum value as 0. Larger D value indicates greater SHG deviation between P/I phases. In Fig. 2a, b, we use heat maps to compare this index for all the three types of layered A₂MZ₄ derivatives. Clearly, the largest SHG deviation between the two phases appears in type-II structures despite their smaller absolute values of ${\chi }_{123}^{\left(2\right)}$ compared to the type-I structures. Type-III structures show the smallest SHG deviation in general, among which no compounds have larger ${\chi }_{123}^{\left(2\right)}$ deviation during the P/I phase change than the original Cu₂MoS₄. It is also worth noting that replacing S with Se can enlarge the SHG divergence between P- and I-phase.

Fig. 2: Deviation Factor D of A₂MZ₄ layered derivatives.

a Deviation Factor D of A₂MS₄. b Deviation Factor D of A₂MSe₄. A deeper red color refers to a larger SHG deviation.

In all the layered A₂MZ₄ derivatives, Cd₂MZ₄ (M from IVA main group), as a typical type-II derivative family, shows the largest SHG deviation between AA and AB stacking type. For example, Cd₂SiS₄ has negative ${\chi }_{123}^{\left(2\right)}$ value in P-phase and positive ${\chi }_{123}^{\left(2\right)}$ value in I-phase, and Cd₂GeSe₄ has a remarkable SHG deviation with D value being 36.74, showing a more promising potential as tunable NLO materials than the original Cu₂MoS₄. Since both P- and I-phases in layered A₂MZ₄ structures are non-centrosymmetric and intrinsically SHG-active, it is greatly desirable to find such change and its origin from near-zero to high ${\chi }^{(2)}$ values in Cd₂GeSe₄. In addition, we calculated the SHG susceptibilities dependent of pump laser energy in the A₂MZ₄ layered materials with the largest incident photon energies less than half of the bandgap energy. The higher photon energy would lead to the strong optical absorption of SHG radiation in the materials, which is out of the scope of our present study. As shown in Supplementary Figs. 10 and 11, the frequency-dependent SHG results also support that Cd₂GeSe₄ maintains larger P/I SHG deviation compared with typical type -I and III candidates like layered Cu₂MoS₄ and layered Al₂CdS₄ throughout the whole transparent region, consistent with the trend in static calculation results.

Microscopic SHG origin in Cd₂GeSe₄

As a representative type-II derivative, Cd₂GeSe₄ has the largest SHG deviation between the P- and I-phases. To investigate the underlying mechanism, the electronic structures and band-resolved SHG contribution in corresponding orbitals for Cd₂GeSe₄ are plotted in Fig. 3. In both phases, SHG-active states span a large area in the conduction bands but are concentrated within a near-gap section (0 ~ −5 eV) in the valence bands. The deeper valence states such as Cd 4d orbitals at −8 ~ −9 eV only play a minor part in the total SHG effect, despite their large DOS distribution. Based on the orbital composition and distribution of ${\chi }_{123}^{\left(2\right)}$, the near-gap section can be divided into three regions (regions I, II and III). Region I includes the highest occupied states which mainly consist of the localized Se 4p orbitals and a small amount of Cd 4d and 5p orbitals, with the energy spanning of −0.90 ~ 0 eV and −0.41 ~ 0 eV in the P- and I-phases, respectively. The states in region II or −2.26 ~ −0.90 eV (−1.72 ~ −0.41 eV) for the P (I) phase are mainly composed of Se 4p and Cd 5p orbitals, most of which attributed to [CdSe₄] motifs. The orbitals in region III, namely −4.64 ~ −2.26 eV (−4.17 ~ −1.72 eV) for the P- (I-) phase, are mainly composed of Ge-Se bonding electrons from Ge 4p and Se 4p orbitals, most of which belong to [GeSe₄] motifs. Band-resolved SHG analysis shows that the remarkably larger negative SHG contribution are produced by the localized Se 4p orbitals in the region I of P-Cd₂GeSe₄, with the ${\chi }_{123}^{\left(2\right)}$ contribution being −9.58 pm·V⁻¹ (−0.64 pm·V⁻¹) for the P- (I-) phase. Meanwhile, the regions II and III produce negative (−2.88 pm·V⁻¹ for P-phase and −3.27 pm·V⁻¹ for I-phase) and positive (2.58 pm·V⁻¹ for P-phase and 2.21 pm·V⁻¹ for I-phase) SHG contribution, respectively. Due to the anti-parallel alignment of [CdSe₄] and [GeSe₄] motifs in Cd₂GeSe₄, we see the two types of tetrahedra counteract with each other and generate SHG polarization in opposite directions.

Fig. 3: Electronic structure and orbital distribution of NLO susceptibility of P/I-Cd₂GeSe₄.

a–c Band structures, DOS and band-resolved SHG contributions of P-Cd₂GeSe₄. d–f Band structures, DOS, and band-resolved SHG contributions of I-Cd₂GeSe₄. g–i ${\chi }_{123}^{\left(2\right)}$-weighted electron densities in regions I ~ III for P-Cd₂GeSe₄. j–l ${\chi }_{123}^{\left(2\right)}$-weighted electron densities in regions I ~ III for I-Cd₂GeSe₄. Red and blue areas in (g–l) represent the orbitals having positive and negative ${\chi }_{123}^{\left(2\right)}$ values, respectively.

For a more intuitive presentation of the three valence regions and their respective SHG contributions, the ${\chi }_{123}^{\left(2\right)}$-weighted electron density schemes of the region I, II, and III orbitals from P- and I-Cd₂GeSe₄ are displayed in Fig. 3g–l. In both structures, the major ${\chi }_{123}^{\left(2\right)}$ contributions with positive values originate from the electron clouds resembling Ge-Se bonding orbitals, while the negative ${\chi }_{123}^{\left(2\right)}$ contributions from the electron clouds around Se and Cd. The highest occupied orbital in the region I generally produce negative ${\chi }_{123}^{\left(2\right)}$ components for both phases. However, in the I-phase a small portion of the region I orbitals weighted with positive ${\chi }_{123}^{\left(2\right)}$ also appears, neutralizing part of the SHG contribution. Accordingly, one may conclude that the SHG divergence between the two phases mainly comes from region I where the Se 4p orbitals dominate SHG effect.

Sliding mechanism for SHG difference

To further unveil how the structural change between two different stacking types brings such giant SHG deviation, we first investigated the cases of monolayer and bilayer AA- /AB-stacked Cd₂GeSe₄ systems with different distances between layers. We find that the SHG difference between AA and AB stacking type decreases as the layer distance becomes larger, and eventually, the bilayer SHG degenerates into the same magnitude as the separate monolayer, as displayed in Fig. 4a. Therefore, the SHG deviation can be attributed to the coupling between layers since the layer distance directly reflects the strength of interlayer interactions.

Fig. 4: Mechanism on the sliding induced NLO regulation in vdW Cd₂GeSe₄.

a Variation of ${\chi }_{123}^{\left(2\right)}$ in the Cd₂GeSe₄ bilayer systems with different layer distances. All the bilayer systems are formed by stacking two identical monolayers to eliminate the influence of structural difference within each layer on the SHG susceptibility. b Charge density difference between adjacent layers in P-/I-Cd₂GeSe₄. The areas where the electron density has been enriched/depleted are highlighted in red/blue color. c Charge depletion (blue)/accumulation (red) and the highest occupied orbitals (transparent yellow area) shown together on the Se atoms facing adjacent layers in P-/I-Cd₂GeSe₄. d Illustration of the three-region model for band composition and SHG contribution in the layered A₂MZ₄ structures.

As an indicator of interlayer charge coupling, the charge density difference (CDD) is implemented to visualize the charge redistribution with respect to the positions of adjacent layers in Cd₂GeSe₄, which is displayed in Fig. 4b. In both phases these redistributed charges are transferred from charge depletion areas which overlap with the Se 4p electrons that constitute the region I (see Fig. 4c) and would thus influence the SHG effect from the corresponding orbitals. Clearly, in the P-phase the accumulated charges are mainly localized near the Se atoms, while in the I-phase these charges are located away from the Se atoms, forming a porous charge area in the middle of two neighboring layers. Figure 4c also shows that the interlayer shared electrons in the I-phase serve as a bridge between Se atoms in adjacent layers, making possibly a direct hybridization and charge transfer via the close-contacted (3.92 Å) interlayer Se…Se pair, i.e., the interlayer Z…Z pair, which is absent in the P-phase where the same Se atoms are separated by a longer distance (4.54 Å) and hardly has any shared charge under the AA stacking type. The interlayer Se…Se overlaps in the I-phase can be regarded as a known covalent-like quasi-bond which has also been observed in other 2D materials^36,37,38,39. This type of quasi-bonds weakens the localized Se-4p orbitals that have negative ${\chi }_{123}^{\left(2\right)}$ contribution, eventually resulting in substantial P/I SHG divergence.

Based on the above analysis, we schematically summarize the valence orbital distribution in layered A₂MZ₄ materials in a 3-region model as illustrated in Eq. 4d and Eq. (2):

$${P}_{\rm{total}}={P}_{[{{\rm{AZ}}_{4}}]}-{P}_{[{{\rm{MZ}}_{4}}]}+{P}_{{\rm{Z}}-{\rm{p}}}\pm \Delta {P}_{{\rm{inter}}}$$

(2)

where ${P}_{{\rm{total}}}$ is the overall SH polarization of a layered A₂MZ₄ material, ${P}_{{\rm{Z}}-{\rm{p}}}$, ${P}_{{[{\rm{AZ}}}_{4}]}$ and ${P}_{{[{\rm{MZ}}}_{4}]}$ are the SHG polarizations from the localized Z-p electrons, [AZ₄] and [MZ₄] motifs, respectively, corresponding to the major SHG-active bands from region I, II, and III. The minus sign indicates that the polarizations from [AZ₄] and [MZ₄] motifs counteract each other. Moreover, the interlayer polarization $\Delta {P}_{{\rm{inter}}}$ is added to describe the SHG change on the localized Z-p orbitals brought by the interlayer charge coupling as discussed above. In addition to Cd₂GeSe₄, the similar 3-region electronic structures and SHG distribution are also observed in other type-I, II and III layered A₂MZ₄ materials, such as Cu₂MoS₄ (Supplementary Fig. 1), Cd₂SiS₄ (Supplementary Fig. 2) and Al₂CdS₄ (Supplementary Fig. 3), validating our theoretical model.

Overall, a larger distance change in the interlayer Z…Z pair between two phases can lead to greater $\Delta {P}_{{\rm{inter}}}$ difference and is beneficial to increase the SHG deviation D (Supplementary Fig. 4). Our further investigation reveals that the generally larger D value in the type-II layered A₂MZ₄ materials than in other derivatives can be ascribed to their highly distorted crystal structures and consequently enlarged Z…Z distance (${d}_{{\rm{z}}\ldots {\rm{z}}}$) change between P- and I-phases. Owing to the larger atom radius of A-site than that of M-site elements (Supplementary Fig. 5), type-II structures show a more prominent structural mismatch between [AZ₄] and [MZ₄] tetrahedrons (e.g., bond angles and bond lengths in Supplementary Table 1 and Supplementary Table 2) compared with typical type-I/III structures. Specifically, the [AZ₄] tetrahedra in all type-II structures are stretched compared to the case of typical type-I/III structures, which is illustrated in Fig. 5. Such distortion in structure then enlarges the location change of Z-site atoms along both horizontal and vertical directions via interlayer sliding which dominate the P/I SHG deviation in type-II structures.

Fig. 5: Illustration of lattice distortion in type-II derivatives.

a Side view. b Top view. c Oblique view. Z1 (red), Z2’ (green), and Z3(pink) represent different Z-site atoms in AA/AB-stacked phase. Z1 and Z3 are in the upper layer and Z2’ in the lower layer. The black dashed line represents the interlayer Z…Z pair which is composed of Z1 and Z2’ (Z3 and Z2’) in AA (AB) stacking. ${d}_{{\rm{z}}\ldots {\rm{z}}}$ represents the distance between the paired Z atoms in adjacent layers. ${d}_{\parallel }$ and ${d}_{\perp }$ are its horizontal and vertical components.

Looking into the XY plane, the Z-site atoms in most type-I and III structures are located near 1/4 (and 3/4) of the face diagonal in a unit cell (Supplementary Table 3). As a result, the Z3 site in AB stacking almost coincides with the original Z1 site in AA stacking, meaning that the interlayer Z1…Z2’ pair in AA stacking and Z3…Z2’ pair in AB stacking are nearly equivalent and the horizontal ${d}_{\parallel }$ component of ${d}_{{\rm{z}}\ldots {\rm{z}}}$ barely changes as shown in Fig. 5a. In all type-II structures, however, the Z atoms deviate from their original position and draw closer to the M site because of the distorted structure, and these two interlayer Z…Z pairs before and after the (1/2,1/2,0) interlayer shift show apparently smaller ${d}_{\parallel }$ distance in AB stacking. At the same time, along the c-axis the vertical ${d}_{\perp }$ components of ${d}_{{\rm{z}}\ldots {\rm{z}}}$ (directly related to the layer distance h in Table 1) are substantially reduced between P- and I-phases in the type-II structures compared with the type-I/III structures. Also, the structural mismatch in the type-II structures results in a flattened lattice, as depicted in Fig. 5b. Both make the type-II structures significantly reduced ${d}_{\perp }$ for I-phases. Therefore, the type-II structures undergo dramatic distance and orientation change for the interlayer Z atom pair (${{d}_{{\rm{z}}\ldots {\rm{z}}}}^{2}={d}_{\parallel }^{2}+{d}_{\perp }^{2}$) in the layer sliding process (Fig. 5c), which makes them exhibit prominent P/I SHG deviation compared with other type structures.

The detailed structural data, including the fractional coordinates in the XY plane for Z1, Z2’, and Z3 site of P/I-A₂MZ₄ and corresponding ${d}_{{\rm{z}}\ldots {\rm{z}}}$, ${d}_{\parallel }$ and ${d}_{\perp }$, are listed in Supplementary Table 3. The above analysis is also supported by the CDD comparison among type-I, II, and III compounds, e.g., Cu₂MoS₄, Cd₂SiS₄, and Al₂CdS₄ (see Supplementary Fig. 6 in supporting information).

Stability and feasibility

With the extremely large D value, it deserves attention to evaluate the stability of Cd₂GeSe₄. First-principles calculations show that the layered Cd₂GeSe₄ has no negative vibration mode in phonon spectra (Supplementary Fig. 7 in supplement materials), indicating its dynamic stability and the potential to be practically synthesized. Also, no prominent difference in E_g between P and I phases is found (Table 1 and Supplementary Fig. 8), meaning that its optical absorption and transmission properties remain stable during the P/I phase change by (1/2, 1/2, 0) layer shift.

In addition, the feasibility of actual phase change is essential for the potential application of Cd₂GeSe₄. By transition state search, the energy barrier is calculated to be 78.87 meV per atom for the P-to-I (AA stacking-to-AB stacking) process driven by interlayer sliding and 110.13 meV per atom for the reverse process, as illustrated in Fig. 6a, which is on the same magnitude as for phase transformation in other bulk vdW systems^40,41, making it possible to be actively switched between these two phases given certain external excitations (e.g., temperature and stress). In a bilayer system, remarkably, the interlayer sliding is allowed under ambient conditions since the required energy barriers are significantly reduced (15.72 or 36.73 meV per atom for the same AA-to-AB or AB-to-AA processes, according to Fig. 6b) and comparable to the thermal excitation energy at room temperature (~26 meV per atom which is roughly estimated by E=k_bT when T = 300 K)¹⁶, despite slightly larger than those in some common 2D vdW materials like h-BN and MoS₂^42,43,44. It is also practical to obtain monolayer or multilayer AA-/AB-Cd₂GeSe₄ by mechanical exfoliation from bulk P/I phases as the calculated exfoliation energies (Fig. 6c) for ≤10 layers are around 10 ~ 20 meV per Ų, on the same level with other representative 2D materials⁴⁵. All these exfoliated finite-layer Cd₂GeSe₄ present a remarkable layer dependence in SHG effect as the disparity between low/high-SHG states enlarges with the increase of layer numbers, see Fig. 6d. Eventually, the ${\chi }_{123}^{\left(2\right)}$ of AA- (AB-) stacked multilayers would converge at approximately −8.7 pm·V⁻¹ ( ~ 0.2 pm·V⁻¹) which is roughly equivalent to that of bulk P (I) phase. Overall, the stacked Cd₂GeSe₄ may serve as an easily tunable SHG-active layered material with a large SHG modulation ratio of >36 times and small exfoliation energy. Our study reveals that, although the interlayer vdW interaction is weak, it has a significant influence on the SHG generation. This finding is also consistent with other recent studies^46,47 in which the interlayer charge transfer in vdW heterostructures can lead to the emergence of SHG signal. Our study stands out for introducing large SHG regulation without the break of inversion symmetry.

Fig. 6: Evaluation on the feasibility to actively tune the NLO property of stacked Cd₂GeSe₄ phases by sliding and to obtain few-layer Cd₂GeSe₄ from bulk phases.

a Evolution of total energy from P- to I-phase along the transition path for bulk Cd₂GeSe₄. b Evolution of total energy from AA to AB stacking along the transition path for bilayer Cd₂GeSe₄. c Exfoliation energies of monolayer and multilayer Cd₂GeSe₄ with AA/AB stacking from corresponding bulk (P/I) phases. d Layer-dependent SHG effect in Cd₂GeSe₄. All finite-layer Cd₂GeSe₄ systems are fully geometrical optimized.

In conclusion, we systematically investigated the tunable SHG effect of A₂MZ₄-type ternary van der Waals layered chalcogenides which can be regulated by P/I phase change featuring a (1/2, 1/2, 0) interlayer shift without changing the non-centrosymmetric structure. Among all 34 investigated materials, Cd₂GeSe₄ has an extremely large regulation ratio in SHG strength (>36 times between high- and low-SHG states), making it a strong contender for practical materials with tunable and even switchable NLO properties. Our theoretical analysis proposes that the occupied SHG-active bands in layered Cd₂GeSe₄ can be divided into 3 regions, and the SHG strength change is rooted in the near-edge region I dominated by Se 4p orbitals. We then find that the different interlayer Se…Se charge coupling on localized Se 4p orbitals, although resulting from the weak vdW interaction, is crucial for the strong SHG deviation in the two stacking types of Cd₂GeSe₄. The AB- stacked I-phase is more favorable for the formation of interlayer Se…Se coupling, weakening the localization of Se 4p orbitals and consequently influencing their SHG contribution. Such a regulation model can further be extended to other substituted A₂MZ₄-type layered chalcogenides with different elemental compositions, introducing a mechanism to adjust optical nonlinearity for vdW layered materials by tuning interlayer charge redistribution.

Methods

First principle calculations

In the first-principles calculations, the energy band gap E_g, density of states (DOS), and static second order NLO susceptibility ${\chi }_{{ijk}}^{(2)}$ of all designed materials are calculated using CASTEP⁴⁸, a software package based on the density functional theory (DFT). The generalized gradient approximation (GGA)⁴⁹ of Perdew−Burke−Ernzerhof (PBE)⁵⁰ functional is chosen and k-point meshes spanning less than 0.04/ų in the Brillouin zone are generated by Monkhorst-Pack method⁵¹. Norm-conserving pseudopotentials (Cu: 3d¹⁰4s¹, Ag: 4d¹⁰5s¹, Mo: 4d⁵5s¹, W: 5d⁴6s², Zn: 3d¹⁰4s², Cd: 4d¹⁰5s², Hg: 5d¹⁰6s², Al:3s²3p¹, Ga: 3d¹⁰4s²4p¹, In: 5s²5p¹ Si: 3s²3p², Ge: 4s²4p², S: 3s²3p⁴, Se: 4s²4p⁴) are used to describe the interactions between the valence electrons and the ionic cores⁵². A dispersion correction by Tkatchenko and Scheffler is adopted to account for the interlayer vdW interactions⁵³. Other parameters are set as default in the ultra-fine option of the software. All structures are optimized using the BFGS method⁵⁴. In optical calculation, a Monkhorst-Pack k-point mesh of 5 × 5 × 5 (5 × 5 × 2) which spans less than 0.04/ų is chosen for the P- (I-) phase. In type-I materials, the number of empty bands used to calculate SHG susceptibility is 200 (400) for P- (I-) phase. In type-II and III materials, the number of empty bands used to calculate SHG susceptibility is 200 for both phases. Convergence test for k-point and empty bands are performed on Cu₂MoS₄ and Cd₂SiS₄ as shown in Supplementary Fig. 12. Due to the underestimation of energy band gap E_g by PBE functional, the HSE06 functional^55,56 is used to obtain the accurate E_g and the scissors operator⁵⁷ is adopted to correct the optical calculations. For type-II materials with large SHG deviations, phonon calculation by linear response method is performed to evaluate their dynamic stabilities⁵⁸.

Calculation of NLO properties

The static SHG susceptibility ${\chi }_{{ijk}}^{(2)}$ is calculated using an expression developed by Lin et al.⁵⁹:

$${{\rm{\chi }}}_{{ijk}}^{(2)}={{\rm{\chi }}}_{{ijk}}^{(2)}\left(\text{VE}\right)+{{\rm{\chi }}}_{{ijk}}^{(2)}\left(\text{VH}\right)+{{\rm{\chi }}}_{{ijk}}^{(2)}({\rm{two}}\, {\rm{bands}})$$

(3)

$${{{\chi }}}_{{ijk}}^{(2)}\left(\text{VH}\right)=\frac{{{\rm{e}}}^{3}}{2{{{\hslash }}}^{2}{m}^{3}}\sum _{v{v}^{{\prime} }c}\int \frac{{{\rm{d}}}^{3}k}{{4{{\pi }}}^{3}}P\left({ijk}\right)\text{Im}[{p}_{v{v}^{{\prime} }}^{i}{p}_{{v}^{{\prime} }c}^{j}{p}_{{cv}}^{k}]\left(\frac{1}{{{{\omega }}}_{{cv}}^{3}{{{\omega }}}_{v{\prime} c}^{2}}+\frac{2}{{{{\omega }}}_{{vc}}^{4}{{{\omega }}}_{{cv}{\prime} }}\right)$$

(4)

$${{{\chi }}}_{{ijk}}^{(2)}\left(\text{VE}\right)=\frac{{{\rm{e}}}^{3}}{2{{{\hslash }}}^{2}{m}^{3}}\sum _{{vcc}{\prime} }\int \frac{{{\rm{d}}}^{3}k}{{4{{\pi }}}^{3}}P\left({ijk}\right)\text{Im}[{p}_{{vc}}^{i}{p}_{{{cc}}^{{\prime} }}^{j}{p}_{{c}^{{\prime} }v}^{k}]\left(\frac{1}{{{{\omega }}}_{{cv}}^{3}{{{\omega }}}_{{vc}{\prime} }^{2}}+\frac{2}{{{{\omega }}}_{{vc}}^{4}{{{\omega }}}_{c{\prime} v}}\right)$$

(5)

$${{{\chi }}}_{{ijk}}^{(2)}\left({{\rm{two}}\, {\rm{bands}}}\right)=\frac{{{\rm{e}}}^{3}}{2\,{{{\hslash }}}^{2}{m}^{3}}\sum _{{vc}}\int \frac{{{\rm{d}}}^{3}k}{{4{{\pi }}}^{3}}P\left({ijk}\right)\,\frac{{\text{Im}}[{p}_{{vc}}^{i}{p}_{{cv}}^{j}({p}_{{vv}}^{k}-{p}_{{cc}}^{k})]}{{{{\omega }}}_{{vc}}^{5}}$$

(6)

Here, i, j and k are Cartesian components. v and v’ denote valence bands, while c and c’ denote conduction bands. $P\left({ijk}\right)$ denotes full permutation and shows the Kleinman symmetry⁶⁰ of the SHG coefficients. The entire SHG coefficient ${\chi }_{{ijk}}^{(2)}$ of a material can be divided into contributions from virtual hole (VH), virtual electron (VE), and two-band processes. The contribution from two-band process is extremely small and can be neglected. The band-resolved SHG analysis is performed to visualize the contribution of the ${\chi }_{{ijk}}^{(2)}$ in a level-by-level basis so that the states giving major contributions to the SH process can be identified⁶¹. The energy-dependent SHG susceptibility is calculated based on similar mechanisms as described by Duan et al.⁶².

Orbital-selective ${\chi }^{(2)}$-weighted electron density schemes are obtained by multiplying the ${\chi }^{(2)}$ contribution and electron density of certain orbitals. With this method, the charge density of selected orbitals irrelevant to SHG are not shown whereas those vital to SHG are displayed in real space. If the orbital electron density $\rho \left({\bf{r}}\right)$ of all j states is expressed as Eq. (7), then the ${\chi }^{(2)}$-weighted electron density ${\rho }_{\chi }\left({\bf{r}}\right)$ of all i states can be expressed as Eq. (8). Here ${\bf{r}}$ represents the position vector and $\varPsi$ represents the wave function. ${\chi }_{i}^{(2)}$ is the ${\chi }^{(2)}$ contribution of a certain i state.

$$\rho \left({\bf{r}}\right)=\sum _{j}{\left|{\varPsi }_{j}\left({\bf{r}}\right)\right|}^{2}$$

(7)

$${\rho }_{\chi }\left({\bf{r}}\right)=\sum _{i}{\chi }_{i}^{(2)}{\left|{\varPsi }_{i}\left({\bf{r}}\right)\right|}^{2}$$

(8)