Ultra-low lattice thermal conductivity and high thermoelectric performance in chiral phonon-protected heterostructures

Abstract

Heterostructures have recently emerged as suitable candidates for energy applications owing to their intriguing high electronic properties and low thermal conductivity. Traditionally, experimentally observable chiral phonons, protected by the crystal hexagonal symmetry in two-dimensional honeycomb lattices and originating primarily at the Brillouin zone corner, play a critical role in achieving ultra-low lattice thermal conductivity (κ_l). It is mainly attributed to the three-fold rotational symmetry-driven circularly polarized chiral phonons that restrict certain scattering processes and reduce thermal resistance. Herein, we identify an innovative mechanism that not only preserves phonon chirality but also minimises lattice thermal conductivity by satisfying certain phonon selection rules and enhancing anharmonic scattering, particularly in the long wavelength limit. This study encompasses first-principles calculations and solutions to the Boltzmann transport equation, along with various thermal transport theories to investigate the structural, thermal and electronic transport properties of group-IV-based dichalcogenides in hexagonal phase MSe₂ (M = Mo and W) and their heterostructures (HS), i.e. MoSe₂/WSe₂ and MoSeTe/WSeTe. Interestingly, phenomena such as membrane effect and hybridization of acoustics and low-lying optical phonons combined with phonon branching and phonon bunching are identified to execute an imperative role in suppressing thermal transport, resulting in a low κ_l value of 2.83 W m⁻¹ K⁻¹ in MoSe₂/WSe₂ and an ultra-low κ_l of 0.5 W m⁻¹ K⁻¹ in MoSeTe/WSeTe HS. The disparity in computing carrier mobility using deformation potential theory (DPT) is addressed herein, and to rectify the inconsistency, we additionally incorporate the Fröhlich interaction, which accounts for longitudinal optical phonons in estimating the carrier mobility accurately. Overall, the present work exploits the intrinsic properties of chiral phonons through breaking the crystal inversion symmetry, providing important insights into the rational design of low-dimensional thermoelectric materials.

---

1. Introduction

The intrinsic chirality of phonons in non-magnetic two-dimensional (2D) materials exhibiting a hexagonal lattice with broken inversion symmetry was predicted theoretically at the corners of the Brillouin zone (BZ).¹ This occurs because phonons in such lattices inherit the threefold rotational symmetry (C₃) of the lattices and their function of motion is invariant under C₃ except for a phase difference, Ĉ₃(u_k) = e^−i(2π/3)lu_k, where u_k is the function of atomic motion and l with possible values of 0, ±1.^1,2 Particularly, phonons at K and K′ valleys are notable due to the absence of time-reversal symmetry, where both their lattice and momentum vectors are preserved.^1–3 Due to the C₃ symmetry of these vectors, the vibrational plane wave is composed of non-degenerate chiral phonons that are not due to the superposition of linear modes; in contrast, it is because of the circular polarization in nonchiral media.^2,4 Herein, the honeycomb structure of the 2D transition metal dichalcogenide (TMD) compound MSe₂ (M = Mo and W) also exhibited a hexagonal symmetry, with atomic arrangement governed by the space group P6̄m2 (D¹_3h).⁵ In these monolayers, each M atom is strongly bonded with six Se atoms, among which three Se atoms are in the upper layer and three Se atoms are in the lower layer, forming a trigonal prism with M at the central position. This particular atomic arrangement results in a single Se–M–Se layer with intralayer covalent bonding, which provides mechanical strength. Moreover, the construction of heterostructures (HS) by vertically stacking different 2D crystals on top of each other via interlayer van der Waals (vdW) interactions is one of the effective ways to break the inversion symmetry compared to their homo-structures such as MoSe₂/MoSe₂ and WSe₂/WSe₂, which can induce phonon chirality in the structure. Synthesising such materials not only induces chirality but also introduces new distinct features such as the emergence of low-lying optical (LLO) phonons,^6,7 phonon bunching,⁸ strong hybridisation of acoustic and optical phonons⁹ and phonon branching (dispersive nature of phonon branches),¹⁰ which can significantly reduce κ_l.

Recently, the emergence of chiral phonons in the WSe₂ monolayer has been identified experimentally through the intervalley scattering of holes during the indirect infrared absorption using transient infrared spectroscopy.¹ Moreover, the breaking of spatial inversion symmetry in MoSe₂/WSe₂ vdW heterostructures induces long-lived indirect interlayer excitons with opposite circular polarizations, as observed under the chiral optical pumping.^11,12 In recent times, Zhang et al. have carried out first-principles calculations to show the potential application of chiral phonons in the optical transition of MoS₂/WS₂ heterostructures.³ Moreover, Pandey et al. shows that chiral phonons in the 1D chain of Ba₃N follow new symmetry-based scattering rules that limit three-phonon interactions in the system.¹³ However, weak chain coupling breaks these scattering rules in the bulk structure, resulting in a lower κ_l value.¹³ Therefore, tailoring the mechanism of chiral phonons specifically by weakening phonon chirality can enhance the three-phonon scattering events. The inherent phonon–phonon interactions arising from such anharmonic systems can undergo the lowest-order anharmonic processes, namely AAA (acoustic–acoustic–acoustic), AAO (acoustic–acoustic–optical), AOO (acoustic–optical–optical) and OOO (optical–optical–optical), which includes the decay of a phonon into two others or the coalescence of two phonons into a third as reported by Ravichandran et al.¹⁴ Recently, the MoS₂/MoSe₂ bilayer heterostructure has shown an intrinsically low lattice thermal conductivity of 25.39 W m⁻¹ K⁻¹ compared to the values of 84 W m⁻¹ K⁻¹ and 59 W m⁻¹ K⁻¹ of their monolayer counterparts MoS₂ and MoSe₂, respectively, highlighting the effective significance of heterostructures in lowering the lattice thermal conductivity and their potential application in thermoelectric devices.¹⁵

Considering the experimental synthesis and technical effectiveness, we used the first-principles-based density functional theory (DFT) calculations to effectively model the MoSe₂/WSe₂ HS by choosing single-layer MoSe₂ and WSe₂ as parent materials. Additionally, we investigated the variation in thermoelectric properties occurring due to the substitution of Te atoms in place of Se, at one side of each layer in the HS. Furthermore, we explored the underlying mechanisms of phonon modes, which enable overcoming the limitation imposed by phonon chirality in minimising κ_l and maximising ZT in both HS compared to their monolayer counterparts. Our results reveal that the MoSe₂/WSe₂ HS exhibits phonon branching, leading to the AOO scattering process, whereas the bunching of acoustic and LLO phonon modes observed in the MoSeTe/WSeTe HS results in the AAA process. Interestingly, the MoSe₂/WSe₂ HS shows chirality due to the breaking of spatial inversion symmetry compared to their homo-structures; however, doping of Te induces distinctive phonon features such as phonon bunching and mass-induced phonon softening, enhancing the three-phonon scattering process resulting in an ultra-low κ_l value of 0.5 W m⁻¹ K⁻¹ at room temperature and a maximum ZT value of 3.37 at 700 K.

2. Computational details

In this study, the first-principles-based density functional theory and density functional perturbation theory (DFPT)¹⁶ calculations were used, as implemented in the Quantum Espresso (QE) package.¹⁷ The projector-augmented wave (PAW)¹⁸ pseudopotentials along with generalized gradient approximation (GGA),¹⁹ parameterized by the Perdew, Burke, and Ernzerhof (PBE)²⁰ method, were used to consider the exchange–correlation interactions to determine the structural and transport properties of both heterostructures. The van der Waals density functional scheme²¹ was used to describe the long-range interlayer interactions of the heterostructures. We consider a large vacuum gap of thickness 20 Å along the z-axis to eliminate the van der Waals (vdW) interactions between periodic layers. For both heterostructures, we used an optimized Monkhorst–Pack k-grid of 12 × 12 × 4 to sample the Brillouin zone (BZ), along with the kinetic energy cutoff of wavefunction and charge density on a plane-wave basis at 60 Ry and 240 Ry, respectively, throughout the calculation to ensure the precision and convergence in the optimized geometry. The convergence threshold criterion for self-consistency in these structures was maintained as 10⁻⁸ Ry. The ab initio determination of phonon transport properties is enabled by obtaining the interatomic force constants (IFCs) and iteratively solving the phonon Boltzmann transport equation (pBTE) as implemented in ShengBTE²² alongside the QE package.¹⁷ The phonon dispersion curves and harmonic (second order) force constants were computed by performing the DFPT calculations using 3 × 3 × 2 coarse grids of q-points with a strict convergence threshold in the self-consistent field (SCF) calculation. The third-order anharmonic IFCs were obtained using the thirdorder.py code,²² which utilises a cut-off to include interactions up to the third nearest neighbours. For both heterostructures, the third-order anharmonic interatomic force constants (IFCs) required for ShengBTE calculations were based on a 3 × 3 × 2 supercell.

In this study, for electronic transport calculations such as density of states (DOS) analysis for both heterostructures, we performed a non-self-consistent field (NSCF) calculation using an increased k-grid of 24 × 24 × 4. Moreover, we incorporated the LOBSTER code²³ as a post-processing tool along with the QE package for evaluating the crystal orbital Hamilton population (COHP)²⁴ analysis and computed the integrated crystal orbital bond index (ICOBI).²⁵ The COHP analysis revealed the interaction between two atomic orbitals of neighbouring atoms of a material in terms of their Hamiltonian matrix element.²⁴ The multiplication of this matrix element with the corresponding density of states gives a quantitative measure of the bonding strength arising from bonding, non-bonding, or antibonding interactions.²⁴ Additionally, the ICOBI values serve as an effective tool for understanding the chemical bonding hierarchy by detecting covalent and ionic bonds in solid-state materials.²⁵ Furthermore, we evaluated the electrical transport coefficient (i.e., Seebeck coefficient), by solving the semi-classical Boltzmann transport theory within the constant scattering time approximation (CSTA), as implemented in the BoltzTraP code.²⁶

3. Results and discussion

3.1. Crystal structure and stability

The honeycomb structure of the two-dimensional (2D) transition metal dichalcogenide (TMD) compound MSe₂, where M = Mo and W, exhibited a hexagonal symmetry with atomic arrangement governed by the space group P6̄m2 (D¹_3h).²⁷ In these monolayers, each M atom is strongly bonded with six Se atoms, among which 3 Se atoms are in the upper layer and 3 Se atoms are in the lower layer, forming a trigonal prism, where M is at the central position. This particular atomic arrangement results in a single Se–M–Se layer with intralayer covalent bonding, which provides mechanical strength. Additionally, the monolayer MX₂ compounds in their 1H phase preserve the mirror symmetry and have a non-centrosymmetric structure due to the lack of inversion symmetry.²⁸ The 1H phase of MSe₂ was synthesized experimentally,^29,30 which consisted of three atoms in the unit cell, where the metal atom was sandwiched between two layers of chalcogen atoms, revealed as a trigonal prismatic coordination.^27,31,32 The corresponding lattice parameter and bond length values of single-layer MoSe₂ and WSe₂ are presented in Table S1 of the ESI,† consistent with the previously reported theoretical values.²⁷ The negligible lattice mismatch of 0.1% is in the allowable error range, computed using ³³ where a₁ and a₂ are the lattice parameters of the monolayers MoSe₂ and WSe₂, respectively, demonstrating that these monolayers are suitable to construct a heterostructure (HS).³³

As shown in Fig. 1a, the MoSe₂ monolayer is vertically stacked on top of the WSe₂ monolayer via interlayer van der Waals interactions, constituting the MoSe₂/WSe₂ HS. Composing single-layer MoSe₂ on WSe₂ along the normal direction breaks the mirror symmetry at the ±K point, which gives rise to non-degenerate phonon modes with definite polarisation. Based on the relative position of the two monolayers, different stacking configurations are considered. Previous studies have demonstrated that the low-energy configuration of such heterostructures is achieved by stacking the metal atom of one layer directly above the chalcogen atoms of the other layer, namely the AA′ stacking.^34,35 In our present study, the Mo and Se atoms of the MoSe₂ layer were exactly aligned over the Se and W atoms of the WSe₂ layer, respectively, to form the MoSe₂/WSe₂ heterostructure with the AA′ configuration, as shown in Fig. 1a.³⁵ This heterostructure also lacked the inversion symmetry similar to its parent materials (MoSe₂ and WSe₂ monolayers), arising due to the AA′ configuration.³⁵ After geometry optimisation, the MoSe₂/WSe₂ HS maintains the hexagonal symmetry similar to its parent materials, with lattice parameters of a = b = 3.32 Å and an interlayer distance (d) of 6.49 Å (the vertical distance from Mo to W atom), which closely aligns with the previously reported theoretical value (d = 6.48 Å).³⁶ The bond lengths between Mo–Se (∼2.54 Å) and W–Se (∼2.54 Å) atoms in the heterostructure are identical and are closely similar to those in the individual MoSe₂ (∼2.54 Å) and WSe₂ (∼2.55 Å) monolayers,²⁷ respectively, as presented in Table S1 of the ESI.† The side and top views of the MoSe₂/WSe₂ HS are illustrated in Fig. 1a and b, respectively.

Fig. 1 Side view of the optimised heterostructures: (a) MoSe₂/WSe₂ and (c) MoSeTe/WSeTe. (b) and (d) Top view of MoSe₂/WSe₂ and tilted top view of MoSeTe/WSeTe heterostructures, respectively. Here, Se1 and Se2 correspond to different Se atoms connected to Mo and W in the MoSe₂/WSe₂ heterostructure.

Furthermore, we modeled the 2H phase of the MoSeTe/WSeTe HS by stacking one layer of Janus MoSeTe on top of the WSeTe monolayer, which preserved the hexagonal symmetry having AA′ stacking analogous to the MoSe₂/WSe₂ heterostructure. In our present study, we designed different configurational spaces with random sampling of Se and Te atoms to explore different atomic arrangements and calculate the DFT energies for those configurations. Our analysis reveals that the structure shown in Fig. 1c represents the most dynamically stable configuration under the harmonic approximation, which is opted for further investigation. The optimized MoSeTe/WSeTe heterostructure has a lattice parameter of a = b = 3.46 Å and an interlayer distance of 6.33 Å. This decrement in the interlayer distance of Te doped system is attributed to the enhanced van der Waals attraction between the sub-layers in MoSeTe/WSeTe (0.39 eV) compared to the MoSe₂/WSe₂ (0.35 eV) heterostructure. Additionally, the bond length of the Mo–Se atom is the same as the bond length of W–Se (∼2.56 Å), while the bond lengths of Mo–Te and W–Te are also found to be identical (∼2.73 Å), and the corresponding values are presented in Table S1 of the ESI.† The increase in lattice parameter and bond lengths in MoSeTe/WSeTe compared to the MoSe₂/WSe₂ HS is attributed to the larger atomic size of Te than that of the Se atom. The schematic representations of the side and tilted top views of the MoSeTe/WSeTe HS are depicted in Fig. 1c and d, respectively. Moreover, to assess the stability and bonding nature of these heterostructures, we evaluated the binding energy, Born criteria of elastic stability³⁷ (Table S2, ESI†), and performed the crystal orbital Hamilton population (COHP) (Fig. S1, ESI†) and Bader charge analyses (Table S3, ESI†), as discussed in ESI† (Note S1).

3.2. Phonon spectra

Over the years, extensive investigations have been carried out on lattice dynamics studies in understanding phonon dispersion, specifically phonon eigen modes at high symmetry points using the DFPT method. The computed phonon spectra of MoSe₂ and WSe₂ monolayers (parent materials) are shown in Fig. 2a and b, which are consistent with previously reported data.^32,38,39Fig. 2c and d illustrate the calculated phonon dispersion curves as a function of frequency for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, respectively, along the high-symmetry Γ–M–K–Γ direction in the irreducible Brillouin zone (BZ). A system has N atoms per unit cell, and the phonon band structure features 3N branches, including three acoustic branches and 3N − 3 optical branches.⁴⁰ Consequently, the HS consists of six atoms in the unit cell, yielding a total of 18 vibrational modes in the phonon dispersion. Importantly, none of the phonon dispersions show imaginary frequencies under the harmonic approximation, indicating the dynamical stability of both heterostructures. Among these 18 vibrational modes, the three lowest energy phonon modes correspond to the acoustic branches, while the remaining modes are classified as optical. The low energy acoustic branches include the ZA mode, characterized by an out-of-plane displacement and a quadratic dispersive behaviour, described as ω ∝ q², whereas the TA and LA modes correspond to in-plane transverse and longitudinal displacements, respectively, exhibiting a linear q dependence given as ω ∝ q.³⁸

Fig. 2 Phonon dispersion bands for (a) MoSe₂ and (b) WSe₂ monolayers, (c) MoSe₂/WSe₂ and (d) MoSeTe/WSeTe, along the high symmetry points in the BZ. Acoustic phonon, low lying optical (LLO) and high energy optical phonon branches are represented in red, green and blue colours, respectively.

Since phonon modes play a crucial role in determining the lattice thermal conductivity of semiconductors, we additionally conducted a detailed analysis of these phonon dispersion curves, shown in Fig. 2c and d, in comparison to their monolayer counterparts. For MoSe₂ and WSe₂ monolayers, there are three atoms in the unit cell, yielding a total of 9 vibrational modes in the phonon spectra. In the case of MoSe₂ monolayer, the phonon modes are arranged in three sets, where the first group are the three highly dispersive acoustic branches within a frequency range of ∼0–152 cm⁻¹, the second group are the three mid-range optical phonon branches confined within the frequency width of ∼160–233 cm⁻¹ and third are the three higher energy phonon bands within ∼275–344 cm⁻¹, as evident in Fig. 2a. However, in the WSe₂ monolayer, the presence of heavier W atoms shifts the three aforementioned higher energy phonon modes toward the mid-range optical region (∼238–298 cm⁻¹), thus resulting in the softening of these phonon modes. Interestingly, the phonon spectra of the WSe₂ monolayer exhibit bunching phonon modes (non-dispersive or branches coming close to one another), specifically observed for the acoustic branches, as shown in Fig. 2b, attributed to the presence of heavier atoms. The bunching of acoustic phonon branches eliminates the dispersive nature of acoustic modes and constrains them to a narrow frequency width of ∼0–120 cm⁻¹, which is 25% lower than that of the MoSe₂ monolayer.

Furthermore, our study reveals that the modelling of HS from their monolayers introduces new distinctive features in the phonon spectra. A notable characteristic in the phonon dispersion of heterostructures compared to its monolayers is the appearance of three low-lying optical (LLO)⁴¹ phonons (as evident in Fig. 2c and d), including the LO mode, with frequencies as low as 36.98 cm⁻¹ for MoSe₂/WSe₂ and 35.60 cm⁻¹ for MoSeTe/WSeTe heterostructures. Recently, Xiaokun Gu et al. have also reported the occurrence of LLO phonon branches in the MoS₂ bilayer, with the flexural ZO mode reaching a minimum of ∼60 cm⁻¹, which have direct evidence in reducing the lattice thermal conductivity by 12% compared to its single layer.⁴² The emergence of such LLO phonons in heterostructures can be attributed to the interlayer van der Waals interactions and the proximity effect.^33,43 The existence and impact of proximity effect in heterostructures demonstrate that the presence of one layer induces effects in another layer, which is absent in the case of monolayer. It has been previously investigated and shown to have significant influence on properties such as superconductivity and magnetism.⁴³ In our present study, the emergence of LLO phonons, including the LO mode whose frequency is found to be ∼76% lower than that of its parent materials at the high-symmetry Γ point (160 cm⁻¹ and 167 cm⁻¹ for MoSe₂ and WSe₂ monolayers, respectively, as shown in Fig. 2a and b), can be attributed to the proximity of the MoSe₂ sub-layer to the WSe₂ sub-layer. Furthermore, we plotted the eigen vectors of low-lying flexural ZO mode for both heterostructures, as shown in Fig. S2a and b (ESI†), where it is clearly visible that the aforementioned proximity effect limits the vibration of the WSe₂ sub-layer (similar to the WSeTe sub-layer in MoSeTe/WSeTe), thereby reducing the phonon frequency in both heterostructures compared to their parent materials.

Additionally, these LLO phonons in the MoSe₂/WSe₂ HS exhibit coupling with the neighbouring acoustic branches along the high symmetry points, as shown in the phonon band structure (Fig. 2c). Particularly, the LLO phonon branches, especially the LO mode in MoSe₂/WSe₂, exhibit significant dispersiveness (phonon branching) reaching a maximum frequency of 163 cm⁻¹. The presence of such dispersive phonon modes, specifically low-energy phonons, resembles monolayer MoSe₂ phonon band dispersion, attributed to the reduced frequency gap between the LLO and higher energy optical phonons, leading to an enhanced coupling between the optical branches compared to the single-layer MoSe₂ and WSe₂ structures. It is well known that the pronounced lattice thermal conductivity observed in MoSe₂ and WSe₂ monolayers is primarily due to the presence of a large frequency gap in the phonon dispersion curve.³⁸ However, the presence of dispersive LLO optical modes in the MoSe₂/WSe₂ HS reduces the frequency gap, contributing to an increased phonon scattering rate, which can simultaneously result in a low lattice thermal conductivity relative to its parent materials. Moreover, the occurrence of LLO modes in the low frequency range induces the acoustic branches of MoSe₂/WSe₂ HS to exhibit strong phonon softening (Fig. 2c) compared to its parent materials, as evident from the phonon dispersion curve shown in Fig. 2a and b, which further suppresses the heat carrying phonons in the system.

However, in the case of MoSeTe/WSeTe HS, the substitution of the heavier chalcogen atom Te in place of lighter Se reduces the dispersiveness of the phonon branches, effectively eliminating the phonon branching phenomenon, as evident in Fig. 2d. Moreover, the acoustic and LLO branches of the MoSeTe/WSeTe HS exhibit the phonon bunching phenomena similar to those observed in the phonon spectra of the WSe₂ monolayer (Fig. 2b). Additionally, the incorporation of Te atoms lowers the frequency (phonon softening) of all phonon modes, including the acoustic and LLO branches introducing distinct energy separation among the phonon groups. One such separation occurs between the LLO (∼30–110 cm⁻¹) and mid-range optical phonons (136–199 cm⁻¹), and another between the mid-range and high-energy optical phonons (216–311 cm⁻¹), attributed to the substitution of heavier Te atoms. The incorporation of heavier Te atoms follows the general diatomic linear chain model, where the zone edge frequency of the acoustic branch is ω_ac ∝ M^−1/2 (M is the mass of heavier element), which causes the relative shift of acoustic spectral frequencies to a lower value in MoSeTe/WSeTe than that in the MoSe₂/WSe₂ HS, which is consistent with the results obtained by Thomas et al.⁴¹ Moreover, the frequency gap between the LLO and mid-range optical modes of the MoSeTe/WSeTe HS primarily arises from the minimal dispersiveness of the LLO branches compared to the MoSe₂/WSe₂ HS. Notably, the bunching nature of acoustic and LLO phonons in the MoSeTe/WSeTe HS results in strong coupling between the acoustic and LLO phonon branches, confining them within a narrower frequency width of ∼108 cm⁻¹, which is lower than that of the MoSe₂/WSe₂ HS (∼161 cm⁻¹). This enhanced coupling between the acoustic phonons themselves and with LLO phonon branches may enhance the phonon scattering mechanisms, which limits the lattice thermal conductivity in the system. Furthermore, in the MoSeTe/WSeTe HS, the coupling of phonon branches varies with the q direction and is found to be stiffer at the high symmetry K point (as evident in Fig. 2d), attributed to the occurrence of the phonon softening phenomena. Henceforth, the strong coupling between acoustic and LLO phonons, combined with pronounced acoustic and optical phonon softening, will significantly limit the group velocity of heat carrying acoustic phonons and enhance the anharmonic three-phonon scattering process in the MoSeTe/WSeTe HS, resulting in a low lattice thermal conductivity in the structure.

To provide a more comprehensive understanding of the atomic contribution, we also plotted the phonon density of states (phonon DOS) as a function of frequency, as illustrated in Fig. S3 of the ESI.† In the case of MoSe₂ monolayer, low-frequency phonon modes are mainly contributed by the vibration of both Mo and Se atoms, whereas higher frequencies are dominated by the vibration arising from the Mo atom. However, in the case of WSe₂ monolayer, the low-frequency range is dominated by the vibration of heavier W atom and higher frequencies are dominated by the vibration from the Se atoms. For the MoSe₂/WSe₂ HS, it is clearly visible that the major proportion of atomic vibrations in the low-frequency region encompassing both the acoustic and LLO branches are primarily arising from the heavier W atom. However, the higher energy optical branches are due to the vibration of lighter Mo and Se atoms. Furthermore, in the case MoSeTe/WSeTe HS, the major contribution to the acoustic region is from the heavier Te atoms along with vibrations arising from the W atom. The high-energy optical modes are due to the rattling atoms Mo and Se, which are found to have a large atomic vibration as observed in the high energy range in the phonon spectra.

3.3. Phonon chirality

Phonon chirality,^1,2,44 defined by the circular polarisation of phonons (S^z_ph), is a distinct feature that emerges at the high symmetry points in the BZ, particularly in systems with a two-dimensional (2D) honeycomb or hexagonal lattice, which exhibits a three-fold rotational symmetry (C₃).¹ In such honeycomb lattices, phonons remain invariant under three-fold discrete rotations about the axis perpendicular to the lattice plane at high symmetry Γ and ±K points.² Moreover, modelling a material with the lack of spatial inversion symmetry is crucial for obtaining non-degenerate phonon modes with definite circular polarisation at a high-symmetry K valley. Henceforth, in the case of MoSe₂/WSe₂ HS, the aforementioned lack of spatial inversion symmetry arises due to the presence of two layers with differing atomic masses (i.e., one layer composed of Mo and the other of W) within the unit cell. A recent study has demonstrated that, when two layers are composed of identical masses (homo-structures such as MoSe₂/MoSe₂ and WSe₂/WSe₂), the spatial inversion symmetry is maintained, resulting in opposite phonon polarizations in the two layers, which cancel out each other, leading to zero net polarisation at the BZ corner.⁴ However, in such cases, the superposition of degenerate modes reveal the circular polarisation at ±K valleys.⁴ Furthermore, by breaking a spatial inversion symmetry through the introduction of layers with different masses, as in the MoSe₂/WSe₂ HS compared to their homo-structures, we can preserve a definite circular phonon polarization in one of the sub-layers, although, the other sub-layer remains static, as evident in Fig. 3, resulting in a net phonon chirality at valley points ±K. The strength of these non-vanishing circular polarisations S^z_ph can be measured in terms of their quantized phonon pseudo angular momentum (PAM), where, S^z_ph > 0 and S^z_ph < 0 correspond to right-handed and left-handed circular polarization of the sub-lattice or heterostructure, respectively.³ Furthermore, S^z_ph = 0 represents the linear vibration of the sub-lattice, resulting in the absence of chirality.³

Fig. 3 Lattice vibrations of phonon modes at the high-symmetry K point for the MoSe₂/WSe₂ HS and the corresponding frequencies of each mode from 1 to 18. Circular polarisation or phonon chirality is represented by ellipses⁴⁹ with red arrows. Forces along the out-of-plane direction (z axis) are represented by purple arrows. The size of each ellipse and the length of each purple arrow are proportional to the amplitudes of the circular and linear vibrations, respectively. The direction of each red arrow denotes the corresponding direction of chirality and the phonon modes are arranged in the increasing order of energy.

Furthermore, the conservation of such quantized phonon PAM contributes an additional selection rule apart from the well-known energy and momentum conservation rules. Hence, in our present study, the introduction of three-fold rotation operator to the phonon wave function ψ_ph yields the expression where l_ph gives the phonon PAM, which has a discrete value of 0, +1, or −1 at the ±K points with a three-fold rotational symmetry.³ Herein, all vibrating sub-layers or sub-lattices contribute to ψ_ph and their vibrations give the same phase change under three-fold rotation, i.e., the same l_ph. Henceforth, we can compute l_ph by any vibrating sub-lattice. In contrast, the sub-lattices that remain static do not have any contribution to ψ_ph and thus cannot be used to obtain l_ph. Interestingly, the term ψ_ph has both a nonlocal contribution from the Bloch phase factor e^iq·r and a local contribution from sub-lattice relative vibrations.³ Henceforth, usually the phonon PAM is obtained as a sum of corresponding nonlocal orbital PAM (l^o) and local spin PAM (l^s) of the vibrating sub-lattice.³ Moreover, l^o depends on atomic in-plane sites and l^s is determined by sub-lattice vibration. The right-handed circular vibration, left-handed circular vibration, and linear vibration along the z-axis are all eigen states of with l^s = 1, −1, and 0, respectively.³ In order to show the conservation of phonon PAM, we also computed the aforementioned right-handed circular vibration, left-handed circular vibration, and linear vibration along the z-axis for both the heterostructures, as evident in Fig. 3 and 4. Here, the circular vibrations are shown in ellipse and the direction of rotation (right and left) is indicated by the red arrow. A recent study has suggested that the occurrence of such phonon chirality limits the three-phonon interaction occurring in the system, attributed to a higher lattice thermal conductivity (κ_l).¹³

Fig. 4 Lattice vibrations of phonon modes at the high-symmetry K point for the MoSeTe/WSeTe heterostructure and the corresponding frequencies for each mode from 1 to 18.

In the MoSe₂/WSe₂ and MoSeTe/WSeTe HS, the time reversal symmetry at the Γ point is preserved, resulting in the vanishing of phonon polarisation (nearly equal to zero).³ Hence, we extend our analysis to the K point, where the phonon has definite circular polarization arising from the non-degenerate chiral phonon modes due to the breaking of time reversal symmetry.³ To provide a comprehensive understanding of how phonon chirality affects the κ_l value of MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, we calculated the phonon eigen vectors and analysed the specific vibrational motion of both layers, as illustrated in Fig. 3 and 4. The phonon chirality on individual atoms is represented by an ellipse, with the amplitude of vibrations indicated by the radius of the ellipse. The different chirality of phonons representing the right- and left-hand phonon rotations is depicted by arrows within the circles. The phonon eigen vectors for each mode at the high-symmetry K point in the MoSe₂/WSe₂ heterostructure closely resemble those of the previously reported MoS₂/WS₂ bilayer heterostructure.³ Moreover, the strengths of phonon chirality determined from the phonon PAM (l_ph) by incorporating the terms l^o and l^s for both heterostructures are presented in Tables 1 and 2. In the case of MoSe₂/WSe₂ HS, a majority of the phonon modes possess a definite circular polarisation indicated by l_ph = ±1, as evident from Table 1. Interestingly, these chiral phonons observed in the MoSe₂/WSe₂ heterostructure differ from those in the MoSeTe/WSeTe heterostructure, probably due to the substitution of Te atoms. As shown in Fig. 3, it is evident that in the case of MoSe₂/WSe₂ HS, either the MoSe₂ or the WSe₂ layer remains stiff with no circular polarisation, while the net chirality arises from other layer, or vice versa, attributed to the lack of spatial inversion symmetry compared to their homo-structures. In contrast, the incorporation of Te atoms in the MoSeTe/WSeTe HS induces opposite circular polarisation in both layers, which will cancel each other and result in a zero net phonon chirality. Furthermore, comparing the phonon PAM values (shown in Tables 1 and 2), it is clear that in the case of MoSe₂/WSe₂ and MoSeTe/WSeTe HS, the phonon chirality remains persistent. Recent investigation by Tribhuwan Pandey et al. has shown that the symmetry-driven chirality in single-chain Ba₃N restricts the phonon–phonon interactions, resulting in high lattice thermal conductivity.¹³ However, the weakening of this chirality by breaking the symmetry in its bulk counterpart enhances the phonon scattering, leading to a threefold reduction in κ_l compared to the single-chain structure.¹³ In contrast, for MoSe₂/WSe₂ and MoSeTe/WSeTe HS, despite the presence of phonon chirality, the effective tuning of phonon features such as emergence of LLO phonon modes, phonon branching and phonon bunching (as discussed earlier) induced by heterostructure formation may significantly enhance the phonon–phonon scattering, probably resulting in low κ_l.

Table 1 PAM corresponding to each phonon mode and phonon frequency at the high-symmetry K point for the MoSe₂/WSe₂ HS. Here, l^o and l^s represent the nonlocal orbital and local spin PAM, respectively. The static sub-lattices have no well-defined l^s, denoted by a dash (—). l_ph corresponds to the phonon PAM of the particular phonon mode, considering both l^o and l^s of the sub-lattices

Modes	Frequency (cm^−1)	l ^oMo,Se1,Se2	l ^oW,Se3,Se4	l ^sMo	l ^sSe1	l ^sSe2	l ^sW	l ^sSe3	l ^sSe4	l ph
1	83	−1	1	—	—	—	1	0	0	−1
2	113.9	−1	1	—	—	—	0	−1	−1	1
3	127.8	−1	1	—	—	—	−1	1	1	0
4	129.7	−1	1	0	1	1	—	—	—	−1
5	132.4	−1	1	−1	0	0	—	—	—	1
6	157.5	−1	1	1	−1	−1	—	—	—	0
7	189.9	−1	1	—	0	0	—	—	—	0
8	195.3	−1	1	—	−1	−1	—	—	—	0
9	195.4	−1	1	—	1	1	—	0	0	−1
10	200.3	−1	1	—	1	1	—	0	0	−1
11	204.7	−1	1	—	—	—	—	0	0	0
12	205.9	−1	1	—	—	—	—	1	1	0
13	240.5	−1	1	—	—	—	−1	1	1	0
14	240.9	−1	1	—	—	—	0	−1	−1	1
15	250.2	−1	1	—	—	—	−1	0	0	1
16	278.1	−1	1	−1	0	0	—	—	—	1
17	280.9	−1	1	0	1	1	—	—	—	−1
18	300.7	−1	1	1	−1	−1	—	—	—	0

---

Table 2 PAM corresponding to each phonon mode and phonon frequency at the high-symmetry K point for the MoSeTe/WSeTe HS

Modes	Frequency (cm^−1)	l ^oMo,Se,Te	l ^oW,Se,Te	l ^sMo	l ^sSe	l ^sTe	l ^sW	l ^sSe	l ^sTe	l ph
1	47.4	−1	1	0	—	1	—	—	−1	0
2	48.8	−1	1	—	—	1	0	—	−1	0
3	52.9	−1	1	—	—	—	−1	0	0	1
4	68.2	−1	1	1	0	0	—	—	—	−1
5	69	−1	1	−1	—	—	0	0	1	0
6	72.3	−1	1	0	0	−1	—	—	—	1
7	146.8	−1	1	—	—	—	1	0	0	−1
8	153.1	−1	1	−1	0	—	—	—	—	0
9	159.9	−1	1	—	1	—	−1	1	1	0
10	165.9	−1	1	—	—	—	—	−1	−1	0
11	166.5	−1	1	—	1	1	—	—	—	0
12	174.1	−1	1	0	−1	—	—	—	—	0
13	217.1	−1	1	—	—	—	0	1	—	0
14	220.9	−1	1	—	—	—	−1	1	—	0
15	226.4	−1	1	—	—	—	1	0	—	0
16	245.7	−1	1	−1	0	—	—	—	—	0
17	254.4	−1	1	0	—	—	—	—	—	0
18	254.4	−1	1	1	−1	—	—	—	—	0

---

3.4. Phonon transport properties

To provide a comprehensive understanding of the phonon–phonon interaction occurring in MoSe₂/WSe₂ and MoSeTe/WSeTe HS, we evaluated the anharmonic scattering rates as a function of frequency for both heterostructures at room temperature, as depicted in Fig. 5a and b. Moreover, the anharmonic three-phonon scattering process is crucial for understanding the phonon dynamics during thermal transport and can be calculated using the following expression:⁴⁰where λ, λ′, and λ′′ denote the three phonon branches involved in the collision process, Γ⁺ stands for the absorption process given by λ + λ′ → λ′′ and Γ⁻ represents the emission process given as λ → λ′ + λ′′. Moreover, all the three-phonon scattering processes satisfy the conservation of momentum and energy given by q ± q′ = q′′ + G, where G is the reciprocal lattice vector, and ω_λ ± ω_λ′ = ω_λ′′, respectively.

Fig. 5 (a) and (b) Phonon scattering rate, (c) and (d) mode-dependent Grüneisen parameters and (e) and (f) group velocity for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, respectively, at 300 K.

The phonon frequency range shown in Fig. 5 is consistent with the dispersion curves as in Fig. 2. It is noteworthy that the scattering rate of the MoSe₂/WSe₂ HS shown in Fig. 5a is much higher than that of the two parent materials (0.1 ps⁻¹ and 1 ps⁻¹ for MoSe₂^38,45 and WSe₂⁴⁵ monolayer, respectively), within the whole frequency range. Previous studies have shown that the presence of high anharmonic scattering rates is one of the fundamental reasons for low lattice thermal conductivity in 1T phase monolayers of HfSe₂, SnSe₂ and ZrS₂.^40,46 Additionally, the appearance of LLO phonons is the causative factor for the enhanced anharmonic scattering process observed in the heterostructures compared to their parent materials. In the case of MoSe₂/WSe₂ HS, it is notable that the scattering rate of LLO phonons (∼12.5 ps⁻¹ at a frequency range of 150 cm⁻¹) is more pronounced than acoustic and optical phonon modes, attributed to the phonon branching phenomena. Such dispersive LLO phonon branches get scattered by both acoustic and high-energy optical phonons due to the reduced frequency gap between them, resulting in a higher scattering rate of LLO modes. Hence, these highly dispersive phonon branches dominate three-phonon scattering, which primarily arises from the AOO process⁴⁷ involving low-energy acoustic and optical modes within the frequency range of 150 cm⁻¹. Moreover, the large optical bandwidth of the MoSe₂/WSe₂ HS further strengthens the AOO process to dominate the scattering rates even at low frequencies.

Furthermore, in the MoSeTe/WSeTe HS, the anharmonic scattering rate is higher across the entire frequency range, as evident from the phonon eigen modes shown in Fig. 4. In the MoSe₂/WSe₂ HS, the presence of chirality resulting from the breaking of spatial inversion symmetry compared to their homo-structure counterparts constrains the anharmonic interactions and reduces the phonon scattering rate in the system. Moreover, the doping of Te induces distinctive phonon phenomena (phonon bunching) in conjunction with mass anisotropy, which will substantially enhance the three-phonon scattering events that can consequently reduce the lattice thermal conductivity in MoSeTe/WSeTe compared to the MoSe₂/WSe₂ HS. The doping of such heavier Te atom induces substantial phonon softening, particularly affecting the acoustic and LLO phonon modes. This mass-induced phonon softening constraints the phonon frequency towards a low-energy region (discussed earlier in phonon spectra using a diatomic linear chain model), facilitating enhanced phonon–phonon interactions and, thereby, increasing the phonon scattering rate. Additionally, in the low frequency range below 35 cm⁻¹, the MoSeTe/WSeTe HS exhibits a pronounced phonon scattering rate of ∼105 ps⁻¹, attributed to the enhanced coupling of acoustic phonons themselves and with LLO modes, resulting from the bunching behaviour of phonon branches and strong acoustic phonon softening, due to the incorporation of Te atoms. Our results show that the bunching nature (non-dispersive) of acoustic phonon branches, combined with strong acoustic phonon softening within a very low frequency range of 75 cm⁻¹, induces phonon–phonon interactions, primarily dominated by the AAA process.⁴⁷ The anharmonic scattering values observed in the MoSeTe/WSeTe HS are higher than those previously reported for the MoS₂–MoSe₂ lateral superlattice, which shows a promising lattice thermal conductivity of ∼14.6 W m⁻¹ K⁻¹ among the 2H phase TMDs.³⁸

To further understand the increased anharmonic scattering rate occurring in MoSeTe/WSeTe compared to the MoSe₂/WSe₂ heterostructure (as discussed above), we performed a microscopic analysis on weighted phase space (W_3ph), a robust factor that describes the allowed phonon–phonon scattering limited by fundamental momentum and energy conservation conditions. Fig. S4 (ESI†) illustrates the weighted phase space as a function of frequency, depicting the available transition probabilities for phonon absorption (λ + λ′ → λ′′) and emission (λ → λ′ + λ′′) processes in both heterostructures. The transition probabilities for phonon–phonon scattering reach maximum values of 0.001 ps⁴ rad⁻⁴ (at 20 cm⁻¹) and 0.08 ps⁴ rad⁻⁴ (at 10 cm⁻¹) for MoSe₂/WSe₂ and MoSeTe/WSeTe HS, respectively. Interestingly, the W_3ph plot for the MoSeTe/WSeTe HS clearly highlights the presence of significantly higher transition probabilities or scattering processes depicted as more dense points, particularly in the low-frequency region below 75 cm⁻¹, attributed to the bunching of acoustic and LLO phonon branches. This bunching nature of phonon modes greatly increases the available phase space for three-phonon scattering processes, particularly the AAA scattering channels. The proximity of phonon branches constrained within a narrow frequency range enables a higher probability of scattering events in the structure, leading to a substantial increase in the anharmonic scattering rate. The presence of significant scattering channels combined with high transition probability values observed in the MoSeTe/WSeTe HS further substantiates the computed high anharmonic scattering rate in the system compared to the MoSe₂/WSe₂ HS. Moreover, the dominant scattering channels available in the MoSeTe/WSeTe HS can substantially suppress the heat carrying phonons, further leading to a low lattice thermal conductivity.

Additionally, the anharmonic scattering rate of the investigated heterostructures shows prominent scattering events in the acoustic and LLO phonon frequency ranges, indicated by highly dense points for both HS, as evident in Fig. 5a and b, compared to the previously reported WSe₂ and WSe₂/WTe₂ bulk structures⁴⁸ as shown in Fig. S5 of the ESI.† These intensified scattering events observed in MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures suggest enhanced phonon–phonon scattering, which suppresses the heat transport and probably contributes to a significantly low κ_l value.

In order to explain the curious behaviour of lattice thermal conductivity in the studied heterostructures, we evaluated the strength of anharmonic interactions occurring in the system by computing the mode-specific Grüneisen parameter (γ) expressed as follows:⁴⁹

where ω_k(q) is the frequency of the phonon mode k, q is the wave vector, and V is the volume. Fig. 5c and d illustrate the mode Grüneisen parameters for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures at 300 K, respectively. For the MoSe₂/WSe₂ HS, the mode Grüneisen parameter in the long wavelength limit reaches a maximum of 33, which is higher than that of the previously reported MoS₂/MoSe₂ bilayer heterostructure.³³ Similar high γ values are distinctive features observable in the 1T phase structures such as SnS₂ (∼30)⁴⁶ and SnSe₂ (∼35)⁴⁶ monolayers, which serve as fundamental factors responsible for suppressing their lattice thermal conductivity as low as 6.4 W m⁻¹ K⁻¹ and 3.8 W m⁻¹ K⁻¹, respectively.⁴⁶ Moreover, the average γ value given as⁵⁰ of MoSe₂/WSe₂ HS was found to be 3.49. Such high Grüneisen parameter observed in the MoSe₂/WSe₂ HS can probably contribute to a very low lattice thermal conductivity. However, the substitution of Te increases the strength of anharmonic interaction in the MoSeTe/WSeTe HS, as depicted in Fig. 5d, compared to the MoSe₂/WSe₂ heterostructure. Herein, the enhanced γ value in the Te-doped system manifests the elevated anharmonicity in the crystal structure, which stems from the phonon softening, phonon bunching and enhanced coupling of acoustic and LLO branches (as discussed earlier). This significant increase in the magnitude of anharmonic scattering strength due to Te doping results in an elevated scattering rate observed in the MoSeTe/WSeTe HS, as illustrated in Fig. 5b. Furthermore, the negative value of γ observed for the acoustic branches of the MoSeTe/WSeTe HS caused upon Te doping has already been shown in various 2D structures,³³ and is primarily driven by the membrane effect,^33,51 along with the hybridisation of acoustic branches with LLO phonons.³³ The phenomena of membrane effect arising due to larger atomic radius associated with Te doping are consistent with the observed Grüneisen values in the MoSeTe/WSeTe HS, suggesting that the expansion of the membrane (increase in lattice parameter) and chemical bond length always tends to increase the compressibility of the phonon modes, especially in the low-frequency region (TA and LA) except for the ZA branch. This enhanced compressibility increases anharmonic interactions at long wavelengths, leading to higher γ values in heavier Te-doped MoSeTe/WSeTe than in MoSe₂/WSe₂. Furthermore, the exceptional behaviour of the ZA branch is attributed to the usual fact that, when the membrane undergoes expansion, the atoms in the layer are less free to move along the out-of-plane direction (z axis).

To uncover the thermal transport mechanisms in both heterostructures, we also computed the phonon group velocity of each mode expressed as³³ where ω, k and q denote the phonon frequency, vibrational mode index and wave vector, respectively. Fig. 5e and f illustrate the frequency-dependent phonon group velocities at 300 K for both MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, respectively. Notably, the MoSe₂/WSe₂ HS possesses the highest phonon group velocity across the entire frequency region compared to the MoSeTe/WSeTe HS, in which a maximum group velocity of 8 km s⁻¹ is observed for the LA mode, while the TA and quadratic ZA modes reach 4.1 km s⁻¹ and 2.5 km s⁻¹, respectively. Moreover, the LLO modes display high group velocities (7.8 km s⁻¹), reaching as high as the longitudinal acoustic mode, arising due to their dispersive nature (as shown in Fig. 2c), allowing the LLO phonons to travel longer distances without scattering. The low group velocities of higher energy optical phonons indicate the dominating nature of acoustic and LLO phonon contribution to the lattice thermal conductivity. However, Te doping suppresses the phonon group velocities of all vibrational modes in the low-energy acoustic region due to increased scattering, visible in Fig. 5f, which results from the strong coupling nature of the acoustic branches themselves as well as with the LLO phonon branches and the enhanced phonon softening in the MoSeTe/WSeTe HS. Furthermore, the incorporation of Te in the MoSeTe/WSeTe HS reduces the frequency of LLO modes (especially, the LO mode as discussed earlier), arising due to phonon softening, which limits the group velocity of the LA branch to a minimum of ∼6 km s⁻¹. Such LLO phonon softening also enhances the coupling between these branches, leading to a low group velocity (5.7 km s⁻¹) for these modes. Consequently, the suppression of group velocity via Te substitution may significantly lower the κ_l value.

The temperature-dependent lattice thermal conductivity (κ_l) was calculated by solving the phonon Boltzmann transport equation (pBTE) in conjunction with the interatomic force constants (IFCs), as implemented in the ShengBTE package.²² The lattice thermal conductivity is given as follows:²²

where N, Ω, f₀ and υ_λ are the number of uniformly spaced q points in the Brillouin zone, volume of the unit cell, Bose–Einstein distribution function depending on the phonon frequency ω_λ, and phonon group velocity, respectively. F^β_λ = τ^β_λ(υ_λ + Δ_λ), where τ^β_λ stands for the phonon lifetime under the single-mode relaxation time approximation (SMRTA) and Δ_λ is the deviated value predicted by SMRTA.

Fig. 6a and b illustrate the lattice thermal conductivity as a function of temperature, calculated using both relaxation time approximation (RTA) and iterative methods. Our results reveal a very low lattice thermal conductivity for the MoSe₂/WSe₂ HS, with a minimum of ∼2.83 W m⁻¹ K⁻¹ using the iterative method and ∼2.42 W m⁻¹ K⁻¹ by RTA, whereas the MoSeTe/WSeTe HS shows an ultra-low value of 0.50 W m⁻¹ K⁻¹ using the iterative method and 0.142 W m⁻¹ K⁻¹ by RTA at room temperature. There is discrepancy in the κ_l values obtained from both methods, as RTA neglects the term Δ_λ, because the approximation considers the excitation of each individual phonon mode ignoring the phonon distribution.⁴⁶ Such assumption works very well only for materials where normal scattering dominates the Umklapp scattering, which is typically applicable for high-thermal conductivity materials.⁴⁶ In contrast, the iterative method, which accounts for the three-phonon scattering, solves the pBTE exactly by including both the Umklapp and normal scattering processes and provides a more accurately converged κ_l value, which is applicable for both high- and low-lattice thermal conductivity materials.⁴⁶ At high temperatures, the mismatch of RTA values in predicting κ_l is less than 3%, as the Umklapp scattering becomes more significant. The computed lattice thermal conductivity values for both heterostructures are significantly lower than that of single-layer MoSe₂ (∼43 W m⁻¹ K⁻¹)^36,38,52 and WSe₂ (∼52 W m⁻¹ K⁻¹)^36,52 at room temperature. A previous study has shown that the creation of the MoS₂–MoSe₂ lateral superlattice increases the anharmonic scattering in the system, resulting in a promising κ_l value of 14.6 W m⁻¹ K⁻¹ at 300 K, which is nearly 1/5 and 1/4 of the two parent materials.³⁸ Similarly, the modelling of the MoSe₂/WSe₂ HS has reduced the lattice thermal conductivity by ∼98% and ∼94% compared to the MoSe₂ and WSe₂ monolayer counterparts, respectively. The strategy of stacking MoSe₂ over the WSe₂ monolayer reduces the acoustic–optical phonon band gap (key factor for higher κ_l in the parent materials as discussed earlier), significantly enhances the three-phonon scattering dominated by the AOO process and consequently helps to achieve a low κ_l value of ∼1.24 W m⁻¹ K⁻¹ at 700 K (as shown in Fig. 6a), similar to those in the 1T phase of the HfSe₂ and SnSe₂ monolayers.^40,46 Moreover, the substitution of Te in the MoSeTe/WSeTe HS amplifies the three-phonon scattering process along with strong anharmonic interactions arising due to the mass anisotropy and distinct phonon features such as phonon softening and phonon bunching in the system, leading to an ultra-low κ_l value of 0.5 W m⁻¹ K⁻¹ at room temperature. The low lattice thermal conductivity achieved by incorporating heavy Te atoms follows the traditional Keyes theory,⁵³ which describes inverse relationship between κ_l and heavier elements. The presence of strongly coupled acoustic and LLO phonons in the low-frequency region, along with the high Grüneisen parameter, further increases the scattering rate dominated by the AAA process in the MoSeTe/WSeTe HS, consequently resulting in a very low κ_l value reaching a minimum of 0.21 W m⁻¹ K⁻¹ at 700 K, which is ∼83% lower than that of the MoSe₂/WSe₂ HS.

Fig. 6 Lattice thermal conductivity using iterative, RTA and Slack method as a function of temperature for (a) MoSe₂/WSe₂ and (b) MoSeTe/WSeTe heterostructures, respectively. (c) and (d) Cumulative lattice thermal conductivity along the x and y directions as a function of frequency for (a) MoSe₂/WSe₂ and (b) MoSeTe/WSeTe heterostructures, respectively.

Furthermore, to provide a comprehensive understanding of the acoustic phonon contribution to κ_l in both heterostructures, we computed the lattice thermal conductivity using the Slack model, that assumes that the lattice thermal resistance arises only from the intrinsic phonon–phonon interaction, and is expressed as follows:^54,55

where A = 3.04 × 10⁻⁶ is the correction factor⁵⁴ for κ_s in units of W m⁻¹ K⁻¹, M_av is the average atomic mass, δ is the volume per atom in cubic angstrom, γ is the average Grüneisen parameter, n is the number of atoms in the primitive unit cell, T is the temperature and θ_D is the Debye temperature given as^46,55 where h, k_B and m are the Planck constant, Boltzmann constant and number of atoms per unit volume, respectively, and ν_a denotes the average group velocity given as⁵⁵. The computed κ_s values for both MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures are plotted as a function of temperature, as shown in Fig. 6a and b, respectively, and the corresponding values are given in Table S5 of the ESI.† For the MoSe₂/WSe₂ HS, the lattice thermal conductivity calculated using the Slack model shows higher values than those calculated by the iterative method across the entire temperature range. By considering solely the acoustic phonon interactions, κ_s for the MoSe₂/WSe₂ HS was found to be as high as 4.58 W m⁻¹ K⁻¹ at 300 K, indicating the significance of optical phonon contribution in the system. Moreover, as we discussed earlier, the three-phonon scattering in the MoSe₂/WSe₂ HS is primarily dominated by the AOO process; consequently, the higher κ_s value in the system is attributed to the negligence of the optical phonon interaction. For the MoSeTe/WSeTe HS, the Slack model predicts a comparatively low lattice thermal conductivity of 0.15 W m⁻¹ K⁻¹ at room temperature, indicating the inaccuracy of the Slack model in computing κ_s for the complex system. The complex crystal structure of the MoSeTe/WSeTe HS induces distinct phonon phenomena including bunching of LLO phonons and strong coupling of LLO phonons with acoustic modes as evident from the phonon spectra, as shown in Fig. 2d. However, the Slack model accounts for the lattice thermal conductivity of materials by considering the scattering arising from the acoustic branches. Consequently, for MoSeTe/WSeTe, the Slack model considers the three-phonon scattering mechanism solely dominated by the AAA processes, driven by the strong coupling of acoustic branches resulting from strong phonon softening, attributed to the substitution of the Te atom. The aforementioned AAA process in the MoSeTe/WSeTe HS leads to a high Grüneisen parameter (especially for the acoustic branches, as evident in Fig. 5d) as well as high Debye temperature (∼201.27 K) in the system, resulting in a very low κ_s value. A similar trend has been observed in a recent study, where the Slack model predicted a low lattice thermal conductivity value of ∼17.6 W m⁻¹ K⁻¹, compared to the iterative method (∼54 W m⁻¹ K⁻¹) for single-layer MoSe₂.⁵² Hence, such low value of κ_s observed for the MoSeTe/WSeTe HS indicates the significant domination of acoustic phonon interaction over the low lattice thermal conductivity.

Furthermore, to understand the specific phonon contribution towards thermal conductivity, we plotted the cumulative lattice thermal conductivity as a function frequency for both the heterostructures, as shown in Fig. 6c and d. In the case of MoSe₂/WSe₂ HS, κ_l is mainly contributed by phonons within the frequency range of <150 cm⁻¹, contributing 91.5% to the total lattice thermal conductivity. However, in the case of MoSeTe/WSeTe HS, the major contribution (88.2%) is constrained within a very low frequency range of <75 cm⁻¹, suggesting the overall significance of low-energy phonon contribution in lattice thermal conductivity.

3.5. Carrier transport properties

Apart from phonon transport properties, to accurately account for the electronic transport properties of the studied heterostructures, we plotted the electronic band structures as well as projected density of states (PDOS) for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, respectively, as illustrated in Fig. S6 (ESI†). Further analyses detailing the orbital contribution to electronic transport properties are elaborated in the ESI† (Note S2). Moreover, to account for the mobility of charge carriers, the deformation potential theory (DPT) provides a formula that incorporates carrier effective mass to accurately calculate the longitudinal acoustic (LA) phonon-limited mobility in 2D semiconductors given as follows:^40,56where μ_LA is the mobility of charge carriers, C_2D is the elastic constant, ħ is the reduced Planck constant, e is the electronic charge, E₁ is the deformation potential constant, m* is the carrier effective mass, k_B is the Boltzmann constant and T is the temperature. Detailed description of carrier transport formalism is provided in the ESI† (Note S3).

The calculated m*, C_2D and E₁ values for electrons and holes in MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures are provided in Table S6 (ESI†). The variations in curvature observed in the band structures of MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures result in distinct effective masses for charge carriers. In the case of MoSe₂/WSe₂ HS, the m* values of electrons and holes were found to be 0.22m₀ and 0.27m₀, respectively, where m₀ represents the bare electron mass. In the MoSeTe/WSeTe HS, the effective mass of electrons (∼0.23m₀) is similar to that of the MoSe₂/WSe₂ HS, while the effective mass of holes shows a steady decline of 0.20m₀. These results suggest that the substitution of Te instead of the Se atom does not significantly alter the curvature of CB, but induces slightly steeper bands in the valence band, contributing to a lower effective mass of holes in the MoSeTe/WSeTe HS. Table S6 (ESI†) presents the elastic constant values for both heterostructures, which are significantly higher than those reported for parent materials (MoSe₂ and WSe₂)⁵⁷ and the previously reported MoS₂/WS₂ HS (∼259.95 J m⁻²).⁵⁸ These higher C_2D values, as shown in Table S6 (ESI†), indicate that these heterostructures exhibit enhanced ability to withstand external or deforming forces. Furthermore, the deformation potential constant was calculated under both tensile and compressive strain conditions. Notably, the magnitude of hole–phonon coupling strength consistently surpasses the electron–phonon coupling strength in both heterostructures, indicating higher mobility for electrons than that for holes. In the case of MoSeTe/WSeTe HS, the substitution of Te in place of the Se atom slightly reduces the E₁ value for electrons (∼2.58 eV), suggesting a weaker acoustic phonon–electron coupling in MoSeTe/WSeTe than that in the MoSe₂/WSe₂ HS (∼2.71 eV). However, the Te substitution significantly increases the deformation potential constant (∼7.36 eV) for holes, suggesting a higher probability of scattering events to occur in MoSeTe/WSeTe than that in the MoSe₂/WSe₂ HS, which could potentially result in a lower hole mobility.

Fig. 7a and b illustrate the acoustic phonon limited mobility as a function of temperature for both heterostructures, and the corresponding values at room temperature are given in Table S6 (ESI†). The computed mobility values indicate that substituting Te for the Se atom results in a significant increase in the mobility of both electrons and holes in MoSeTe/WSeTe compared to the MoSe₂/WSe₂ heterostructure. Despite similar effective masses of electron in both the heterostructures, the pronounced increase in the electron mobility of MoSeTe/WSeTe (∼93 664.03 cm² V⁻¹ s⁻¹) compared to the MoSe₂/WSe₂ HS (∼49 157.20 cm² V⁻¹ s⁻¹) is attributed to elevated C_2D and E₁. The computed μ_LA value for both the heterostructures is higher than that of their constituent monolayers⁵⁷ and previously reported Janus 2D materials.⁵⁷ Moreover, to understand the underlying mechanism of higher mobility in heterostructures, we computed the relaxation time for both the HS using the equation derived as,⁵⁹ and the calculated relaxation time values for both MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures at 300 K are presented in Table S6 (ESI†). The observed relaxation time for both electrons (∼12.25 ps) and holes (∼1.73 ps) in MoSeTe/WSeTe are significantly higher than those of the MoSe₂/WSe₂ HS, indicating the reduced scattering events in the system, which allows the carriers to travel a longer distance without any hinderance, resulting in high carrier mobility in the MoSeTe/WSeTe HS.

Fig. 7 (a) p-Type (holes) and (b) n-type (electrons) LA limited mobility, (c) p-type and (d) n-type LO limited mobility and (e) p-type and (f) n-type total mobility as a function of temperature for both heterostructures.

In our study, the obtained DPT mobilities (Fig. 7a and b) are not quite promising, as only the scattering of charge carriers by LA phonons is considered and the contribution from optical phonons is completely neglected.⁵⁷ In such scenario, considering the scattering arising from longitudinal optical (LO) phonons is inevitable to precisely determine the total mobility of the materials. The LO mode-resolved phonon mobility is given as follows:⁵⁷

where M_M and M_X are the masses of the metal and chalcogen atoms, respectively. Moreover, A is the cross-sectional area of the unit cell, t is the effective thickness, ε is the optical dielectric constant, ħ is the reduced Planck constant, ω is the phonon frequency, e is the electronic charge, n_ħω is the Bose–Einstein distribution function, m* is the carrier effective mass and Z_MB is the Born effective charge of the metal atom. The μ_LO values as a function of temperature for both the MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures computed using eqn (6) are illustrated in Fig. 7c and d, respectively. The computed μ_LO values for both heterostructures are very low compared to μ_LA, emphasizing the fact that the mobility of these heterostructures is mostly dominated by LO phonon scattering. Interestingly, the computed μ_LO value of MoSe₂/WSe₂ is higher than that of the MoSeTe/WSeTe HS for both electrons and holes across the entire range of temperatures. All the parameters used for calculating the LO phonon limited mobility using the Fröhlich interaction are summarized in Table S7 of the ESI.†

Furthermore, the precise total mobility (μ_total) limited by both LA and LO phonons was computed using Matthiessen's rule derived as,⁵⁷ and Fig. 7e and f represent the total mobility as a function of temperature. The results for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures indicate that μ_total is primarily limited by the LO phonon scattering, consistent with previous studies on 2D TMDs.⁵⁷ The computed μ_total values of HS is very close to the LO phonon limited mobility, henceforth, the omission of LO phonon scattering is the pivotal factor in the overestimation of DPT mobility by nearly an order of magnitude two in MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, aligning with previously reported study.⁶⁰ Interestingly, although the MoSeTe/WSeTe HS exhibits higher μ_LA, the incorporation of scattering from LO phonons results in a lower μ_total value for electrons (∼221.61 cm² V⁻¹ s⁻¹) and holes (∼91.41 cm² V⁻¹ s⁻¹) in MoSeTe/WSeTe than that of the MoSe₂/WSe₂ HS. This indicates that the elevated LO phonon scattering impedes the charge carrier transport in the MoSeTe/WSeTe HS.

In our present study, we employed the versatile software tool BoltzTraP²⁶ to compute the temperature-dependent Seebeck coefficient (S). The Seebeck coefficient was obtained by solving the semi-classical Boltzmann transport equation within the constant scattering time approximation (CSTA), as implemented in the BoltzTraP code.²⁶ In this theory, the Seebeck coefficient is expressed as follows:²⁶

where α and β are the tensor indices, and Ω, μ, e, f₀ and T are the cell volume, chemical potential, electronic charge, Fermi distribution function, and temperature, respectively. Moreover, σ_αβ are the transport distribution tensor elements calculated using the Fourier interpolation of the electronic band structure.

The temperature-dependent Seebeck coefficients for p- and n-type MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures are shown in Fig. 8a and b, and the corresponding values are presented in Table S8 of the ESI.† Our findings reveal that the n-type Seebeck coefficient of the MoSe₂/WSe₂ HS (∼1316.53 μV K⁻¹) is significantly higher than that of the MoSeTe/WSeTe heterostructure (∼671.56 μV K⁻¹) at 300 K, attributed to the presence of heavy bands that result in an elevated carrier effective mass in the MoSe₂/WSe₂ HS. A similar pattern is observed for the p-type S values, where the MoSe₂/WSe₂ HS (∼1302.52 μV K⁻¹) exhibits a value notably higher than that of the MoSeTe/WSeTe HS (∼612.74 μV K⁻¹). The computed S values for the MoSe₂/WSe₂ HS are nearly twice as high as that of the MoSeTe/WSeTe structure, indicating a significantly higher thermopower. The obtained S values for the heterostructures are higher than those reported for the parent materials (MoSe₂ and WSe₂ monolayers).³² Notably, the S values for both p- and n-type heterostructures show a steady decline with the increase in temperature, which can be attenuated to the excitation of more carriers at higher temperatures. This trend of declining Seebeck values with the increase in temperature is consistent with the previously reported data.⁵

Fig. 8 Seebeck coefficient: (a) p-type (holes) and (b) n-type (electrons) as a function of temperatures using the BoltzTraP code for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures. (c) and (d) show computed electrical conductivity as a function of temperature for p-type and n-type for both HS.

Furthermore, we evaluated the electrical conductivity,⁵⁹σ = neμ, which is a transport coefficient property directly related to the mobility of charge carriers, to accurately compute the thermoelectric figure of merit, ZT. Here σ, n, e and μ represent the electrical conductivity, carrier concentration, charge of electron, and mobility of charge carriers, respectively. Fig. 8c and d illustrate the calculated electrical conductivity as a function of temperature for the p-type and n-type in the MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures. The optimisation of carrier concentration in conjunction with the Seebeck coefficient value plays a pivotal role in the fine-tuning of electrical conductivity to achieve more precise and accurate ZT values. In our study, we carefully optimised the carrier concentration in the range of 10¹⁷–10¹⁸ cm⁻³ and 10¹⁸–10¹⁹ cm⁻³ for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, respectively, aligning it with the previously reported data.³² The calculated electrical conductivity values show an increasing trend with temperature, as evident from Fig. 8c and d. The electrical conductivity values in the MoSeTe/WSeTe HS are significantly higher for both p-type (∼49.38 S cm⁻¹) and n-type (∼87.97 S cm⁻¹) than those of the MoSe₂/WSe₂ heterostructure at 700 K. The pronounced electrical conductivity observed for the MoSeTe/WSeTe HS over the entire temperature range is due to the enhanced carrier concentration, arising due to the low Seebeck coefficient in the system compared to the MoSe₂/WSe₂ HS.

3.6. Figure of merit

Using the phonon and electron transport coefficients calculated above, we computed the dimensionless figure of merit, ZT, for both p- and n-type semiconductors as a function of temperature in the MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures, as illustrated in Fig. 9a and b. The prevalent n-type nature of both the heterostructures is clearly visible with ZT values reaching a maximum of 0.72 for the MoSe₂/WSe₂ heterostructure and 3.37 for the MoSeTe/WSeTe heterostructure at 700 K, as evident in Fig. 9. The ZT value of the MoSeTe/WSeTe HS shows an increment across the entire temperature range, surpassing that of the MoSe₂/WSe₂ structure. Such pronounced ZT value in the MoSeTe/WSeTe structure is attributed to the overall low κ_l value in the system, arising due to distinct phonon features such as phonon bunching and mass-induced phonon softening combined with strong coupling of low-energy phonons. Moreover, such low κ_l in the MoSeTe/WSeTe heterostructure helps to overcome the limitation of low S, resulting in superior thermoelectric efficiency of the system. Furthermore, for the p-type MoSeTe/WSeTe heterostructure, the optimised ZT value (0.77) is obtained at 500 K, but shows a slight decline at high temperatures, attributed to the very low S. The incorporation of Te in the MoSeTe/WSeTe heterostructure enhances the ZT value at higher temperatures, outperforming the previously reported thermoelectric materials such as MoSe₂ (∼0.25),³² WSe₂ (∼0.8),³² and MoS₂/WS₂ heterostructure (∼0.65).⁵⁸ Such high ZT observed in the MoSeTe/WSeTe HS is primarily due to the emergence of low κ_l. In the MoSe₂/WSe₂ heterostructure, such low κ_l is due to the presence of highly dispersive LLO phonons, which reduce the acoustic–optical phonon band gap, whereas in the MoSeTe/WSeTe HS, the substitution of Te increases the three-phonon scattering process and enhances the magnitude of the anharmonic scattering strength due to mass anisotropy and phonon bunching, resulting in a low κ_l value and consequently high ZT.

Fig. 9 Thermoelectric figure of merit, ZT: (a) p-type (holes) and (b) n-type (electrons) as a function of temperature for MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures. (c) Thermoelectric efficiency (η) for both p-type and n-type heterostructures.

Furthermore, we also computed the thermoelectric efficiency of both heterostructures using the following expression:⁶¹

where the average ZT value is given as , where T_h and T_c denote the temperatures of the hot and cold sides of the thermoelectric device, respectively. For both the heterostructures, we considered T_h and T_c as the highest (700 K) and lowest (300 K) temperatures maintained for computing the ZT values. The n-type efficiency of the MoSeTe/WSeTe HS reaches a maximum of 20.3%, almost two times higher than the efficiency of the MoSe₂/WSe₂ heterostructure (10.1%), as illustrated in Fig. 9c. The prominent n-type nature of the materials compared to the p-type counterparts is revealed from their low efficiency, especially observed for the MoSe₂/WSe₂ HS (0.33%), whereas the MoSeTe/WSeTe HS exhibits an efficiency of 9.9%. The superior efficiency observed for the MoSeTe/WSeTe HS can be attributed to the ultra-low κ_l value in the material.

4. Conclusion

In summary, we have studied the thermal and electronic transport properties of MoSe₂/WSe₂ and MoSeTe/WSeTe heterostructures using the first-principles-based DFT calculations. Our results revealed that these heterostructures are dynamically, energetically, and mechanically stable based on their phonon band structure, binding energy and Born criteria of elastic constant, respectively. The low lattice thermal conductivity in both heterostructures compared to their monolayer counterpart is due to the distinctive features observed in the phonon spectra. The emergence of LLO phonon modes and strong coupling between the acoustic and LLO branches are the common features observed in the phonon spectra of both heterostructures. Additionally, the MoSe₂/WSe₂ HS exhibited the phonon branching phenomena, which led to the enhanced AOO scattering process, resulting in a low κ_l value compared to the parent material. In the case of Te substitution, strong phonon softening results in the bunching of acoustic and optical phonon modes creating a window for enhanced AAA scattering processes. Furthermore, the Te doping has increased the three-phonon scattering process in the system compared to the MoSe₂/WSe₂ HS, which is clearly analysed from the scattering rate and Grüneisen parameter arising due to distinctive phonon features and mass anisotropy in the respective system. Furthermore, the transition probability of three-phonon scattering channels assessed from the weighted phase space volume also shows an increased probability of scattering in the MoSeTe/WSeTe HS. All these distinct phonon dynamic features are responsible for the observed low κ_l value of 2.83 W m⁻¹ K⁻¹ in MoSe₂/WSe₂ and ultra-low κ_l of 0.5 W m⁻¹ K⁻¹ in the MoSeTe/WSeTe HS. Moreover, we have incorporated the LO phonon scattering by considering the Fröhlich interactions, which is an important factor for computing the accurate total mobility. By considering both phonon and electronic transport properties, n-type MoSeTe/WSeTe shows high ZT values of 3.37 (with an efficiency of 20.3%) and 0.72 (with an efficiency of 10.1%) for the MoSe₂/WSe₂ HS at 700 K revealing their prominent n-type character.

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

All data supporting this study are included in the article and its ESI.†

Acknowledgements

JP and RKB acknowledge SERB-DST (EEQ/2022/000325) for financial support.

References

1. H. Zhu, J. Yi, M. Y. Li, J. Xiao, L. Zhang, C. W. Yang, R. A. Kaindl, L. J. Li, Y. Wang and X. Zhang, Science, 2018, 359, 579–582.
2. L. Zhang and Q. Niu, Phys. Rev. Lett., 2015, 115, 115502.
3. W. Zhang, A. Srivastava, X. Li and L. Zhang, Phys. Rev. B, 2020, 102, 174301.
4. X. Xu, W. Zhang, J. Wang and L. Zhang, J. Phys.: Condens. Matter, 2018, 30, 225401.
5. G. Özbal, R. T. Senger, C. Sevik and H. Sevinçli, Phys. Rev. B, 2019, 100, 085415.
6. X. Gu, B. Li and R. Yang, J. Appl. Phys., 2016, 119, 085106.
7. J. J. Ma, J. J. Zheng, X. L. Zhu, P. F. Liu, W. D. Li and B. T. Wang, Phys. Chem. Chem. Phys., 2019, 21, 10442–10448.
8. A. Kundu, Y. Chen, X. Yang, F. Meng, J. Carrete, M. Kabir, G. K. H. Madsen and W. Li, Phys. Rev. Lett., 2024, 132, 116301.
9. W. Li, J. Carrete, G. K. Madsen and N. Mingo, Phys. Rev. B, 2016, 93, 205203.
10. M. J. Matthews, M. A. Pimenta, G. Dresselhaus, M. S. Dresselhaus and M. Endo, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, R6585.
11. B. Miller, A. Steinhoff, B. Pano, J. Klein, F. Jahnke, A. Holleitner and U. Wurstbauer, Nano Lett., 2017, 17, 5229–5237.
12. A. T. Hanbicki, H. J. Chuang, M. R. Rosenberger, C. S. Hellberg, S. V. Sivaram, K. M. McCreary, I. I. Mazin and B. T. Jonker, ACS Nano, 2018, 12, 4719–4726.
13. T. Pandey, C. A. Polanco, V. R. Cooper, D. S. Parker and L. Lindsay, Phys. Rev. B, 2018, 98, 241405.
14. N. K. Ravichandran and D. Broido, Phys. Rev. X, 2020, 10, 021063.
15. J. J. Ma, J. J. Zheng, X. L. Zhu, P. F. Liu, W. D. Li and B. T. Wang, Phys. Chem. Chem. Phys., 2019, 21, 10442–10448.
16. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515–562.
17. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, J. Phys.: Condens. Matter, 2009, 21, 395502.
18. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
19. X. Hua, X. Chen and W. A. Goddard, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 16103–16109.
20. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
21. J. Klimes, D. R. Bowler and A. Michaelides, J. Phys.: Condens. Matter, 2009, 22, 022201.
22. W. Li, J. Carrete, N. A. Katcho and N. Mingo, Comput. Phys. Commun., 2014, 185, 1747–1758.
23. S. Maintz, V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Comput. Chem., 2016, 37, 1030–1035.
24. V. L. Deringer, A. L. Tchougréeff and R. Dronskowski, J. Phys. Chem. A, 2011, 115, 5461–5466.
25. P. C. Müller, C. Ertural, J. Hempelmann and R. Dronskowski, J. Phys. Chem. C, 2021, 125, 7959–7970.
26. G. K. H. Madsen and D. J. Singh, Comput. Phys. Commun., 2006, 175, 67–71.
27. G. Özbal, R. T. Senger, C. Sevik and H. Sevinçli, Phys. Rev. B, 2019, 100, 085415.
28. P. Zheng, W. Wei, Z. Liang, B. Qin, J. Tian, J. Wang, R. Qiao, Y. Ren, J. Chen, C. Huang, X. Zhou, G. Zhang, Z. Tang, D. Yu, F. Ding, K. Liu and X. Xu, Nat. Commun., 2023, 14, 1–7.
29. E. Liu, J. van Baren, Z. Lu, T. Taniguchi, K. Watanabe, D. Smirnov, Y. C. Chang and C. H. Lui, Nat. Commun., 2021, 12, 6131.
30. A. Ramasubramaniam, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 115409.
31. D. Han, X. Yang, M. Du, G. Xin, J. Zhang, X. Wang and L. Cheng, Nanoscale, 2021, 13, 7176–7192.
32. S. Kumar and U. Schwingenschlögl, Chem. Mater., 2015, 27, 1278–1284.
33. J. J. Ma, J. J. Zheng, X. L. Zhu, P. F. Liu, W. D. Li and B. T. Wang, Phys. Chem. Chem. Phys., 2019, 21, 10442–10448.
34. W. Zhang, A. Srivastava, X. Li and L. Zhang, Phys. Rev. B, 2020, 102, 174301.
35. J. He, K. Hummer and C. Franchini, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 075409.
36. B. Amin, N. Singh and U. Schwingenschlögl, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 075439.
37. R. C. Andrew, R. E. Mapasha, A. M. Ukpong and N. Chetty, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 125428.
38. G. Ding, J. He, G. Y. Gao and K. Yao, J. Appl. Phys., 2018, 124, 165101.
39. P. M. Jipin, T. Choudhary and R. K. Biswas, Phys. Scr., 2024, 99, 025301.
40. G. Chen, W. Bao, Z. Wang and D. Tang, Phys. Chem. Chem. Phys., 2023, 25, 9225–9237.
41. I. O. Thomas and G. P. Srivastava, J. Appl. Phys., 2018, 123, 135703.
42. X. Gu, B. Li and R. Yang, J. Appl. Phys., 2016, 119, 085106.
43. N. Tilak, M. Altvater, S. H. Hung, C. J. Won, G. Li, T. Kaleem, S. W. Cheong, C. H. Chung, H. T. Jeng and E. Y. Andrei, Nat. Commun., 2024, 15, 8056.
44. Y. Pan and F. Caruso, npj 2D Mater. Appl., 2024, 8, 42.
45. J. Zhang, X. Li, K. Xiao, B. G. Sumpter, A. W. Ghosh and L. Liang, J. Phys.: Condens. Matter, 2019, 32, 025306.
46. A. Shafique, A. Samad and Y. H. Shin, Phys. Chem. Chem. Phys., 2017, 19, 20677–20683.
47. N. K. Ravichandran and D. Broido, Phys. Rev. X, 2020, 10, 021063.
48. T. Choudhary, J. Peter and R. K. Biswas, ACS Appl. Energy Mater., 2024, 7, 5149–5162.
49. J. Gu, L. Huang and S. Liu, RSC Adv., 2019, 9, 36301–36307.
50. P. F. Liu, T. Bo, J. Xu, W. Yin, J. Zhang, F. Wang, O. Eriksson and B. T. Wang, Phys. Rev. B, 2018, 98, 235426.
51. L. F. Huang and Z. Zeng, J. Appl. Phys., 2013, 113, 083524.
52. R. Farris, O. Hellman, Z. Zanolli, D. Saleta Reig, S. Varghese, P. Ordejón, K. J. Tielrooij and M. J. Verstraete, Phys. Rev. B, 2024, 109, 125422.
53. R. W. Keyes, Phys. Rev., 1959, 115, 564.
54. E. J. Skoug, J. D. Cain and D. T. Morelli, Appl. Phys. Lett., 2010, 96, 181905.
55. G. A. Slack, J. Phys. Chem. Solids, 1973, 34, 321–335.
56. J. Bardeen and W. Shockley, Phys. Rev., 1950, 80, 72–80.
57. L. Cheng and Y. Liu, J. Am. Chem. Soc., 2018, 140, 17895–17900.
58. X. Zhao, G. Tang, Y. Li, M. Zhang and Y. Nie, ACS Appl. Electron. Mater., 2021, 3, 2995–3004.
59. R. K. Biswas and S. K. Pati, ACS Appl. Energy Mater., 2021, 4, 2081–2090.
60. S. P. Keshri, S. K. Pati and A. Medhi, J. Chem. Phys., 2023, 159, 144704.
61. J. Yuan, M. Li and H. Wang, Phys. Rev. B, 2024, 109, 235202.

---

Footnote

† Electronic supplementary information (ESI) available: The data used in and generated by this study. See DOI: https://doi.org/10.1039/d5cp02035b

---