High-throughput bandstructure simulations of van der Waals hetero-bilayers formed by 1T and 2H monolayers

Two-dimensional materials exhibit unique properties from quantum confinement and are promising for future opto/nanoelectronic devices. Beyond individual monolayers, vertically stacked van der Waals heterostructures (vdWHs) enable nearly arbitrary combinations of layers due to weak interlayer bonding, opening a vast design space. However, the combinatorial variety of materials, stackings, and supercell geometries makes conventional high-throughput electronic structure calculations impractical. Standard GGA-DFT underestimates gaps, while higher-level methods are too costly, and lattice mismatch often requires large supercells. This work addresses the research challenge of efficiently predicting full band structures of arbitrary vdWHs for rapid screening. The proposed approach assumes interlayer interactions are dominated by dielectric screening, enabling construction of vdWH band structures from monolayer data with minimal computational cost, and is applied to hexagonal XY₂ (1T/2H) monolayers.

Prior efforts include the Computational 2D Materials Database (C2DB), which compiles monolayer properties and supports band alignment and band gap estimation for prototype heterostructures. Extensive experimental and theoretical studies have focused on specific systems (e.g., graphene, h-BN, MoS₂, MoSe₂, WS₂, WSe₂), revealing phenomena such as tunable interlayer excitons and twist-angle-dependent coupling. Existing high-throughput approaches typically predict lineup and gap but not full heterostructure band dispersions across the Brillouin zone. This work advances beyond prior databases by computing complete heterostructure band structures via band folding of monolayer energy surfaces mapped over the irreducible Brillouin zone.

- Materials selection and database: 93 hexagonal XY₂ monolayers (1T and 2H polytypes) were selected from MaterialsCloud, including dichalcogenides and dihalides, with some metals on the outer layers. For each monolayer (ML), mechanical and electronic properties were computed. - DFT setup: All-electron FHI-AIMS with numeric atom-centered orbitals (tight settings). Structural relaxations with GGA-PBE plus Tkatchenko–Scheffler vdW corrections (Hirshfeld partitioning). k-point sampling scaled with lattice constants. Vacuum ~30 Å along z. For bulk references, kz sampling set by c lattice parameter. - Binding and mechanical properties: Binding energy per area Eb = (Ecrystal − EML)/S using relaxed bulk and ML energies. Strain response obtained by applying biaxial strain in −3% to +3% (0.5% steps), relaxing ions, fitting E(ε) = E0 + αε²; 2D elastic modulus as = α/s. - Electronic structure (monolayers): GGA band structures first; HSE06 with spin–orbit coupling (SOC) applied for nonmetallic, nonmagnetic MLs to obtain accurate band gaps (direct and indirect), band-edge energies (Ev, Ec), and effective mass tensors (second derivatives m−1 = ∂²E/∂ki∂kj computed via three-point stencil near band extrema). SOC reduces gaps by ~0.17 eV on average (max 0.52 eV for PtTe₂). Band-edge strain coefficients kv, kc computed at GGA for efficiency. - Selection for heterostructures: Excluded magnetic MLs and those becoming magnetic at HSE06. Retained 38 nonmetallic, nonmagnetic MLs for vdWH construction. - Heterostructure construction: Coincidence lattice method to find minimal-area common supercells for bilayers with relative twist angles 0–30° (0.1° steps), 2% length tolerance, 0.05° angle tolerance, targeting strains on each ML <2%. Supercell indices {n} define real and reciprocal superlattices; supercell area and BZ scale with N. - Band folding and non-interacting approximation: Compute monolayer electron energy surfaces over a uniform 40×40 grid spanning the 2D irreducible Brillouin zone (IRBZ). Approximate vdWH bands by superimposing ML bands with appropriate BZ folding into the heterostructure supercell. This relies on the assumption that interlayer interaction primarily induces electrostatic screening and small band-edge shifts; Anderson-like lineup holds for type-I/II stacks. - Validation with explicit bilayer DFT: For 45 bilayers (constructed from 10 representative MLs: 2H MoS₂, MoSe₂, WS₂, WSe₂; 1T GeI₂, HfS₂, PtSe₂, SnS₂, SnSe₂, ZrS₂), performed PBE+vdW structural optimization and band calculations, classifying alignment types and comparing band-edge positions to ML values (accounting for strain-induced shifts) to validate the non-interacting model. - High-throughput screening: Generated 703 hetero-bilayers from the 38 MLs (38×37/2), computed band structures along supercell M–K–Γ paths via folding, identified band alignment type (I/II/III), direct/indirect gaps, and band offsets ΔEv, ΔEc. Special analyses for primitive-cell vdWHs (≤2% lattice mismatch) and Γ–Γ systems where both edges lie at Γ.

- Monolayer database: Among 93 MLs, GGA finds 49 with nonzero gaps and 44 metals; 23 show finite spin moments (some gapped, metallic, or half-metallic). After HSE06+SOC and excluding magnetism, 38 MLs are nonmetallic and nonmagnetic. SOC reduces band gaps by ~0.17 eV on average (max 0.52 eV for PtTe₂). Effective masses generally below me; Mo/W dichalcogenides have low, similar electron and hole masses. - Binding and structure trends: Interlayer distance d and thickness h increase with heavier chalcogen/halide; halides thicker on average. Binding energy Eb increases for smaller d and heavier chalcogen/halide. Metal-terminated MLs exhibit very large Eb (>100 meV/Ų) except OTI₂ and STI₂. Elastic moduli as span ~1.4–11.1 eV/Ų (broad spread), dichalcogenides stiffer than halides; 2H stiffer than 1T when both exist. - Validation set of 45 bilayers (explicit PBE+vdW): 8 type-I, 35 type-II, 2 type-III. Band-edge energies in vdWHs lie close to ML values after accounting for strain-induced shifts; remaining deviations attributed to dielectric screening. Non-interacting superposition (Anderson-like rule) holds well for type-I/II; type-III may exhibit charge transfer beyond the model. - High-throughput heterostructures (703 total): 349 type-I, 292 type-II, 62 type-III. Band gaps for type-II vdWHs span ~0 eV (e.g., YbI₂/SnSe₂) to 5.5 eV (MgBr₂/ZnCl₂). Most band offsets within ~1 eV; ΔEv tends to be larger than ΔEc; many cases with ΔEc ≈ ΔEv often correlate with small gaps. Overall gap distribution for type-I and type-II stacks peaks at 1.4–1.6 eV; few structures exceed 3 eV. - Primitive-cell, direct-gap type-II vdWHs: 13 bilayers with minimal unit cell have direct or pseudo-direct gaps, including 2H pairs MoS₂/WS₂, MoSe₂/WSe₂, MoTe₂/WTe₂, and ten 1T halide combinations (e.g., PtTe₂/GeBr₂ gap 0.38 eV with ΔEv 2.72 eV; YbI₂/HgBr₂ gap 1.28 eV with large offsets; high-gap systems such as MgBr₂/CdCl₂). Several 1T-based combinations (e.g., YbI₂/GeI₂, YbI₂/CdI₂) have gaps ~1.6 eV with sizable offsets. - Direct-gap via band folding with small supercells (≤3×3): Multiple 1T/2H and mixed stacks achieve direct gaps due to folding (e.g., YbI₂/HfS₂ Eg 0.17 eV with large offsets; PtO₂/NiO₂ Eg 1.39 eV with offsets >2 eV), offering tunability relevant to sensing. - Γ–Γ class: For MLs with VBM and/or CBM at Γ, Γ–Γ transitions persist irrespective of stacking orientation/size, simplifying fabrication constraints; numerous Γ–Γ type-II combinations identified with a range of gaps and offsets, including wide-gap halide pairs.

The study demonstrates that for type-I and type-II vdWHs, band edges of the heterostructure remain close to those of constituent monolayers once strain-driven shifts are considered, validating an Anderson-like alignment and supporting a non-interacting superposition approach. This enables efficient, accurate estimation of full heterostructure band dispersions via band folding without explicitly simulating large, strained supercells at high computational cost. Mapping 703 bilayers reveals broad tunability of band gaps (near-IR to UV) and offsets critical for device design (e.g., excitonic, photodetector, photovoltaic, and catalytic applications). Identification of primitive-cell direct-gap type-II stacks, including lesser-studied 1T halide bilayers with sizable offsets and mid-gap energies (~1.6 eV), expands candidate materials beyond well-known Mo/W dichalcogenides. The Γ–Γ class underscores cases where device-relevant transitions are robust to stacking variations. Overall, the findings address the challenge of navigating the combinatorial space of vdWHs and provide actionable guidance for experimental synthesis and application-specific selection.

The work introduces and validates a high-throughput, low-cost method to compute full band structures of vdW hetero-bilayers by folding monolayer band energy surfaces, assuming weak interlayer electronic hybridization dominated by dielectric screening. A curated database of 93 XY₂ monolayers provided mechanical and electronic descriptors; 38 nonmagnetic semiconducting MLs yielded 703 bilayers with systematic identification of band alignment types, gaps (≈0–5.5 eV), and offsets. Notably, several primitive-cell, direct-gap type-II heterostructures were uncovered, including ten 1T halide pairs beyond widely studied 2H Mo/W systems, and robust Γ–Γ cases. The approach generalizes to arbitrary 2D symmetries and multilayer stacks, offering a navigation map for the vast vdWH composition space. Future directions include extending to multilayer/multicomponent stacks, incorporating beyond-electrostatic interlayer interactions where hybridization is significant, and benchmarking against many-body (GW/BSE) heterostructure calculations and experiments.

- Core approximation assumes interlayer electronic interaction is limited to classical electrostatic screening; strong interlayer hybridization and charge transfer (notably in type-III alignments) are not fully captured by the non-interacting model. - Magnetic monolayers and cases becoming magnetic at higher-level theory were excluded; results do not cover magnetic vdWHs. - Band-edge strain coefficients were computed at the GGA level (for cost), introducing some inaccuracy relative to HSE06. - GGA underestimates gaps; while corrected with HSE06 for MLs, heterostructure validation calculations were at PBE+vdW level. - vdW functional choices lead to quantitative differences in binding energies across methods; strain limited to ≤2% and minimal-area supercells may not capture all low-energy stacking registries. - Non-uniform treatment of SOC (perturbative) and absence of explicit many-body (GW/BSE) corrections for vdWHs could affect absolute gap values and excitonic effects.