Long-Range Machine Learning of Electron Density for Twisted Bilayer Moiré Materials

TL;DR

The paper tackles the computational bottleneck of ab initio electronic structure calculations for moiré 2D materials by extending SALTED, a density-based Gaussian-process framework, with long-range LOVV descriptors to capture interlayer electrostatics. The approach enables predictive extrapolation to large twisted superlattices (thousands of atoms) from training data on small bilayers, supporting downstream properties such as band structures, SOC effects, and domain-wall fields. Key findings show meV-scale accuracy for low-energy bands across graphene, hBN, and TMDCs, with robust performance where locality-based descriptors fail, and substantial speedups over fully converged DFT. The work provides a transferable, first-principles-inspired workflow for exploring moiré physics and designing quantum materials with long-range interlayer interactions.

Abstract

Moiré superlattices in two-dimensional (2D) materials exhibit rich quantum phenomena, but ab initio modelling of these systems remains computationally prohibitive. Existing machine learning methods for accelerating density-functional theory (DFT) can target the prediction of different quantities and often rely on the locality assumption. Here we train a Gaussian process regression SALTED model exclusively on the electron densities of small displaced bilayer structures and then extrapolate electron density prediction to the large supercells required to describe small twist angles between these bilayers. We show the necessity of long-range descriptors to yield reliable band structures and electrostatic properties of large twisted bilayer structures, when these are derived from predicted densities. We demonstrate that the choice of descriptor determines the distribution of residual density errors, which in turn affects the downstream electronic properties. We apply our models to twisted bilayer graphene, hexagonal boron nitride, and transition metal dichalcogenides, focusing on the model's capacity to predict complex phenomena, including flat band formation, bandwidth narrowing, domain-wall electric fields, and spin-orbit coupling effects. Beyond moiré materials, this approach provides a general methodology for electronic structure prediction in large-scale systems with substantial long-range phenomena related to non-local geometric information.

Figures (34)

• Figure 1: Moiré band predictions across diverse twisted bilayer materials. (A) Low-energy moiré band error vs system size for twisted bilayers. SALTED models used here are those with the best band structure prediction performance for each material and descriptor, obtained by following the model optimisation workflow, see \ref{['sec:train-val']} and Supporting Information S2. SOAP fails catastrophically for TiS2 (errors $> 40 \mathrm{meV}$, off-scale) and MoS2 (errors $> 100 \mathrm{meV}$, off-scale), thus not shown here. (B) Comparisons between DFT-converged and SALTED-predicted moiré band structure for graphene, hBN, TiS2, ZrS2, and MoS2, respectively. These five representative twisted bilayer systems have small twist angles $<4^\circ$ and contain $> 1300$ atoms per superlattice. Only the graphene prediction is based on SOAP; the others are based on LOVV. Energy alignment: band gap centred at zero for all materials, with gap central region omitted for hBN, ZrS2, and MoS2. The predictions accurately match DFT band structures across all materials, demonstrating the framework's capability to handle various 2D systems.
• Figure 2: Behaviour of band gaps and band widths with decreasing twist angle. (A) Highlighted flat band(s) in the band structure, with materials' names and twist angles. (B) Band width vs moiré superlattice size (equivalently, twist angle). (C) Band gap vs moiré superlattice size. LODE descriptor predictions for MoS2 are not shown because the prediction error was too large. All structures use commensurate supercells with parallel unit cells in each layer. HVB: highest valence band. LCB: lowest conduction band.
• Figure 3: Prediction of band structure of twisted bilayer TMDCs with spin-orbit coupling. Band structures for (A) TiS2 and (B) ZrS2 at a twist angle $\theta \approx 5.09^\circ$. Grey lines and cyan lines show scalar-relativistic DFT without and with perturbative SOC, respectively. Purple dotted lines show SALTED predictions with perturbative SOC constructed from the predicted density. Bands are aligned to the conduction band region to emphasise the valence band splitting indicated by the red double arrows. The excellent agreement between SALTED+SOC and DFT+SOC demonstrates that SALTED-predicted electron density accurately captures the orbital character necessary for perturbative SOC calculations.
• Figure 4: SALTED predictions for relaxed TBG. (A) SALTED predicted TBG band structure (maroon lines) at $1.61^\circ$, and $1.05^\circ$, compared to full DFT results (blue diamonds). The band structures are predicted by models based on SOAP with GPR regularisation $\eta=10^{-3}$, guided by the band structure metric. (B) The relaxation effect of TBG ($\theta \approx 2.13^\circ$, 2884 atoms/cell). The upper layer deregistration is shown with arrows (in-plane displacement) and colour map (out-of-plane: red for upward, blue for downward, grey for minimal). The energetically unfavourable AA-stacking exhibits counter-clockwise curl with upward displacement (red), while favourable AB-stacking regions show downward displacement (blue). These relaxations minimise energy by reducing AA-stacking area with increased interlayer distance and expanding AB-stacking area with decreased interlayer distance. (C) The upper band gap (between flat bands and the conduction bands above, shown in yellow in (A)) at $\Gamma$ point against twist angles.
• Figure 5: Twisted bilayer hBN relaxation patterns and electric fields. (A) Relaxed TB-hBN structure ($\theta \approx 1.89^\circ$, 3676 atoms/cell) showing atomic position deregistration. The upper layer deregistration is shown with arrows (in-plane displacement) and colour map (out-of-plane: red for upward, blue for downward, grey for minimal). (B) Schematic view of the relaxed TB-hBN structure showing stacking regions (AA, AB, BA) and predicted in-plane electric field 1.5 nm above the bilayer centre. Blue and light red dots represent B and N nuclei respectively, with red arrows indicating the electric field direction and magnitude. The AB-BA domain wall is marked with a yellow line, and the star indicates the in-plane location of the electric field values shown in (c). (c) In-plane electric field intensity perpendicular to the AB-BA domain boundary as a function of twist angle, evaluated at various heights above the bilayer centre. SALTED predictions (SOAP-based model, $\eta=10^{-6}$, optimised for density accuracy) show excellent agreement with full DFT results. The characteristic relaxation length represents reaching equilibrium for the domain wall.