Kagome electronic states in gradient-strained untwisted graphene bilayers

Abstract

Moiré superlattices in twisted homo-bilayers have revealed exotic electronic states, including unconventional superconductivity and correlated insulating phases. However, their fabrication process often introduces moiré disorders, hindering reproducibility and experimental control. Here, we propose an alternative approach using gradient strain to construct moiré superlattices in untwisted bilayer graphene (gs-BLG). Through force-field and first-principles calculations, we show that gs-BLG exhibits kagome-like interlayer-spacing distributions and strain-tunable kagome electronic bands. The competition between interlayer coupling and in-plane strain relaxation leads to distinct structural deformations, giving rise to three forms of diatomic kagome lattices: subtle, pronounced, and distorted. kagome electronic bands are identified near the Fermi level in their band structures. Modulating strain gradients enables tailoring bandwidths and signs of hopping parameters of these kagome bands, providing a versatile platform for studying exotic electronic phases. Our findings establish gradient strain as an alternative to twist engineering, opening an avenue for exploring emergent electronic phases in graphene-based systems.

---

Wei Ji

Wei Ji is a Chair Professor of Physics in the School of Physics, Renmin University of China. He received his PhD degree from the Institute of Physics, CAS in 2008. Prior to joining Renmin University of China, he spent four years at McGill University as a visiting scholar and then a postdoctoral fellow. His research focuses on surface and interface modeling of emerging low-dimensional functional materials. Congratulations on the 10th anniversary of Nanoscale Horizons. We extend our sincere wishes for its continued success and its ongoing impactful contributions to the field of nanotechnology in the years ahead.

---

New concepts

Moiré superlattices, existing in both untwisted and twisted heterostructures as well as twisted homo-structures, have gained increasing attention for their potential to host exotic electronic states. However, the formation and characterization of moiré superlattices in untwisted homo-structures remain largely unexplored, particularly for their geometric and electronic properties. In this work, for untwisted homo-structures, we propose and demonstrate a novel theoretical strategy for constructing moiré superlattices using in-plane strain gradients (gradient strain) in bilayer graphene. By using combined force-field relaxation methods and large-scale density functional theory calculations, we identify three kagome-like patterns in fully relaxed geometries and reveal that the associated kagome electronic bands near the Fermi level are tunable via in-plane strain conditions. These kagome-related properties are rationalized by correlating local stacking-order-dependent interlayer couplings with both in-plane and out-of-plane structural distortions. Furthermore, we establish a connection between hexagonal lattices observed in ridge bilayers and kagome lattices in fully relaxed bilayers via their line graph representations. Our work extends the family of materials hosting novel physical phases to untwisted homo-bilayers, opening exciting avenues for strain-engineering of quantum electronic states in low-dimensional materials.

1. Introduction

Two-dimensional (2D) moiré superlattices can be formed by stacking 2D materials, either homogeneous or heterogeneous, using strategies such as twisting or mismatched lattices. Numerous combinations of lattice-mismatched 2D hetero-bilayers and tunable twist angles diversify the periodicity and symmetry of moiré superlattices.¹ Electronic flat bands and their resulting strong correlation effects were observed in twisted bilayer graphene^2–4 and transition metal chalcogenides.^5,6 A wide spectrum of resulting intriguing states, including unconventional superconductivity,^7–9 correlated insulating states,^10–13 Wigner crystals¹⁴ and rich topological features^15–18 were found in twisted bilayers. Moreover, the slowly varying moiré potential in moiré superlattices significantly promotes the formation of moiré excitons, including both inter- and intra-layer excitons, as well as moiré trions.^19–21

However, even slight structural perturbations can significantly vary the electronic properties of these moiré superlattices.^22–24 For instance, superconducting states observed in magic-angle twisted bilayer graphene may vanish even under fluctuations of 0.1° for the twisting angle or 1 nm for the superlattice periodicity.^25,26 Additionally, the introduction of twisting in most experimental assembly techniques inevitably induces strain fluctuations, referred to as moiré disorders.^22,27–29 The twisting-induced moiré disorders severely limit the reproducibility of experiments and thus the alignment with theoretical findings.^30–34 To address this challenge, researchers employed minimal twist angles (approaching untwisted)³⁵ and large twist angles near 30 degrees,³⁶ combined with strain, to reduce moiré disorder. These strategies facilitated the observation of electronic flat bands and other novel electronic states in the samples. However, exploring exotic electronic states in untwisted moiré superlattices with minimal moiré disorder offers a more straightforward approach. For instance, fractional Chern insulating states were predicted in untwisted but periodically strained monolayer graphene.³⁷ Moreover, a one-dimensional flat band was experimentally observed in an untwisted Bi(110)/SnSe(001) heterostructure.³⁸ Graphene is the strongest 2D material ever measured, which can endure a reversible tensile strain up to 25%,^39–42 showing its potential for forming moiré superlattices implemented by strain. By applying a bubble structure, bilayer graphene exhibits a strain difference of approximately 2% between the two layers.^43–47

Inspired by experimental observations, here, we theoretically construct gradient-strained bilayer graphene (gs-BLG) moiré superlattices by introducing an in-plane strain difference in untwisted bilayer graphene. Using force-field structural relaxations and first-principles density functional theory (DFT) calculations, we reveal kagome-like interlayer-spacing patterns and resulting kagome electronic bands (flat bands and flattened Dirac bands) near the Fermi level. We first present the structural models of gs-BLGs, highlighting differences among the four local stacking configurations. Structural relaxations reveal distinct interlayer-spacing maps forming subtle, pronounced, and distorted diatomic kagome lattices. Our DFT calculations verify the emergence of kagome bands, with wavefunction visualizations validating their real-space distribution. Furthermore, we demonstrate that strain gradients can be used to tune kagome band dispersion, offering an approach to engineering electronic properties. Finally, we establish a direct correlation between out-of-plane corrugation and kagome patterns, demonstrating gradient strain as an effective strategy for constructing moiré superlattices and realizing novel electronic states.

2. Results and discussion

2.1. Geometric structure of gradient-strained bilayer graphene

A gs-BLG moiré superlattice could be constructed by stacking a 3N × 3N and a (3N − 1) × (3N − 1) graphene monolayer supercells together. A series of gs-BLG moiré superlattices can be constructed in varying strain rates by adjusting the supercell size (N) (Fig. S1, ESI†). Fig. 1a shows a gs-BLG model, without structural relaxations, for N = 12, in which the gs-BLG consists of a 35 × 35 (gray, top layer) and a 36 × 36 (red, bottom layer) graphene supercell. Its electronic band structures and real-space distributions of wavefunctions are briefly discussed in Fig. S2 (ESI†). Each gs-BLG moiré superlattice contains four distinct high-symmetry stacking configurations, namely AA, AB, BA and SP, as highlighted in Fig. 1a, and they are shown in detail in Fig. 1b. In the AA stacking (blue in Fig. 1), each atom of the top layer is positioned directly over a carbon atom of the bottom layer, a half-unit-cell lateral shift results in the AB (violet) or BA (orange) stacking. Moreover, there is a transition stacking configuration between AB and BA, referred to as the SP (saddle point) stacking (green). The optimal layer spacings and relative energies of these stacking configurations, together with their atomic models, are depicted in Fig. 1b.

Fig. 1 (a) Top- (upper panel) and side-views (lower panel) of atomic structures of an in-plane gradient strained bilayer graphene moiré superlattice. The top layer (in grey) consists of 3N − 1 × 3N − 1 (here, N = 12) graphene unit cells and the bottom layer has 3N × 3N graphene unit-cells. The top and bottom layer are subjected to potential strains to ensure they have the same in-plane lattice constants for the bilayer supercell. The supercell is marked by the black dashed diamond. Within the diamond, four high-symmetry stacking orders are denoted by differently colored circles, with inset texts showing their labels. (b) Optimal interlayer spacing and relative energy profile of four high-symmetry stacking orders along the red dotted arrow shown in (a). Atomic models are presented for these four stacking orders inside the circles with associated colors.

In the gs-BLG models, the top and bottom layers share the same superlattice constant. The in-plane strains of the two layers can take three primary combinations, determined by the supercell lattice constant: (1) the top layer stretches to conform to the non-strained bottom layer (SN-gs-BLG); (2) the top- (bottom) layer stretches (compresses) to accommodate an averaged lattice constant (SC-gs-BLG); (3) the bottom layer compresses to match the non-strained top layer (NC-gs-BLG). Two additional combinations, where both layers are simultaneously stretched (SS-gs-BLG) or compressed (CC-gs-BLG), also exist and show geometry and electronic structures similar to SN- and NC-gs-BLGs, respectively. Therefore, this work focuses on the first three combinations.

Fig. 2(a–c) show the fully relaxed atomic structures of these three gs-BLGs for N = 12, exhibiting distinct in-plane strain values: 2.86% stretching for the top layer in SN-gs-BLG (case 1, Fig. 2a), ±1.44% averaged strain for SC-gs-BLG (case 2, Fig. 2b), and −2.78% compression for the bottom layer in NC-gs-BLG (case 3, Fig. 2c). Relaxed SN- and SC-gs-BLGs structures feature patterns combining hexagonal (AB and BA regions, violet and orange circles) and triangular (AA regions, blue circles) lattices, as depicted in Fig. 2a and b. However, NC-gs-BLG displays a three-pointed star pattern centered at the AA stacking regions (blue circles). The SP stacking regions (green circles) shifted toward the AA regions, causing the original AB and BA regions to rotate and expand, maximizing their areas (Fig. 2c).

Fig. 2 Atomic structure of gradient-strain bilayer graphene. (a) Top views of the atomic structures for: SN-gs-BLG under 2.86% stretched strain applied to the top layer, (b) SC-gs-BLG under ±1.44% stretched (compressive) strain applied to the top (bottom) layer, and (c) NC-gs-BLG under -2.78% compressive strain applied to the bottom layer. (d)–(f) Projections of the corresponding interlayer spacing of the three, and the white in color bar corresponds to the optimal interlayer spacing of the SP stacking. Black dashed diamonds and differently colored circles represent the supercell and the four stacking orders in the same scheme used in Fig. 1a.

These structural variations are linked to out-of-plane distortions, represented with out-of-plane local atomic displacement (ΔL_Oop), which establish distinct interlayer-spacing (Δd) patterns. Fig. 2(d–f) presents the maps of interlayer spacing for the three gs-BLGs, highlighting the relative out-of-plane distortions. While interlayer spacing is closely tied to stacking order, different in-plane strain combinations significantly reshape spatial distribution patterns. Regions with the smallest interlayer spacing are represented in dark red in Fig. 2(d–f), form subtle diatomic (SN-gs-BLG, Fig. 2d), pronounced diatomic (SC-gs-BLG, Fig. 2e), and distorted diatomic (NC-gs-BLG, Fig. 2f) kagome lattices, delineated by green dashed triangles and lines. Additional visualizations and discussion of these structures and out-of-plane structural corrugations are provided in Fig. S3 (ESI†).

Interlayer-spacing patterns arise from a competition between in-plane elastic (strain) energy and interlayer coupling energy. Both SN- and SC-gs-BLGs exhibit similar interlayer-spacing patterns (Fig. 2d and e), primarily ascribed to biaxial in-plane stretching. While out-of-plane deformation is minimal under in-plane stretched strains, the distribution of in-plane deformation depends on interlayer interactions. In strongly coupled AB and BA stacking domains (red in Fig. 2d and e), local lattice mismatch is minimized to preserve interlayer coupling strength. Conversely, weakly coupled AA stacking domains (in blue) concentrate most in-plane deformation, efficiently relieving stress with minimal interlayer interaction breakdown cost. This redistribution of structural deformation also leads to small out-of-plane corrugations (<0.4 Å) around the AA regions facilitating strain release. Fig. S4 (ESI†) illustrates the in-plane distortion [local atomic displacements (ΔL_IP)] and out-of-plane distortion distributions, supporting this analysis.

Under pure compression in NC-gs-BLG, the in-plane elastic (strain) energy increases substantially, causing the interlayer coupling energy of the original AB and BA domains to be insufficient to balance the heightened strain energy. Thus, out-of-plane distortions become significantly amplified (see Fig. S4, ESI†), more effectively relieving compressive stress. The interlayer spacing in the AA regions increases from 3.4 Å to over 6.5 Å (Fig. 2f). Moreover, strongly coupled AB and BA domains (red in Fig. 2f) expand their areas through in-plane rotations, namely AB domains rotate counterclockwise, while BA domains rotate clockwise, further lowering the total energy of compressed NC-gs-BLG. These rotations drive the SP stacking regions (green circles) to converge, from the midpoints of AB and BA domains, toward the AA regions (Fig. 2c and f). A distorted diatomic kagome lattice forms by connecting regions of small interlayer spacing (dark red regions), as delineated by dark green dashed triangles and lines in Fig. 2f.

2.2. Electronic properties of gradient-strained bilayer graphene

Given the formation of three diatomic kagome lattices, kagome bands are anticipated in the band structures of these gs-BLGs, as confirmed by our DFT calculations. Fig. 3a presents the layer-projected band structure of SN-gs-BLG, highlighting kagome bands (green dashed lines) characterized by a Dirac point (DP1) and a flat band (FB1). A lower-energy nominal FB (residing near the Fermi level at the Gamma point) appears obscured due to perturbations from neighboring bands, which, together with bands FB1 and DP1, exhibit partial characteristics of diatomic kagome bands. The diatomic kagome band features become pronounced in SC-gs-BLG, aligning with the interlayer-spacing map, as illustrated in Fig. 3b. The previously observed DP1 and FB1(in green) persist, while new DP2 and FB2 bands (in brown) emerge, exhibiting quintessential electronic band features of a diatomic kagome lattice.⁴⁸ Fig. 3c presents the band structures of NC-gs-BLG, which reveals two sets of kagome bands, namely KB1/DP1 (in green) and KB2/DP2 (in brown). These six bands represent characteristic bands for a distorted diatomic kagome lattice with opposite signs of hopping parameters. In the absence of spin–orbit coupling (SOC), the states at the Dirac point of KB1 and at the Γ-point of KB2 and DP2 exhibit small gaps due to inversion symmetry breaking. Additionally, orbital projections further reveal that the kagome bands primarily originate from C p_z orbitals, with minor contributions from in-plane p orbitals, as detailed in Fig. S5 (ESI†).

Fig. 3 Electronic properties of gs-BLGs. (a)–(c) Layer atomic projection energy band diagrams for SN-gs-BLG, SC-gs-BLG, and NC-gs-BLG. The colorbar on the right side of (c) indicates the weight distribution between the top and bottom layer in the projection band. Gray-blue colored lines represent bands combinedly contributed from both the top layer and the bottom layer, while the light blue and orange lines indicate dominance by the top layer and bottom layer. Darker blue (red) lines, signify nearly exclusive contribution from the top (bottom) layer. The dashed lines highlight one or two sets of kagome bands. (d)–(g) Plots of |φ|² for the top layer (TL, (d) and (f) and bottom layer (BL, (e) and (g) for the flat band (dark green dot) and the Dirac point (dark gold dot) shown in (a). (h) and (i) and (j) and (k) Plots of |φ|² for the BL (h), (i) and (k) and TL (j) for FB1 (dark green dot), DP2 (bright green dot) and FB2 (brown dot) shown in (b). (l) and (m) and (n) and (o) Plots |φ|² for the TL (l) and (n) and BL (m) and (o) for FB1 (dark green dot) and FB2 (brown dot) shown in (c). The isosurface value of the wavefunctions are 1 × 10⁻⁵ e Å⁻³. Black dashed diamonds and differently colored circles represent the supercell and the four stacking orders in the same scheme used in Fig. 1a.

The real-space wavefunction norm squares |φ|² of these kagome bands further illustrate their intrinsic connection to the kagome lattices. To distinguish the wavefunction distribution properties more clearly, Fig. 3(d–g) present layer resolved |φ|² of the kagome bands in SN-gs-BLG at the K point of the first Brillouin zone (highlighted with solid dots in Fig. 3a). For FB1, the |φ|² in the top layer (Fig. 3d) is predominantly localized in the AB/BA stacking domains, exhibiting hexagonal lattice symmetry. In contrast, the |φ|² in the bottom layer (Fig. 3e) is concentrated in the SP stacking domains (Fig. 3e), forming an electronic kagome lattice, thereby confirming the presence of the kagome flat band in SN-gs-BLG. In contrast, the layer-resolved |φ|² distribution of DP1 (Fig. 3f and g) reveals that the state in the top layer (Fig. 3f) is concentrated within the AA domains, while their interstitial region forms a hexagonal pattern. The state in the bottom layer (Fig. 3g) exhibits a qualitatively complementary pattern, in which the regions with higher charge densities (shown in darker red) form an obscured electronic kagome lattice. This layer-resolved wavefunction behavior emphasizes the critical interplay between stacking domains and the resulting electronic structure of these kagome bands.

For SC-gs-BLG (Fig. 3h–k), the |φ|² in the bottom layer of FB1 (Fig. 3h) exhibits a similar pattern to that of SN-gs-BLG (Fig. 3e), while the distribution of DP1 is analogous to that in SN-gs-BLG and is therefore not shown. For FB2, the |φ|² in the bottom layer (Fig. 3i) shows a delocalized feature mapped with an obscured kagome pattern, which is consistent with its larger dispersion compared to that of FB1. The DP2 states projected on the top (Fig. 3j) and bottom (Fig. 3k) layers exhibit similar features in that the wavefunctions are localized around the AA regions. However, the wavefunction observed at the SP regions is complementary in the two layers.

The wavefunction distributions of FB1 and FB2 in NC-gs-BLG (Fig. 3l–o) are particularly interesting because the distribution of bands on layers is locked. In particular, the projections of FB1 on the top layer (Fig. 3l) and FB2 on the bottom layer (Fig. 3o) share the same distorted kagome lattice feature. In comparison, FB1 on the bottom layer (Fig. 3m) and FB2 on the top layer (Fig. 3n) exhibit another distorted kagome lattice feature. The rotation operations are opposite in direction in these two lattices. Although not strictly identical, these distorted kagome lattices are similar to the geometry and the interlayer-spacing map. These facts indicate that the emergent kagome bands originate from particular distributions of electronic states, driven by structural relaxation-induced local distortions. Moreover, the electronic states and kagome band structures observed in the N = 12 case remain robust within the range of 10 ≤ N ≤ 20. The geometric and electronic structures for the N = 11 case, another analogous example to the N = 12 case, are presented and discussed in Fig. S9 (ESI†) to illustrate the robustness of these kagome features.

Interlayer charge redistributions are also related to the AB and BA stacking domains. Fig. S6 (ESI†) displays the differential charge density (DCD) of these three gs-BLGs, illustrating that interlayer charge redistribution, primarily occurring within AB and BA domains, grows progressively stronger under compressive strain. Furthermore, we plotted the electronic structures of SS-gs-BLG and CC-gs-BLG in Fig. S7 (ESI†), where robust Dirac and flat bands are also observed in their band structures. Given these two additional cases, we reveal a trend for bandwidth variation of kagome bands under gradient strains. Fig. 4a shows monotonic variations in the bandwidths of the Dirac and flat bands as the lattice shrinks. Specifically, we observe that a pure stretched strain favors more linear Dirac and flatter flat bands. As the degree of compressive strain freedom increases, the energy dispersion of Dirac and flat bands converges. This behavior suggests a potential analogy to the angle dependence observed in twisted bilayer graphene.⁴⁹ By applying appropriate gradient strain over larger-sized samples, the Dirac bands near the Fermi level may potentially be fully flattened, paving the way for exotic correlated electronic states.

Fig. 4 (a) Energy dispersion of the flat (black squares) and Dirac (red dots) bands for the five gs-BLGs. The red, yellow, and blue backgrounds represent configurations of pure tension, both tension and compression, and pure compression, respectively. (b) Out-of-plane distortion for top layer of SN-gs-BLG. (c) Out-of-plane distortion for the bottom layer of SC-gs-BLG and (d) the illustration of the kagome line graph framework. (e) Out-of-plane distortion for the bottom layer of NC-gs-BLG and (f) the illustration of the distorted kagome lattice line graph framework, where the cores (grey) are located on the honeycomb sublattice and the linkers (green) on the distorted kagome sublattice. The red and green dashed lines outline the (distorted) kagome lattices. Black dashed diamonds and differently colored circles represent the supercell and the four stacking orders in the same scheme used in Fig. 1a.

2.3. Correlation of electronic states and out-of-plane distortions

Next, we establish the correlation between the kagome electronic states and the out-of-plane structural distortion. Here, the out-of-plane distortions are quantified by projected lengths along the z-axis of C–C bonds (see Fig. S4 for details, ESI†). Positive and negative projection values indicate that the distortion-induced atomic displacements are aligned and the z-axis is anti-aligned. Although the out-of-plane distortion is negligible (<0.01 Å) in SN-gs-BLG, the zero out-of-plane displacement regions (white in Fig. 4b) already form a kagome lattice, as highlighted with red ellipses. Each zero-displacement region spans the two smallest interlayer-spacing regions marked in Fig. 2d, exhibiting its spatial anisotropy in the real-space. Moreover, these regions are also spatially corresponding to the observed charge density of kagome electronic states shown in Fig. 3e and g.

The introduction of compression in SC-gs-BLG leads to larger out-of-plane displacements (>0.10 Å), as shown in Fig. 4c. Unlike the SN-gs-BLG case, the zero-displacement regions are connected into a continuous hexagonal pattern, which covers all smallest interlayer-spacing regions spatially. This hexagonal pattern is analogue to the hexagon lattice shown by gray lines in Fig. 4d. The wavefunctions are visualized in Fig. 3h and i indicate the kagome electronic states are distributed around the middle points of the hexagonal edges (the SP regions). These SP regions represent the nodes of the line graph of the original hexagonal lattice, in other words, a kagome lattice.⁵⁰

The out-of-plane displacements further enlarge to over 0.5 Å in NC-gs-BLG and form a chiral pattern of the out-of-plane zero-displacement regions (white in Fig. 4e), sufficiently relieving in-plane compressive strains. This chiral zero-displacement pattern covers most of the rotated AB and BA domains, consistent with the distorted diatomic kagome lattice observed in the strongly varied interlayer-spacing distribution map. It is also analogous to anti-trichiral structures already used in mechanistic metamaterials,^51–54 as reproduced by gray curved lines shown in Fig. 4f. Like those in SC-gs-BLG, the kagome electronic states are also distributed around the SP regions. Nevertheless, these SP regions converge to the AA regions, not at the middle points of the AB and BA regions, unlike those in SC-gs-BLG. The line graph (Fig. 4f) of the chiral zero-displacement pattern outlines a distorted kagome lattice where the SP regions (green in Fig. 4f) are located at the nodes of the lattice. Moreover, Fig. S8 (ESI†) depicts the maps of local strains for SS- and CC-gs-BLGs.

3. Conclusions

These findings indicate that the interlayer coupling drives out-of-plane structural distortions under in-plane gradient strain, giving rise to three different diatomic kagome patterns. The SN-gs-BLG and SC-gs-BLG primarily release in-plane strain differences through in-plane distortions, sufficient to accommodate the imposed strains. As a result, out-of-plane distortions remain minimal, favoring the interlayer interactions while maintaining the same symmetry as the in-plane distortions, ultimately forming undistorted diatomic kagome lattices. In NC-gs-BLG, however, the in-plane compressive strain is significantly large, making it unlikely to be fully released through in-plane distortions alone. Instead, the out-of-plane distortions dominate, nearly two orders of magnitude larger than the in-plane ones, which compete with interlayer stacking interactions, forming a chiral distorted kagome lattice. The emergence of these strain-induced distortions provides strong evidence for the realization of kagome lattice structures in fully relaxed gs-BLGs.

In summary, we have demonstrated that untwisted bilayer graphene under gradient strains (gs-BLG) can form moiré superlattices that exhibit kagome electronic states. Through detailed force-field calculations for structural relaxations and first-principles calculations for electronic properties, we established that the competition between interlayer coupling and in-plane strain relaxation governs the structural deformations, ultimately leading to the formation of three diatomic kagome lattices under different gradient-strain conditions. The electronic structure of these strained bilayers reveals kagome bands featuring Dirac and flat bands, with their bandwidths tunable via strain modulation. Notably, we observed that compressive strain enhances out-of-plane distortions, which not only significantly influence the electronic states but also induce structural chirality, potentially enabling further tailoring of kagome-related electronic properties. Recent experimental advances have demonstrated that gradient strains in multilayer materials can be introduced up to approximately 3.3% through the use of hetero-substrates for top and bottom layers,⁵⁵ bubble-induced curvature,⁵⁶ and the use of flexible polymers⁵⁷ or high-thermal-expansion metal substrates.⁵⁸ This upper limit exceeds the gradient strain needed for the N = 12 case (∼2.9%). A recent theoretical study has also explored relevant work.⁵⁹ Thus, as a feasible tool for engineering electronic structures in 2D materials, the introduction of gradient strain provides a relatively robust and controllable approach to constructing moiré superlattices, which mitigates the limitations imposed by twist-induced moiré disorders. Our results suggest a promising avenue for the experimental exploration of strain-engineered correlated electronic states in untwisted graphene bilayers and beyond, with implications for novel quantum states.

4. Methods

All structural relaxations of the homogenously strain graphene bilayers were performed using the implementation of a force field, facilitated by the large atomic/molecular massively parallel simulator (LAMMPS).⁶⁰ The C–C interactions within each graphene layer were described using a second-generation reactive empirical bond order (REBO) potential,⁶¹ while the interlayer van der Waals interactions were accounted for via the Kolmogorov–Crespi (KC) potential.⁶² All atoms were fully relaxed until the residual force per atom was less than 1.0 × 10⁻⁴ eV Å⁻¹. The lattice constant of the graphene primitive cell calculated by LAMMPS is 2.46 Å, in agreement with the experiment. The local in-plane strain is defined as , where r represents the position coordinates of a C–C bond in the fully relaxed structure, l is the fully relaxed length of the C–C bond, l⁰_CC = 1.42 Å is the equilibrium bond length of freestanding graphene, and α denotes the angle between the bond and the x–y plane. For simplicity, the in-plane strain is visualized through in-plane geometrical displacement in Fig. S4 (ESI†). Analogously, the local out-of-plane (shear) strain is formulated as . Its corresponding out-of-plane displacement is defined as , providing a direct metric of out-of-plane distortion.

Density functional theory calculations were performed using the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE) for the exchange–correlation potential, linear combined atomic orbitals (LCAO) method, and a single-zeta polarized (SZP) atomic orbital basis set as implemented in the RESCU package.⁶³ The real space mesh resolution was set to 0.35 Bohr, and the convergence criteria for electronic energy and charge density are both set to 10⁻⁵ to ensure good convergence. Two k-meshes of 3  ×  3  ×  1 and 1  ×  1  ×  1 were adopted for calculations on the 11 × 11 + 12 × 12 and 35 × 35 + 36 × 36 supercells, respectively. The models have a vacuum thickness up to about 36 Å in the z-direction to avoid interactions due to periodicity at the surface. For the layer-resolved wavefunction, we obtain the bilayer partial charge density at the K point of the first Brillouin zone, subtract the total charge density of either the top or bottom layer, and visualize the charge accumulation region to represent the layered wavefunction.

Author contributions

Wei Ji and Xianghua Kong conceived and supervised this project. Zeyu Liu, Xianghua Kong, Zewen Wu performed the first principle calculations; Linwei Zhou, Jingsi Qiao, Cong Wang and Shu Ping Lau participated in data analysis and discussion; Zeyu Liu, Xianghua Kong, Wei Ji and Shu Ping Lau wrote the manuscript with inputs from the all authors.

Conflicts of interest

The authors declare no competing financial interest.

Data availability

All data needed to evaluate the conclusions in the paper are present in the paper and/or the ESI.†

Acknowledgements

We thank Drs. Kui Gong, Yibin Hu, and Yin Wang (all from HZWTECH) for helpful discussions. We gratefully acknowledge the financial support from the National Natural Science Foundation of China (Grants No. 52461160327, 92477205, 12474173 and 12104313), the National Key R&D Program of China (Grant No. 2023YFA1406500), the Research Grants Council of Hong Kong (Grant No. CRS_PolyU501/24), the Department of Science and Technology of Guangdong Province (No. 2021QN02L820) and the Shenzhen Science and Technology Program (Grant No. RCYX20231211090126026, the Stable Support Plan Program 20220810161616001), the Fundamental Research Funds for the Central Universities, and the Research Funds of Renmin University of China (Grants No. 22XNKJ30). Calculations were performed at the Physics Lab of High-Performance Computing (PLHPC) and the Public Computing Cloud (PCC) of Renmin University of China.

Notes and references

1. D. M. Kennes, M. Claassen, L. Xian, A. Georges, A. J. Millis, J. Hone, C. R. Dean, D. Basov, A. N. Pasupathy and A. Rubio, Nat. Phys., 2021, 17, 155–163.
2. Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras and P. Jarillo-Herrero, Nature, 2018, 556, 43–50.
3. 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.
4. S. Lisi, X. Lu, T. Benschop, T. A. de Jong, P. Stepanov, J. R. Duran, F. Margot, I. Cucchi, E. Cappelli and A. Hunter, Nat. Phys., 2021, 17, 189–193.
5. M. H. Naik and M. Jain, Phys. Rev. Lett., 2018, 121, 266401.
6. Z. Zhang, Y. Wang, K. Watanabe, T. Taniguchi, K. Ueno, E. Tutuc and B. J. LeRoy, Nat. Phys., 2020, 16, 1093–1096.
7. Y. Saito, J. Ge, K. Watanabe, T. Taniguchi and A. F. Young, Nat. Phys., 2020, 16, 926–930.
8. H. Isobe, N. F. Yuan and L. Fu, Phys. Rev. X, 2018, 8, 041041.
9. M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi and A. Yazdani, Nature, 2021, 600, 240–245.
10. X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi and G. Zhang, Nature, 2019, 574, 653–657.
11. Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, R. C. Ashoori and P. Jarillo-Herrero, Nature, 2018, 556, 80–84.
12. Y. Choi, J. Kemmer, Y. Peng, A. Thomson, H. Arora, R. Polski, Y. Zhang, H. Ren, J. Alicea and G. Refael, Nat. Phys., 2019, 15, 1174–1180.
13. L. Wang, E.-M. Shih, A. Ghiotto, L. Xian, D. A. Rhodes, C. Tan, M. Claassen, D. M. Kennes, Y. Bai and B. Kim, Nat. Mater., 2020, 19, 861–866.
14. E. C. Regan, D. Wang, C. Jin, M. I. Bakti Utama, B. Gao, X. Wei, S. Zhao, W. Zhao, Z. Zhang and K. Yumigeta, Nature, 2020, 579, 359–363.
15. F. Wu, T. Lovorn, E. Tutuc, I. Martin and A. MacDonald, Phys. Rev. Lett., 2019, 122, 086402.
16. K. P. Nuckolls, M. Oh, D. Wong, B. Lian, K. Watanabe, T. Taniguchi, B. A. Bernevig and A. Yazdani, Nature, 2020, 588, 610–615.
17. M. J. Park, Y. Kim, G. Y. Cho and S. Lee, Phys. Rev. Lett., 2019, 123, 216803.
18. Y. Choi, H. Kim, Y. Peng, A. Thomson, C. Lewandowski, R. Polski, Y. Zhang, H. S. Arora, K. Watanabe and T. Taniguchi, Nature, 2021, 589, 536–541.
19. K. L. Seyler, P. Rivera, H. Yu, N. P. Wilson, E. L. Ray, D. G. Mandrus, J. Yan, W. Yao and X. Xu, Nature, 2019, 567, 66–70.
20. C. Jin, E. C. Regan, A. Yan, M. Iqbal Bakti Utama, D. Wang, S. Zhao, Y. Qin, S. Yang, Z. Zheng and S. Shi, Nature, 2019, 567, 76–80.
21. E. Marcellina, X. Liu, Z. Hu, A. Fieramosca, Y. Huang, W. Du, S. Liu, J. Zhao, K. Watanabe and T. Taniguchi, Nano Lett., 2021, 21, 4461–4468.
22. S. Turkel, J. Swann, Z. Zhu, M. Christos, K. Watanabe, T. Taniguchi, S. Sachdev, M. S. Scheurer, E. Kaxiras and C. R. Dean, Science, 2022, 376, 193–199.
23. A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero and E. Zeldov, Nature, 2020, 581, 47–52.
24. C. N. Lau, M. W. Bockrath, K. F. Mak and F. Zhang, Nature, 2022, 602, 41–50.
25. R. Bistritzer and A. H. MacDonald, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 12233–12237.
26. Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu and P. Jarillo-Herrero, Science, 2021, 372, 264–271.
27. A. Uri, S. Grover, Y. Cao, J. A. Crosse, K. Bagani, D. Rodan-Legrain, Y. Myasoedov, K. Watanabe, T. Taniguchi and P. Moon, Nature, 2020, 581, 47–52.
28. S. Grover, M. Bocarsly, A. Uri, P. Stepanov, G. Di Battista, I. Roy, J. Xiao, A. Y. Meltzer, Y. Myasoedov and K. Pareek, Nat. Phys., 2022, 18, 885–892.
29. M. Kapfer, B. S. Jessen, M. E. Eisele, M. Fu, D. R. Danielsen, T. P. Darlington, S. L. Moore, N. R. Finney, A. Marchese and V. Hsieh, Science, 2023, 381, 677–681.
30. P. A. Pantaleón, T. Low and F. Guinea, Phys. Rev. B:Condens. Matter Mater. Phys., 2021, 103, 205403.
31. T. A. de Jong, T. Benschop, X. Chen, E. E. Krasovskii, M. J. de Dood, R. M. Tromp, M. P. Allan and S. J. Van der Molen, Nat. Commun., 2022, 13, 70.
32. N. Nakatsuji and M. Koshino, Phys. Rev. B, 2022, 105, 245408.
33. Y. H. Kwan, G. Wagner, T. Soejima, M. P. Zaletel, S. H. Simon, S. A. Parameswaran and N. Bultinck, Phys. Rev. X, 2021, 11, 041063.
34. C. N. Lau, M. W. Bockrath, K. F. Mak and F. Zhang, Nature, 2022, 602, 41–50.
35. Q. Zheng, C.-Y. Hao, X.-F. Zhou, Y.-X. Zhao, J.-Q. He and L. He, Phys. Rev. Lett., 2022, 129, 076803.
36. S. Li, X. Shi, J. Li, C. He, T. Ouyang, C. Tang and J. Zhong, J. Appl. Phys., 2025, 137, 083902.
37. Q. Gao, J. Dong, P. Ledwith, D. Parker and E. Khalaf, Phys. Rev. Lett., 2023, 131, 096401.
38. Y. Li, Q. Yuan, D. Guo, C. Lou, X. Cui, G. Mei, H. Petek, L. Cao, W. Ji and M. Feng, Adv. Mater., 2023, 2300572.
39. C. Lee, X. Wei, J. W. Kysar and J. Hone, Science, 2008, 321, 385–388.
40. Y. Wei and R. Yang, Natl. Sci. Rev., 2019, 6, 324–348.
41. E. Cadelano, P. L. Palla, S. Giordano and L. Colombo, Phys. Rev. Lett., 2009, 102, 235502.
42. X. Wei, B. Fragneaud, C. A. Marianetti and J. W. Kysar, Phys. Rev. B:Condens. Matter Mater. Phys., 2009, 80, 205407.
43. E. Khestanova, F. Guinea, L. Fumagalli, A. Geim and I. Grigorieva, Nat. Commun., 2016, 7, 12587.
44. A. V. Tyurnina, D. A. Bandurin, E. Khestanova, V. G. Kravets, M. Koperski, F. Guinea, A. N. Grigorenko, A. K. Geim and I. V. Grigorieva, ACS Photonics, 2019, 6, 516–524.
45. J. Zabel, R. R. Nair, A. Ott, T. Georgiou, A. K. Geim, K. S. Novoselov and C. Casiraghi, Nano Lett., 2012, 12, 617–621.
46. Y. Huang, Y.-H. Pan, R. Yang, L.-H. Bao, L. Meng, H.-L. Luo, Y.-Q. Cai, G.-D. Liu, W.-J. Zhao and Z. Zhou, Nat. Commun., 2020, 11, 2453.
47. D. Lloyd, X. Liu, J. W. Christopher, L. Cantley, A. Wadehra, B. L. Kim, B. B. Goldberg, A. K. Swan and J. S. Bunch, Nano Lett., 2016, 16, 5836–5841.
48. S. Zhang, M. Kang, H. Huang, W. Jiang, X. Ni, L. Kang, S. Zhang, H. Xu, Z. Liu and F. Liu, Phys. Rev. B, 2019, 99, 100404.
49. G. Yu, Y. Wang, M. I. Katsnelson and S. Yuan, Phys. Rev. B, 2023, 108, 045138.
50. S. Zhang, M. Kang, H. Huang, W. Jiang, X. Ni, L. Kang, S. Zhang, H. Xu, Z. Liu and F. Liu, Phys. Rev. B, 2019, 99, 100404.
51. R. Hamzehei, S. Rezaei, J. Kadkhodapour, A. P. Anaraki and A. Mahmoudi, Mech. Mater., 2020, 142, 103291.
52. W. Wu, W. Hu, G. Qian, H. Liao, X. Xu and F. Berto, Mater. Des., 2019, 180, 107950.
53. A. Alderson, K. L. Alderson, G. Chirima, N. Ravirala and K. M. Zied, Compos. Sci. Technol., 2010, 70, 1034–1041.
54. A. Alderson, K. L. Alderson, D. Attard, K. E. Evans, R. Gatt, J. N. Grima, W. Miller, N. Ravirala, C. W. Smith and K. Zied, Compos. Sci. Technol., 2010, 70, 1042–1048.
55. F. Yao, D. Rossi, I. A. Gabrovski, V. Multian, N. Hua, K. Watanabe, T. Taniguchi, M. Gibertini, I. Gutiérrez-Lezama, L. Rademaker and A. F. Morpurgo, Nat. Commun., 2024, 15, 10377.
56. P. Ding, J. Yan, J. Wang, X. Han, W. Yang, H. Chen, D. Zhang, M. Huang, J. Zhao, S. Yang, T.-T. Xue, L. Liu, Y. Dai, Y. Hou, S. Zhang, X. Xu, Y. Wang and Y. Huang, Nano Lett., 2024, 24, 8208–8215.
57. Y. Jin, Q. Ren, J. Liu, Y. Zhang, H. Zheng and P. Zhao, Exp. Mech., 2022, 62, 761–767.
58. Y. Que, W. Xiao, X. Fei, H. Chen, L. Huang, S. X. Du and H. J. Gao, Appl. Phys. Lett., 2014, 104, 093110.
59. F. Escudero, A. Sinner, Z. Zhan, P. A. Pantaleón and F. Guinea, Phys. Rev. Res., 2024, 6, 023203.
60. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
61. D. W. Brenner, O. A. Shenderova, J. A. Harrison, S. J. Stuart, B. Ni and S. B. Sinnott, J. Phys.: Condens. Matter, 2002, 14, 783.
62. A. N. Kolmogorov and V. H. Crespi, Phys. Rev. B, 2005, 71, 235415.
63. V. Michaud-Rioux, L. Zhang and H. Guo, J. Comput. Phys., 2016, 307, 593–613.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5nh00307e

‡ These authors contributed equally to this work.

---