How can we engineer electronic transitions through twisting and stacking in TMDC bilayers and heterostructures? a first-principles approach

Abstract

Layered two-dimensional (2D) materials exhibit unique properties not found in their individual forms, opening new avenues for material exploration. This study examines MX₂ transition metal dichalcogenides (TMDCs), where M is Mo or W, and X is S, Se or Te. These materials are foundational for the creation of hetero- and homo-bilayers with various stacking configurations. Recent interest has focused on twisted homogeneous bilayers, as critical twist angles can significantly alter material properties. This work highlights MX₂ TMDC bilayers with twisted angles that form Moiré patterns, essential to understanding the behaviors of these materials. We performed first-principles calculations using Density Functional Theory (DFT) with range-separated hybrid functionals on 30 combinations of six MX₂ materials with two stacking configurations, revealing that the building blocks and stacking arrangements influence the stability of the heterostructure and the band gap energy (E_g). In particular, the MoTe₂/WSe₂ heterostructure, shifted by 60°, exhibits a direct band gap, indicating potential for novel applications. Our investigation of homobilayers included fully relaxed and low-strain scenarios, examining various stacking styles and twisting angles. Under low-strain conditions, MoS₂, WS₂, and WSe₂ can exhibit direct or indirect band gaps at specific twist angles. Additionally, MoS₂ can transition between semiconductor and conductor states, showcasing diverse electronic properties. Critical twist angles, specifically 17.9° and its corresponding angles (42.1°, 77.9° and 102.1°), in twisted WS₂ and WSe₂ bilayers create symmetric Moiré patterns, leading to direct band gaps. The magnitude of the band gap energy can be tuned by varying the twist angles, which also affect the flatness of the electronic band. Like conventional stacking, most twisted TMDC bilayers exhibit favorable interlayer interactions but with more tailorable characteristics. Using heterostructures and controlled twist angles is a powerful approach in material engineering, enabling the manipulation of various electronic behaviors in advanced materials.

---

1 Introduction

Layered two-dimensional (2D) materials exhibit unique characteristics distinct from their bulk counterparts.¹ The presence of weak van der Waals (vdW) interactions between these layers allows stacking in non-equilibrium configurations, unveiling a plethora of fascinating phenomena.^2,3 This versatility is not limited to heterostructures alone; homogeneous 2D materials are also at the forefront of scientific inquiry because of their ability to manipulate a wide range of properties. The field, abundant with a diverse array of building blocks and customization strategies, still holds many secrets yet to be uncovered. In navigating this vast landscape, first-principles calculations emerge as a crucial tool, providing efficient and insightful guidance in ongoing research within this domain.

Within the vast landscape of van der Waals (vdW) heterostructures, a variety of layered materials such as graphene, boron nitride (BN),⁴ transition metal dichalcogenides (TMDCs),^5,6 and metal carbides and nitrides (MXenes) have been extensively utilized as fundamental building blocks. In this context, even simple bilayer configurations of TMDCs emerge as a robust platform for engineering applications given their shared properties and symmetries. In particular, in certain bilayer TMDC heterostructures, the generation of direct E_g has been observed, paving the way for groundbreaking advancements in optoelectronics.^7,8 The role of symmetry is pivotal in determining the electronic characteristics of these materials. Transitioning TMDCs from their bulk form to 2D structures often involves breaking the inversion symmetry, which can lead to fascinating changes in the structure of the band, including the appearance of spin degeneracy.^9–11 Moreover, the hybridization of the density of states from different layers fosters interlayer charge transfer. This charge transfer plays a crucial role in the creation of an electric dipole moment, significantly contributing to the unique electronic properties of these heterostructures.¹²

Beyond heterostructures, homogeneous layered transition metal dichalcogenides (TMDCs) have been shown to display a wide range of electronic behaviors, heavily influenced by their geometric arrangements.^13,14 For example, bilayer WS₂ has been observed to alternate between direct and indirect E_g depending on its stacking configuration. However, comprehensive studies exploring the full gamut of TMDC material families across various stacking geometries remain limited. Additionally, interlayer coupling in 2D layered materials is sensitive not only to stacking types but also to twisted angles.^15,16

The discovery of exceptional superconductivity in twisted bilayer graphene (TBG) at a specific “magic angle” of 1.1° has sparked significant interest in the study of twisted bilayer structures.^17,18 Since then, this research has expanded, exploring the topological and electronic properties of TBG and similar systems.^19–21 At critical twist angles, symmetric superlattices can form, often characterized by Moiré patterns which create spatial confinement and enhance interlayer coupling. This leads to a flattened band structure and modulates the electronic properties, a phenomenon also observed in double-bilayer graphene.^22,23

Given their hexagonal symmetry and van der Waals (vdW) interlayer interactions, TMDCs are postulated to offer similar opportunities for property manipulation through twisting of homobilayers.^24,25 The optoelectronic behavior of MoS₂, for example, has been closely linked to its twist angles.²⁶ Unlike TBG, where flat bands occur within a narrow range of angles, TMDCs are expected to exhibit flat bands across a wider spectrum of twisted angles.^27–29 This makes TMDCs particularly versatile for engineering band structures. Therefore, twisting the layered configurations of homogeneous TMDC bilayers presents a promising approach to effectively modify their intrinsic properties.

Twisted TMDC bilayers demonstrate distinct electronic behaviors based on moiré angle. Early experiments^30,31 established that smaller moiré patterns exhibit band structure modifications without topological properties. Understanding these simpler systems is crucial as they provide a pathway to engineering topological states through twist angle manipulation. At larger periodicities (θ ≈ 1 – 4°), the systems develop rich topological phenomena, including fractional Chern insulators and quantum spin Hall states.^32–36 This angle-dependent transition from trivial to topological states makes TMDCs a versatile platform for studying and controlling electronic topology. We constrain the system size to ensure accurate treatment of long-range exchange interactions while maintaining computational feasibility to limit the size of the moiré lattices and avoid the accumulation of electron correlation errors inherent in range-separated hybrid functionals like HSE06, particularly due to self-interaction errors and the inability to fully satisfy the generalized Koopmans condition in certain systems.^37,38

In this study, MX₂ transition metal dichalcogenides (TMDCs), where M is Mo or W and X is S, Se, or Te, were selected as the foundational elements for the construction of hetero-bilayers and homo-bilayers, with the aim of uncovering novel phenomena. A significant observation is that most bilayers demonstrate thermodynamically favorable interlayer interaction. This work delves into how the choice of building blocks influences such van der Waals interaction. Initially, we analyzed the combinations of these TMDCs in different stacking configurations, focusing on their E_g and band structures. This was followed by a comparison between two types of stacking in homobilayers and their corresponding monolayers, setting the stage for a discussion on twisted homobilayers. In the context of heterobilayers, it was found that the selection of building blocks and stacking types critically influences the electronic properties. A notable discovery was the direct E_g in the MoTe₂/WSe₂ heterobilayer. For homobilayers, while stacking types present a method to adjust properties, twisting emerges as a more dynamic variable for engineering interlayer couplings. Within specific ranges of twisted angles, homobilayers exhibiting symmetric Moiré patterns can induce regular shifts in E_g, toggle between types of E_g, and affect the flatness of electronic bands. This research paves the way for further exploration of innovative materials, with potential applications in energy conversion, electronics, and optoelectronics, marking a significant advancement in the field of material science.

2 Results and discussion

2.1 Outline

The foundational elements of this research are succinctly summarized in Fig. 1, which sketches how the qualified materials are built, screened, and simulated. In detail, six TMDCs were first selected as the building blocks of the bilayers. A high-throughput coding was implemented to construct numerous preliminary bilayers by continuously rotating angles and applying strain to fit minimum unit cells. Due to structural stability, the need to maintain accurate electron correlation, and the limitation of computational power, we collected the preliminary bilayers with low intrinsic strain, a computable number of atoms, or high symmetry. In particular, employing range-separated hybrid functionals places increased demands on computational resources; thus, selecting moderately sized moiré structures helps ensure both computational tractability and sufficient accuracy of electron correlation. The geometries and electronic properties of the final candidates were computed using a hybrid DFT formulation utilizing range-separated hybrid functionals. The ensuing discussion is divided into different categories of layered materials. These include conventionally stacked hetero- and homo-bilayers, as well as twisted homo-bilayers.

Fig. 1 Workflow diagram illustrating the computational study of transition metal dichalcogenides. Starting from the initial structural information (in this case MoTe₂), the workflow progresses by inputting primitive cells into supercell-forge (a modified version of the supercell-core program). Selected moderate-sized structures undergo symmetry analysis to ensure computational tractability, followed by submission to an ab initio electronic structure workflow utilizing Crystal23 for detailed analysis.

2.2 Hetero-bilayers

In this work, we systematically explored TMDC heterostructures as a result of their remarkable material properties and structural versatility. Using six MX₂ TMDCs (M = Mo, W; X = S, Se, Te) as building blocks, we formed 15 bilayer heterostructures. These were further categorized into two stacking types: 0° and 60° shifts, resulting in 30 unique structures. Each structure was subjected to comprehensive geometric optimization and property analysis. Stability, a key factor, was assessed using eqn (1), which calculates the binding energy by comparing the energy of the bilayers with those of their constituent monolayers. According to Fig. S1,† most heterobilayers exhibit a negative binding energy, indicating a greater energetic favorability than separate monolayers. In particular, WTe₂ bilayers tend to have a positive binding energy, suggesting an unfavorable interlayer interaction. In contrast, bilayers with MoTe₂ show stronger binding, which highlights how the choice of building blocks significantly influences stability. The stacking method slightly alters the energy, as seen in the band gap analysis of these heterostructures in Fig. 2(a). Although the band gaps vary with stacking in all heterostructures, those with 1T-MoS₂ consistently exhibit conductive properties, as Fig. S2† illustrates, underscoring the impact of bonding types on electronic properties. Remarkably, MoTe₂/WSe₂ can toggle between indirect and direct band gaps (1.04 eV). The band gap of 2D materials can transition between indirect and direct as a result of variations in their stacking configurations. This phenomenon arises because of the high sensitivity of the electronic properties of 2D materials to their atomic arrangement and interlayer interactions. Key factors contributing to this behavior include interlayer coupling, symmetry, and strain. The interlayer coupling significantly alters the band structure, with different stacking types modifying the strength and nature of these interactions, thereby affecting the band gap.³⁹ The symmetry of the stacking arrangement influences the overlap of electronic orbitals between layers, where certain configurations enhance or reduce the overlap, leading to a direct or indirect band gap.⁴⁰ Additionally, stacking 2D materials can introduce strain or distortion, which can modify the band structure and potentially convert the band gap between indirect and direct.⁴⁰Fig. 2(b) and (c) also presents the band structures and density of states (DOS) of MoTe₂/WSe₂ with 0° and 60° of rotation, respectively. The detailed orbital projection DOS is presented in Fig. S14.† Both MoTe₂ and WSe₂ contribute to the band edges at 0° of stacking, while MoTe₂ and WSe₂ dominate the edge of the valence and conduction band respectively at 60° of stacking. These findings offer valuable information on the construction of TMDC heterostructures and the tailoring of their properties for electronic applications.

Fig. 2 Comparative analysis of band gap shown for 15 bilayer TMDC heterostructures (a), each evaluated at either a 0 or 60-degree shift. In this plot, the 1H phase of MoS₂ is utilized, whereas Fig. S2† incorporates the 1T phase of MoS₂, highlighting the variations in band gap due to different structural phases of the same material. Band structures and the density of states of MoTe₂/WSe₂ bilayers (b) without rotation, and (c) with rotation. The inserted arrows pinpoint the top of the valence band and the bottom of the conduction band.

2.3 Homo-bilayers

TMDC heterostructures and homogeneous bilayers offer versatile avenues for engineering material properties. This discussion first addresses fully optimized TMDC bilayers, followed by those with slight strain. We considered six MoX₂ TMDCs, each with stacking types A–A and A–B. In A–A stacked bilayers, viewed from above, the cations of each monolayer overlap, whereas in A–B stacking, the cations of one layer overlap with the anions of the other, as depicted in Fig. 3(a)–(d). Compared to their monolayer counterparts, the A–A and A–B stacked bilayers of TMDC exhibit unique electronic properties, as summarized in Fig. 3(e). Key observations include: first, the band gaps of these MX₂ TMDCs range from 0.8 to 2.7 eV. The bonding types significantly influence these properties. For example, both monolayer and bilayer forms of 1T-MoS₂ consistently exhibit conductive behavior, as shown in Fig. S3,† a characteristic attributable to its distinctive structure and interlayer coupling, different from 1H-TMDC. Second, the band gap typically follows the trend: A–B stacking bilayer <A–A stacking bilayer < monolayer. Third, monolayers of 1H-MoS₂, MoSe₂, and WS₂ present direct band gaps, unlike others. This evolution from monolayer to bilayer demonstrates significant electronic changes, hinting at the vast potential inherent in multilayered TMDCs for material engineering.

Fig. 3 Top and side views of transition metal dichalcogenides (TMDCs) bilayers are shown. Panels (a and b) display the A–A and A–B stacking arrangements in the 1T phase of TMDCs, while panels (c and d) depict the A–A and A–B stacking in the 1H phase of TMDCs. Panel (e) shows the comparison of E_g for six fully optimized TMDCs, showcasing their properties in monolayer form as well as in A–A and A–B stacked bilayer configurations. Hollow symbols are used to indicate structures with an indirect E_g, and filled symbols represent those with a direct E_g. It's noteworthy that this plot utilizes 1H-MoS₂, in contrast to Fig. S3,† which presents the corresponding data for the 1T phase of MoS₂, thus providing a comparative understanding of the influence of structural phases on E_g.

In this study, we shifted our focus from conventional A–A and A–B stacking to the more intricate relationships between structure and property in homogeneous TMDC layers, facilitated by varying twisted angles. To refine the scope of our investigation, we initially screened twisted bilayers based on strain magnitude, setting a threshold of less than 0.001 post-formation. The comprehensive trends of strain versus twisted angles for all six TMDCs are depicted in Fig. S4.† These trends are further elaborated in Fig. S5,† which zooms in on regions showcasing geometries that meet the strain criteria. This screening approach is based on the hypothesis of an inverse relationship between structural stability and strain. Subsequently, the low-strain structures underwent symmetrization, and those with more than 40 atoms per unit cell were excluded to optimize computational efficiency. The final candidates selected, adhering to these criteria, are highlighted in Fig. S5,† marking a strategic narrowing of the potential structures for detailed analysis. This process not only aids in the identification of structurally viable configurations but also sets the stage for a deeper exploration of the nuanced interplay between twisted angles and material properties in TMDC bilayers.

Simultaneously, Moiré patterns emerged prominently among these selected candidates. Subsequently, these structures underwent symmetrical optimization while maintaining fixed unit cell sizes, preserving both low strain and the captivating Moiré patterns capable of eliciting novel material properties. For example, consider bilayer WSe₂: Fig. 4 vividly illustrates the evolution of Moiré patterns at twisted angles from 0° to 60°. Consequently, the band structures exhibit notable variations, with pronounced flattening due to Moiré patterns' influence on interlayer coupling interactions. This enhanced flatness facilitates energy conversion by mitigating the impact of crystal momentum. Furthermore, the interlayer interaction of these twisted bilayers was meticulously examined, with Fig. 5 showing the binding energy as a function of twisted angles.

Fig. 4 Geometries and band structures of bilayer WSe₂ at various twisted angles are shown for (a) 0°, (b) 17.9°, (c) 27.8°as well as MoTe₂ (d–f) for the corresponding angles. Each panel provides both a top view and a side view of the bilayer geometries (shown in the bottom left of each section). For a broader perspective, analogous analyses for other materials including 1H-MoS₂, 1T-MoS₂, MoSe₂, WS₂, and WTe₂ are comprehensively presented in Fig. S7–S12.† This extensive representation allows for a detailed comparison across a range of materials and twisted angles, offering a deeper understanding of the relationship between geometric configurations and electronic properties in these layered structures.

Fig. 5 Variation in binding energy as a function of twisted angles for six twisted bilayer cases: (a) MoS₂, (b) MoSe₂, (c) MoTe₂, (d) WS₂, (e) WSe₂, and (f) WTe₂. The red and gray dashed lines denote the binding energies for A–A and A–B stacked homobilayers, respectively, providing a reference for comparison against the fluctuating binding energies observed in the twisted bilayers. The vertical green bars represent the approximate range for the “Magic Angles” ∼1–4°.

In essence, the binding energy of the twisted bilayers exhibits a strong correlation with their geometry, often oscillating between the values observed for A–A and A–B stacked bilayers. Because the A–B-stacked bilayers have a more negative binding energy compared to the A–A-stacked ones, the A–B-stacked interaction is more thermodynamically favorable. A more negative binding energy implies that the system is in a relatively low-energy state, requiring more energy to separate the layers. This phenomenon also occurs in the context of bilayer graphene due to a more favorable interlayer interaction in stacking of A–B⁴¹ In particular, except for the WTe₂ bilayer with a 0° twist, the majority of twisted homobilayers display a negative binding energy within the values of the analogies of stack A–A and A–B. This indicates that most bilayer twistronics in our work can possess a relatively favorable interlayer interaction even compared to conventional stackings. Intriguingly, 1T-MoS₂ is predicted to experience a substantial reduction in binding energy at specific twist angles. These findings underscore the role of relatively favorable van der Waals interactions in enabling nonequilibrium stacking configurations, providing insights into the construction of layered twistronics. Furthermore, Fig. 6 visually shows how the band gaps of various TMDC homobilayers varies with twisted angles, revealing different behaviors among different bilayers.

Fig. 6 Variation in band gap across a range of twisted angles for six twisted bilayer materials: (a) MoS₂, (b) MoSe₂, (c) MoTe₂, (d) WS₂, (e) WSe₂, and (f) WTe₂. The graph captures the changes in band gap as a function of the twisting angle in these bilayer structures. The vertical green bars again represent the approximate range for the “Magic Angles”.

For example, the behavior of 1T-MoS₂ is particularly intriguing, as it exhibits a transition between a conductor and a semiconductor depending on the twisted angle, as shown in Fig. 6(a). Interestingly, at a twist angle of 92.2°, 1T-MoS₂ even achieves a direct band gap with a value of 0.232 eV. In contrast, MoSe₂ appears to be less responsive to changes in the band gap-twist relationship, maintaining an indirect bandgap within a limited range of variation (1.79–1.81 eV) at different twist angles, as shown in Fig. 6(b). In contrast, the other four TMDCs (MoTe₂, WS₂, WSe₂, and WTe₂) exhibit symmetric band gap-twist relationships, as illustrated in Fig. 6(c)–(f), with respect to 60°twist angles. These TMDCs exhibit maximum band gap values at 0° twist and minimum band gap values at 60° twist. In particular, both WS₂ and WSe₂ undergo a transition from indirect to direct band gap at specific twist angles, namely 17.9° and its corresponding angles (42.1°, 77.9° and 102.1°). These findings highlight that the band gap values of these four twisted TMDCs fluctuate between the band gap values observed in their stacked A–A and A–B bilayers, as demonstrated in the fully optimized bilayer section.

In general, the tunability of these band gaps is evident across a range of twist angles, as: 1H-MoS₂: 1.87–2.16 eV; 1T-MoS₂: 0.00–0.80 eV; MoSe₂: 1.79–1.81 eV; MoTe₂: 0.81–1.49 eV; WS₂: 1.62–2.28 eV; WSe₂: 1.55–2.10 eV; WTe₂: 0.92–1.63 eV. Specifically, 0.2 to 0.3 eV of differences in band gap values can be achieved between a twisting step. Even a difference of 0.1 eV in material band gaps can be significant in various properties and applications.⁴² Not only does a tiny change in band gap affect the conductivity but also influences the optical features. For example, a difference of 0.1 eV can enhance the efficiency of sensors and can change the color of light-emitting diodes (LEDs).^43,44 Apart from the band gap value, the change in band gap types also has important implications, even if the band gap values are the same. Materials with direct band gaps can efficiently emit and absorb light, whereas materials with indirect band gaps are poor, since the electron transition requires a phonon to conserve momentum.⁴⁵ On the other hand, a direct band gap typically has a shorter carrier lifetime because the recombination of electrons and holes can occur more readily. However, the indirect band gap has a longer carrier lifetime, which can be advantageous for applications such as transistors.⁴⁵ In short, twisting bilayers offer a sophisticated approach to fine-tuning electronic behavior that can be valuable for different applications. The presence of highly symmetric Moiré patterns is believed to induce diverse electronic interactions in low-strain homogeneous bilayer materials.

3 Methods

The methodology of our study is divided into two segments: the establishment of prototypes and ab initio calculations.

3.1 Establishment of TMDC bilayers

The heterostructures were generated using SUPERCELL-FORGE, an extension of the SUPERCELL-CORE program,⁴⁶ which allows the rapid generation and evaluation of heterostructures at various twist angles. Although the supercell core handles the matching of the in-plane lattice using the Zur and McGill method,⁴⁷SUPERCELL-FORGE incorporates additional functionality for the optimization of the out-of-plane geometry. The interlayer separation is determined through experimental data-based modeling with the option of applying a machine learning force field for geometry optimization. FINDSYM⁴⁸ is integrated to analyze the underlying symmetries of the generated structures, identifying space groups, lattice parameters and atomic positions within specified tolerances. The program evaluates candidate structures based on multiple criteria, including lattice strain, unit cell size, and total atom count, enabling efficient screening of potential heterostructures. For prototypes of fully relaxed heterostructures and homogeneous bilayers, we confine the rotational angles to either 0° or 60°, in alignment with the hexagonal symmetry characteristic of the TMDCs studied. We also considered multiple unit cell sizes for each angle, ensuring that the total number of atoms remained below 30 to accommodate computational limitations.

The program SUPERCELL-FORGE was applied to generate the initial models of low-strain yet symmetrized homogeneous bilayers, the candidate structures were selected taking into account structural stress less than 0.001 (Fig. S4 and S5†) over a range of angles with increments of 0.1° from 0° to 120°. These structures were filtered further via the program to fall within certain integer multiples of the lattice parameters, 5, and were restricted to 40 atoms per unit cell, balancing the need for detailed structural representation with computational feasibility. The FINDSYM tolerance parameters within the program were set to a lattice tolerance of 0.001 and 0.01 for the atomic positions. The interlayer distance was determined by a model that included known multilayered TMDC materials. These restrictions allow for materials which are computationally tractable without the use of excessive computational resources and limit the need for considering long-range electron correlation which is outside of the scope of our current DFT-based methods. Finally, the resulting symmetrized materials were allowed to fully relax through the SCF process within Crystal23.

After such a series of selections, scenarios like twistronics with distorted phases of TMDCs or with small rotational angles are excluded from the targets of this work. These can lead to considerable supercells with low symmetries incalculable by the DFT at the level of hybrid functionals. Examples of such supercells containing around 500 atoms per cell are illustrated in Fig. S13.† In contrast to studies on bilayer graphene at the magic angle, which employ a combination of DFT at more approximate levels of functionals alongside other approaches, such as tight-binding models or molecular dynamics,^49,50 the DFT approach with hybrid functionals requires massive computational resources for systems containing hundreds of atoms due to the algorithmic scaling of DFT-like methods with the number of orbitals/atoms in the material's asymmetric unit.

3.2 Computational settings

In the second phase of our study, the previously established prototypes were subjected to geometric optimization and property analysis using unrestricted and periodic DFT. We used the hybrid functional HSE06,^51,52 as implemented in the Crystal23 code,⁵³ and integrated it with D3 dispersion correction⁵⁴ and Becke–Jonesson damping.⁵⁵ For Gaussian basis sets, the triple-zeta valence with polarization (TZVP)⁵⁶ was selected. During geometric optimization, we performed full optimization for both heterostructures and homogeneous bilayers. For the unconventional homo-bilayers with exotic rotational angles, the geometric optimization was performed with fixed lattice constants. The optimization process adhered to stringent convergence criteria, with an RMS force threshold of 1.54 × 10⁻² eV Å⁻¹, a maximum force of 2.31 × 10⁻² eV Å⁻¹, an RMS displacement of 6.35 × 10⁻⁴ Å, and a maximum displacement of 9.53 × 10⁻⁴ Å. To extract detailed electronic information in the property calculations, we projected the first Brillouin zone (∼2π × 1/60 Å⁻¹) onto a dense Pack–Monkhorst k-mesh grid (k_a × k_b × k_c). The mesh grid used much higher shrinking factors (with a·k_a ≥ 60, b·k_b ≥ 60, c·k_c ≥ 60) compared to those used in the geometric optimization phase (where a·k_a ≥ 40, b·k_b ≥ 40, c·k_c ≥ 40).⁵⁷ Also, band structures were computed with the total number of k points along the path set to 1000. These DFT approaches ensured the acquisition of high-resolution data that is essential for a comprehensive understanding of the electronic properties of these complex systems. Apart from the DFT algorithm settings, we executed the computation with 32 processors (Intel Xeon Gold 6148 Processor) on one node. The CPU time requires a geometric optimization ranging from 2 to 16 million seconds.

The calculation of the binding energy per atom (E_BE) for the bilayers in our study is performed using the following equation:

In this equation, E_BL denotes the total energy per unit cell of the bilayer, while E_ML1 and E_ML2 represent the total energies per unit cell of the first and second monolayer building blocks, respectively. The terms N_BL, N_ML1, and N_ML2 correspond to the number of atoms in the unit cell of the bilayer and the first and second monolayers, respectively.

Conclusions

In summary, this work demonstrates how manipulating the stacking and twist angles in bilayer TMDCs (MX₂, with M = Mo, W and X = S, Se, Te) unlocks a broad spectrum of electronic and thermodynamic behaviors. By combining systematic screening of low-strain supercells and range-separated hybrid DFT calculations, we map out how weak van der Waals interactions can be exploited to fine-tune the band gap magnitudes, conductivity, and direct–indirect gap transitions across both heterobilayers and homobilayers. In particular, the MoTe₂/WSe₂ heterobilayer undergoes a striking transformation from indirect to direct band gaps, underscoring its potential in energy-conversion devices. Meanwhile, 1T-MoS₂ exhibits a unique conductor–semiconductor switching at certain twist angles, revealing intriguing prospects for electronic switching applications.

The emergence of Moiré patterns at specific twist angles plays a central role in flattening the band structure, thus enhancing carrier interactions and, in some cases, facilitating the design of more efficient optical and electronic devices. Early experiments^30,31 in twisted TMDC systems provided key insights into how smaller Moiré patterns can induce band structure modifications without immediate topological features. More recent developments indicate that larger periodicities can lead to the formation of fractional Chern insulators and correlated states in twisted bilayer MoTe₂.^58,59 By systematically varying angles up to 60°, we identify critical twists—e.g., around 17.9° in WS₂ and WSe₂ —that induce direct band gaps, expanding the design possibilities for light emission and sensing technologies. Moreover, our binding energy analysis confirms that MoTe₂-based bilayers exhibit robust interlayer coupling, while WTe₂-based bilayers are relatively less stable, demonstrating the intricate interplay between material composition and thermodynamic favorability.

In the future, to fully realize the potential of twisted TMDC bilayers, several critical directions must be pursued. Spin–orbit coupling effects must be explicitly included, as they significantly influence the band topology and electronic properties of these systems. Additionally, the limitations of range-separated hybrid functionals like HSE06, particularly the self-interaction and Koopmans' condition errors,³⁷ necessitate the exploration of higher-level functionals and beyond-DFT methods, including exact diagonalization,^60,61 embedded dynamical mean-field theory (eDMFT),⁶² and quantum Monte Carlo approaches.^63,64 Developing a robust band unfolding procedure for Crystal23 is crucial to examine large Moiré lattices, allowing a detailed analysis of electronic properties. Furthermore, rigorous benchmarking of DFT-based methods for twisted TMDCs with varying Moiré lattice parameters will provide a reliable framework for evaluating and predicting their properties.

By illuminating how stacking configurations, twist angles, and material choices converge to shape electronic properties, this study sets the stage for designing next-generation 2D materials. Whether for energy-harvesting, optoelectronics, or quantum applications, twisted TMDC bilayers represent a versatile and rapidly evolving platform poised to reveal new phases of matter and to enable a range of transformative device technologies.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Author contributions

Conceptualization: YL, WPC, JLMC. Methodology: YL, WPC, JLMC. Software: YL, WPC, JLMC. Validation: YL, JLMC. Formal analysis: YL, JLMC. Investigation: YL, JLMC. Resources: JLMC. Writing—original draft preparation: YL. Writing—review and editing: YL, WPC, JLMC. Visualization: YL, WPC. Supervision: JLMC. Project administration: JLMC. Funding Acquisition: JLMC. All authors have read and agreed to the published version of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported in part by computational resources and services provided by the Institute for Cyber-Enabled Research at Michigan State University.

References

1. D. Voiry, A. Mohite and M. Chhowalla, Chem. Soc. Rev., 2015, 44, 2702–2712.
2. C. Long, Y. Dai, Z.-R. Gong and H. Jin, Phys. Rev. B, 2019, 99, 115316.
3. Z. Shi, X. Wang, Y. Sun, Y. Li and L. Zhang, Semicond. Sci. Technol., 2018, 33, 093001.
4. C. R. Dean, A. F. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim and K. L. Shepard, et al. , Nat. Nanotechnol., 2010, 5, 722–726.
5. P. K. Sahoo, S. Memaran, Y. Xin, L. Balicas and H. R. Gutiérrez, Nature, 2018, 553, 63–67.
6. C. Jin, E. C. Regan, A. Yan, M. Iqbal Bakti Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng and S. Shi, et al. , Nature, 2019, 567, 76–80.
7. H. Terrones and M. Terrones, J. Mater. Res., 2014, 29, 373–382.
8. K. Wang, B. Huang, M. Tian, F. Ceballos, M.-W. Lin, M. Mahjouri-Samani, A. Boulesbaa, A. A. Puretzky, C. M. Rouleau and M. Yoon, et al. , ACS Nano, 2016, 10, 6612–6622.
9. A. M. Jones, H. Yu, J. S. Ross, P. Klement, N. J. Ghimire, J. Yan, D. G. Mandrus, W. Yao and X. Xu, Nat. Phys., 2014, 10, 130–134.
10. J. Sung, Y. Zhou, G. Scuri, V. Zólyomi, T. I. Andersen, H. Yoo, D. S. Wild, A. Y. Joe, R. J. Gelly and H. Heo, et al. , Nat. Nanotechnol., 2020, 15, 750–754.
11. D. Xiao, G.-B. Liu, W. Feng, X. Xu and W. Yao, Phys. Rev. Lett., 2012, 108, 196802.
12. L. Li and M. Wu, ACS Nano, 2017, 11, 6382–6388.
13. X. Zhang, H. Nan, S. Xiao, X. Wan, X. Gu, A. Du, Z. Ni and K. Ostrikov, Nat. Commun., 2019, 10, 598.
14. S. Kim, A. Konar, W.-S. Hwang, J. H. Lee, J. Lee, J. Yang, C. Jung, H. Kim, J.-B. Yoo and J.-Y. Choi, et al. , Nat. Commun., 2012, 3, 1011.
15. B. Tang, B. Che, M. Xu, Z. P. Ang, J. Di, H.-J. Gao, H. Yang, J. Zhou and Z. Liu, Small Struct., 2021, 2, 2000153.
16. Y.-C. Lin, R. Torsi, R. Younas, C. L. Hinkle, A. F. Rigosi, H. M. Hill, K. Zhang, S. Huang, C. E. Shuck and C. Chen, et al. , ACS Nano, 2023, 17, 9694–9747.
17. Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras and P. Jarillo-Herrero, Nature, 2018, 556, 43–50.
18. Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi and E. Kaxiras, et al. , Nature, 2018, 556, 80–84.
19. F. Wu, T. Lovorn, E. Tutuc, I. Martin and A. MacDonald, Phys. Rev. Lett., 2019, 122, 086402.
20. M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young and C. R. Dean, Science, 2019, 363, 1059–1064.
21. K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald and E. Tutuc, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 3364–3369.
22. R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 12233–12237.
23. Y. Cao, D. Rodan-Legrain, O. Rubies-Bigorda, J. M. Park, K. Watanabe, T. Taniguchi and P. Jarillo-Herrero, Nature, 2020, 583, 215–220.
24. K. Liu, L. Zhang, T. Cao, C. Jin, D. Qiu, Q. Zhou, A. Zettl, P. Yang, S. G. Louie and F. Wang, Nat. Commun., 2014, 5, 4966.
25. A. A. Puretzky, L. Liang, X. Li, K. Xiao, B. G. Sumpter, V. Meunier and D. B. Geohegan, ACS Nano, 2016, 10, 2736–2744.
26. S. Huang, X. Ling, L. Liang, J. Kong, H. Terrones, V. Meunier and M. S. Dresselhaus, Nano Lett., 2014, 14, 5500–5508.
27. L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai and B. Kim, et al. , Nat. Mater., 2020, 19, 861–866.
28. E. S. Morell, J. Correa, P. Vargas, M. Pacheco and Z. Barticevic, Phys. Rev. B:Condens. Matter Mater. Phys., 2010, 82, 121407.
29. Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc and B. J. LeRoy, Nat. Phys., 2020, 16, 1093–1096.
30. E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang, K. Yumigeta, M. Blei, J. D. Carlström, K. Watanabe, T. Taniguchi, S. Tongay, M. Crommie, A. Zettl and F. Wang, Nature, 2020, 579, 359–363.
31. Y. Tang, L. Li, T. Li, Y. Xu, S. Liu, K. Barmak, K. Watanabe, T. Taniguchi, A. H. MacDonald, J. Shan and K. F. Mak, Nature, 2020, 579, 353–358.
32. H. Park, J. Cai, E. Anderson, Y. Zhang, J. Zhu, X. Liu, C. Wang, W. Holtzmann, C. Hu, Z. Liu, T. Taniguchi, K. Watanabe, J.-H. Chu, T. Cao, L. Fu, W. Yao, C.-Z. Chang, D. Cobden, D. Xiao and X. Xu, Nature, 2023, 622, 74–79.
33. K. Kang, B. Shen, Y. Qiu, Y. Zeng, Z. Xia, K. Watanabe, T. Taniguchi, J. Shan and K. F. Mak, Nature, 2024, 628, 522–526.
34. Z. Ji, H. Park, M. E. Barber, C. Hu, K. Watanabe, T. Taniguchi, J.-H. Chu, X. Xu and Z.-X. Shen, Nature, 2024, 635, 578–583.
35. K. Kang, B. Shen, Y. Qiu, Y. Zeng, Z. Xia, K. Watanabe, T. Taniguchi, J. Shan and K. F. Mak, Nature, 2024, 628, 522–526.
36. H. Park, J. Cai, E. Anderson, Y. Zhang, J. Zhu, X. Liu, C. Wang, W. Holtzmann, C. Hu, Z. Liu, T. Taniguchi, K. Watanabe, J.-H. Chu, T. Cao, L. Fu, W. Yao, C.-Z. Chang, D. Cobden, D. Xiao and X. Xu, Nature, 2023, 622, 74–79.
37. S. T. Holmes, R. J. Iuliucci, K. T. Mueller and C. Dybowski, J. Chem. Theory Comput., 2015, 11, 5229–5241.
38. V. Ivády, I. Abrikosov, E. Janzén and A. Gali, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 87, 205201.
39. Y. Oh, S. Song and J. Bae, Int. J. Mol. Sci., 2024, 25, 13104.
40. N. V. Vinh, D. Lu and K. Pham, Nanoscale Adv., 2025, 790–799.
41. S. Zhou, J. Han, S. Dai, J. Sun and D. J. Srolovitz, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 92, 155438.
42. M. Grundmann, Physics of Semiconductors, Springer, 2010, vol. 11, pp. 401–472.
43. S. M. Sze, Y. Li and K. K. Ng, Physics of Semiconductor Devices, John wiley & sons, 2021, pp. 697–714.
44. J. Singh, Semiconductor Devices: Basic Principles, John Wiley & Sons, 2000, pp. 460–502.
45. L.-D. Yuan, H.-X. Deng, S.-S. Li, S.-H. Wei and J.-W. Luo, Phys. Rev. B, 2018, 98, 245203.
46. T. Necio and M. Birowska, AIP Adv., 2020, 10, 105105.
47. A. Zur and T. C. McGill, J. Appl. Phys., 1984, 55, 378–386.
48. H. T. Stokes and D. M. Hatch, J. Appl. Crystallogr., 2005, 38, 237–238.
49. S. Carr, D. Massatt, S. Fang, P. Cazeaux, M. Luskin and E. Kaxiras, Phys. Rev. B, 2017, 95, 075420.
50. N. Leconte, S. Javvaji, J. An, A. Samudrala and J. Jung, Phys. Rev. B, 2022, 106, 115410.
51. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
52. A. V. Krukau, O. A. Vydrov, A. F. Izmaylov and G. E. Scuseria, J. Chem. Phys., 2006, 125, 224106.
53. R. Dovesi, A. Erba, R. Orlando, C. M. Zicovich-Wilson, B. Civalleri, L. Maschio, M. Rérat, S. Casassa, J. Baima and S. Salustro, et al. , Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2018, 8, e1360.
54. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
55. A. D. Becke and E. R. Johnson, J. Chem. Phys., 2005, 123, 154101.
56. D. Vilela Oliveira, J. Laun, M. F. Peintinger and T. Bredow, J. Comput. Chem., 2019, 40, 2364–2376.
57. H. J. Monkhorst and J. D. Pack, Phys. Rev. B:Condens. Matter Mater. Phys., 1976, 13, 5188.
58. J. Yu, J. Herzog-Arbeitman, M. Wang, O. Vafek, B. A. Bernevig and N. Regnault, Phys. Rev. B, 2024, 109, 045147.
59. Y. Jia, J. Yu, J. Liu, J. Herzog-Arbeitman, Z. Qi, H. Pi, N. Regnault, H. Weng, B. A. Bernevig and Q. Wu, Phys. Rev. B, 2024, 109, 205121.
60. D. Medvedeva, S. Iskakov, F. Krien, V. V. Mazurenko and A. I. Lichtenstein, Phys. Rev. B, 2017, 96, 235149.
61. Y. Lu and M. W. Haverkort, Eur. Phys. J.:Spec. Top., 2017, 226, 2549–2564.
62. G. Kotliar and D. Vollhardt, Phys. Today, 2004, 57, 53–59.
63. K. Haule, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 75, 155113.
64. B. Bauer, L. D. Carr, H. G. Evertz, A. Feiguin, J. Freire, S. Fuchs, L. Gamper, J. Gukelberger, E. Gull, S. Guertler, A. Hehn, R. Igarashi, S. V. Isakov, D. Koop, P. N. Ma, P. Mates, H. Matsuo, O. Parcollet, G. Pawłowski, J. D. Picon, L. Pollet, E. Santos, V. W. Scarola, U. Schollwöck, C. Silva, B. Surer, S. Todo, S. Trebst, M. Troyer, M. L. Wall, P. Werner and S. Wessel, J. Stat. Mech., 2011, 2011, P05001.

---

Footnote

† Electronic supplementary information (ESI) available: Additional comprehensive plots for binding energy, electronic properties, and stain analysis; and structural CIF files for optimized 2D-bilayers at different twisted angles. See DOI: https://doi.org/10.1039/d5na00112a

---