Engineering Macroscopic Quantum Coherence in Disordered Moiré Metamaterials
Engineering
Macroscopic Quantum Coherence in Disordered Moiré Metamaterials
Author: Rowan Brad Quni-Gudzinas
Contact: rowan.quni@outlook.com ORCID:
ISNI: 0000000526456062
DOI: 10.5281/zenodo.18014173 Date:
2025-12-24 Version: 1.0.1
Abstract: The study of moiré materials is defined by
a central paradox: the experimental observation of robust, macroscopic
quantum phases, such as superconductivity and topological insulation, in
devices known to possess significant atomistic disorder. This work
bridges this “disorder-topology gap” by establishing quantum percolation
as the fundamental mechanism for the emergence of global coherence. We
develop a Stochastic Time-Dependent Ginzburg-Landau (TDGL) model that
incorporates a physically realistic, multi-modal disorder potential.
This potential includes both a smooth, spatially correlated twist-angle
field, representing mesoscopic strain domains, and a sparse field of
sharp pinning sites, representing atomic vacancies. The simulation,
performed on a high-resolution 128x128 lattice, is analyzed using an
objective Hoshen-Kopelman cluster-finding algorithm. Our results
quantitatively map the entire transition, from the initial nucleation of
isolated superconducting “islands” in magic-angle regions to their
proximity-driven growth and coalescence. We identify a sharp,
topological percolation threshold at a simulation time of t=1.25, where
a sample-spanning network first forms at an average superfluid density
of 0.31. This globally coherent state stabilizes despite a significant
twist-angle disorder bandwidth of 0.15°. This work reframes disorder
from a nuisance to be eliminated into a key engineering parameter,
providing a predictive framework for designing the next generation of
robust, scalable moiré quantum devices.
Keywords: Moiré Materials, Quantum Percolation,
Ginzburg-Landau Theory, Straintronics, Disordered Superconductivity,
Twistronics, Computational Physics
1.0 INTRODUCTION & PROBLEM
STATEMENT
1.1 The Moiré
Metamaterial: From Observation to Design
The field of moiré physics has undergone a profound maturation,
transitioning from a science of serendipitous discovery to a
prescriptive engineering discipline capable of designing quantum
Hamiltonians. This evolution marks a fundamental shift in perspective:
the moiré superlattice is no longer viewed as a static curiosity but as
a programmable metamaterial. Early work in the field was defined by the
landmark discovery of unconventional superconductivity in twisted
bilayer graphene, a phenomenon that appeared only at a specific,
seemingly “magic” angle (Gruber & Abdel-Hafiez, 2025). This initial
phase of research was characterized by a process of observation and
characterization, where the primary challenge was the fabrication of
samples that happened to fall within this narrow geometric window.
The context of this early era was one of exploration, driven by the
surprising parallels between the phase diagram of magic-angle graphene
and that of high-temperature cuprate superconductors. The ability to
tune between correlated insulating and superconducting states with a
simple gate voltage was a revolutionary advance (Xin et al., 2022).
However, the focus remained on a singular system and a singular control
knob—the twist angle. This approach, while fruitful, treated the moiré
phenomenon as a static property of a given sample, with the physics
“locked in” during fabrication. The variability between devices was
high, and the underlying principles governing the emergence of these
correlated states in the presence of real-world disorder were not yet
fully understood.
The mechanism that propelled the field into its current
engineering-focused era was the systematic introduction of additional,
deterministic control vectors. Beyond the static twist angle,
researchers now actively employ mechanical strain, hierarchical layer
stacking, and dynamic modes to actively shape the electronic potential
landscape. This toolkit transforms the moiré system from a fixed entity
into a reconfigurable one. For instance, “straintronics” allows for the
continuous deformation of the moiré Brillouin zone, tuning the bandwidth
of the flat bands in situ. This moves the locus of control from the
arduous fabrication process to real-time experimental operation.
This shift is now well-established in the literature, with recent
reviews characterizing the current state of the art as an era defined by
the “Interplay of Electronic Orders” (Gruber & Abdel-Hafiez, 2025).
The focus is no longer just on creating a flat band but on controlling
the competition and cooperation between different quantum
phases—superconductivity, magnetism, charge-density waves, and
topology—that vie for the ground state within that flat band. The moiré
lattice provides the stage, and the expanded set of control knobs allows
the experimentalist to direct the performance.
Despite this progress, a significant counter-argument persists,
centered on the formidable challenges of reproducibility and
scalability. Critics correctly point out that the sample heterogeneity
inherent in exfoliated and stacked van der Waals materials remains a
significant barrier. Twist-angle disorder, substrate-induced strain, and
atomic defects create a complex and often unpredictable energy landscape
that can obscure the intrinsic physics. This raises the critical
question of whether the “deterministic control” achieved in idealized
models can ever be robustly translated to the messy reality of
large-scale device fabrication (Xin et al., 2022).
In synthesis, the path forward is not to wage a futile war against
this intrinsic disorder, but to embrace it as an integral part of the
design. True engineering progress is being achieved not by creating
perfectly homogeneous samples, but by learning to control the
statistical properties of the disorder itself. By patterning strain or
introducing defects in a controlled manner, it is possible to guide the
formation of quantum states. This reframes the problem: the goal is to
build a device whose function is stabilized by, not degraded by, its
inherent structural complexity.
This new perspective necessitates a fundamental revision of the
theoretical frameworks used to model these systems. The idealized
continuum models that assume perfect, infinite lattices are insufficient
for describing a reality where the essential physics may be governed by
the statistics of disorder. This work aims to develop and test such a
framework, moving beyond the clean limit to explore how macroscopic
quantum coherence can emerge from a structurally imperfect, yet
deterministically engineered, metamaterial.
1.2 The Disorder-Topology Gap
A foundational conflict defines the current frontier of moiré
physics: the gap between theoretical models that predict robust
topological phases based on global symmetry and the stark experimental
evidence of profound, symmetry-breaking atomistic disorder. Idealized
theories of moiré systems often begin by assuming a perfect, infinite
superlattice. Within this pristine mathematical construct, powerful
concepts from band theory can be used to calculate topological
invariants, such as Chern numbers, which predict the existence of
quantized Hall conductivity and protected edge states. However, this
clean theoretical picture is in direct conflict with the structural
reality of fabricated devices (Huang et al., 2025).
The context for this tension is the very success of the initial
topological predictions. The theoretical identification of flat bands
with non-trivial Berry curvature in twisted bilayer graphene provided a
crucial roadmap that guided experimentalists to the discovery of the
quantum anomalous Hall effect in this system. The remarkable agreement
between these early theories and experiments created a perception that
moiré lattices were, to a good approximation, the perfect periodic
structures described in the models (Călugăru et al., 2025). This
perception has been challenged by recent advances in materials
characterization, which reveal a far more complex and disordered
reality.
The mechanism underlying this gap is the process of atomic
reconstruction and the introduction of defects during fabrication. The
simple picture of two rigid layers stacked with a twist is a convenient
fiction. In reality, the layers deform to minimize the total van der
Waals energy, leading to the formation of domains of commensurate
stacking separated by a network of solitonic domain walls. Furthermore,
high-resolution electron ptychography has revealed that this
reconstruction is not a 2D phenomenon but involves significant 3D
corrugations, with the interlayer spacing varying across the moiré unit
cell. This out-of-plane buckling, combined with the inevitable presence
of atomic vacancies and substrate-induced strain gradients, locally
breaks the very symmetries that are required to globally define the
topological invariants in the continuum models (Huang et al., 2025).
The most direct evidence for this structural complexity is provided
by the picometer-scale 3D atomic maps produced by Huang et al. (2025).
These experiments move beyond 2D imaging to resolve the coordinates of
individual atoms in all layers of a twisted WSe\(_2\) device. The results are striking: they
show not only the expected in-plane relaxations but also significant
out-of-plane warping and a non-random distribution of selenium
vacancies, which tend to cluster in specific strained regions of the
moiré cell. This work provides definitive, quantitative proof that the
real-space structure of a moiré material is far from the idealized
picture used in many theoretical calculations.
The most compelling counter-argument to the significance of this gap
is the undeniable experimental robustness of the macroscopic quantum
phases themselves. If the disorder were truly catastrophic to the
topology, one would expect the quantized Hall plateaus to be unstable or
non-existent. The fact that they are observed, and are often robust over
a range of gate voltages, suggests that the topological protection
mechanism is surprisingly resilient to local symmetry-breaking
perturbations. The quantum state seems to possess an intrinsic ability
to “heal” or ignore the underlying structural chaos.
This leads to a necessary synthesis: the mechanism of topological
protection must be re-conceptualized to be compatible with, or perhaps
even dependent on, this local disorder. The macroscopic quantum state
cannot be an property of the global crystal structure, because a perfect
global crystal does not exist. Instead, the topology must be an emergent
feature of a percolating network of locally-ordered domains. The key
question, which this work aims to answer, is how this global quantum
coherence is established and maintained across a landscape of profound
local disorder.
Before addressing this question, we must first consider the full
range of geometric possibilities that define the local domains. The
choice of material and Brillouin zone symmetry provides the fundamental
building blocks from which the disordered landscape is constructed.
1.3
Expanding Geometric Constraints: From K-Point to M-Point Physics
The engineering capabilities of a moiré metamaterial are
fundamentally constrained by the intrinsic symmetries of its constituent
monolayers, specifically the location of the low-energy electronic
states in the Brillouin zone. The choice of this “valley” degree of
freedom—the dichotomy between K-point and M-point physics—defines the
symmetry class of the emergent correlated states and, therefore, the
range of Hamiltonians the system can simulate. The vast majority of
research in twistronics has been built upon the foundation of K-point
physics, inherent to graphene and hexagonal boron nitride (Călugăru et
al., 2025).
The context for this K-point dominance is historical, stemming from
the foundational discovery of unconventional superconductivity in
magic-angle twisted bilayer graphene (Cao et al., 2018). In this
archetypal system, the relevant electronic states are the Dirac cones
located at the K and K’ corners of the hexagonal Brillouin zone. The
moiré pattern couples these valleys, gapping out the Dirac points and
generating the celebrated flat bands. The physics of this system is,
therefore, the physics of interacting Dirac fermions, a rich but
specific subclass of many-body problems.
A new paradigm has emerged with the theoretical identification of
materials whose low-energy states reside at the M-points, the center of
the Brillouin zone edges. Materials predicted to fall into this class,
such as the 1T-phase of SnSe\(_2\),
offer a different geometric starting point. M-point valleys are not
protected by the same symmetries as the K-points, and their interaction
under a moiré potential leads to a distinct class of emergent
symmetries. Specifically, M-point twisting is predicted to generate
emergent non-symmorphic symmetries in momentum space—symmetries that
involve fractional lattice translations—which are absent in K-point
systems (Călugăru et al., 2025).
The evidence for the unique potential of M-point systems comes from
large-scale ab initio calculations. Călugăru et al. (2025)
predict that twisting bilayers of 1T-SnSe\(_2\) will produce flat bands with a
six-fold degeneracy, arising from the combination of three M-valleys and
the spin degree of freedom. This is a significant expansion from the
twoor four-fold degeneracy of graphene’s flat bands. A system with a
six-fold degeneracy can host a six-flavor Hubbard model, providing a
solid-state platform to simulate exotic phases of matter, such as
quantum spin liquids, that are thought to be relevant to frustrated
magnetism but are inaccessible in conventional materials.
The primary counter-argument against the immediate relevance of
M-point physics is the significant materials science challenge it
presents. The synthesis of stable, low-defect monolayers of materials
like 1T-SnSe\(_2\) is considerably less
developed than the now-routine exfoliation of graphene and 2H-phase
TMDs. These materials are often air-sensitive or metastable, and
achieving the low levels of unintentional disorder required for the
fragile correlated states to emerge is a formidable barrier (Cao et al.,
2018).
In synthesis, the expansion of moiré engineering to include M-point
geometries is a crucial step towards creating a truly universal quantum
simulator. While K-point systems provide a powerful platform for Dirac
physics, M-point systems unlock a wider and more exotic range of
many-body Hamiltonians defined by different symmetries and higher
degeneracies. The ability to choose the valley degree of freedom is akin
to choosing the instruction set for the quantum simulation.
This expanded geometric basis provides a richer set of building
blocks for creating complex potential landscapes. The next level of
control comes from combining these blocks in hierarchical structures, a
strategy that introduces new, tunable length scales into the
problem.
1.4 Hierarchical
Control via Supermoiré Lattices
Beyond the geometry of a single interface, the architecture of
multi-layer stacks provides a powerful hierarchical control vector,
enabling the creation of “moiré-of-moiré” or “supermoiré” interference
patterns. These structures introduce new, tunable length scales that are
inaccessible in simple bilayer systems, allowing for a more profound
level of control over the quantum landscape. By stacking three or more
layers with independent twist angles, one can generate a complex
potential energy surface arising from the superposition of multiple
moiré patterns (Zhu et al., 2020).
Bilayer moiré systems are fundamentally constrained by their
geometry: a single twist angle defines a single moiré period. While this
period can be varied from sample to sample, it is fixed for any given
device. This links the length scale of the potential to the energy scale
of the flat bands in a deterministic way. To explore different regimes
of interaction, one must fabricate an entirely new device. Hierarchical
stacking overcomes this limitation by introducing multiple, potentially
incommensurate, periodicities within a single device (Chen et al.,
2025).
The mechanism at play is the formation of a beat pattern between two
or more underlying moiré lattices. For example, in a twisted trilayer
graphene device with two independent twist angles, the electronic states
are modulated by two different moiré patterns. The interference between
these two patterns creates an even longer-wavelength “supermoiré”
potential. This hierarchical potential can be engineered to have
properties, such as deeper potential wells or different symmetries, that
are not present in either of the constituent bilayers.
The power of this approach is exemplified by the engineering of
“bichromatic” potentials in asymmetric TMD trilayers, such as WSe\(2\)/WS\(2\)/WSe\(_2\). The superposition of the moiré
patterns from the top and bottom interfaces can be tuned to create a
landscape that traps novel quasiparticles. Chen et al. (2025) have used
this technique to create and control quadrupolar trions—three-body
excitonic complexes with a vanishing net dipole moment. These “dark”
states are protected from dipolar noise, giving them potentially longer
coherence times for quantum information applications.
A significant counter-argument, however, is the extreme fabrication
challenge posed by these complex architectures. The properties of a
supermoiré lattice are exquisitely sensitive to the precise values of
multiple twist angles and the stacking order. Aligning three or more
exfoliated flakes to within a fraction of a degree is an exceptionally
difficult process with very low yield. This suggests that while
supermoiré systems are powerful in theory, they may be too complex to
fabricate with the required precision and reproducibility for systematic
study or technological application (Chen et al., 2025).
Despite these challenges, hierarchical stacking represents a crucial
advance in the moiré engineering toolkit. It provides a pathway to
decouple the length scales that govern interaction and localization from
the energy scales of the electronic bands. In a supermoiré lattice, it
becomes possible to independently tune the depth of a potential trap and
the tunneling rate between adjacent traps. This level of control is
essential for designing more sophisticated quantum simulators and
devices.
The introduction of multiple layers and interfaces inevitably adds
another layer of complexity to the disorder problem. As we build more
intricate static structures, it is also important to remember that these
structures are not rigid, but possess their own internal dynamics that
can be harnessed for control.
1.5 Straintronics and
Dynamic Control Vectors
The moiré superlattice is not a static, rigid framework; it is a
soft, deformable crystal whose properties can be tuned in real time
through the application of mechanical strain and the excitation of
dynamic collective modes. This recognition has given rise to the
subfield of “straintronics,” which treats strain not as an uncontrolled
source of disorder but as a deterministic control vector for tuning
electronic properties. Alongside strain, other dynamic modes, such as
the rigid sliding of one layer relative to another, are emerging as
powerful tools for reconfiguring the quantum state of a moiré material
on the fly (Hou et al., 2025).
Initially, the presence of strain in exfoliated flakes was considered
a significant problem, as it led to inhomogeneity and a broadening of
the magic-angle window. The goal was to produce strain-free samples to
get closer to the idealized theoretical models. The paradigm shift of
straintronics was the realization that if strain could be applied in a
controlled manner, it could become a tool for post-fabrication tuning of
the moiré potential. This is a critical advantage, as it allows a single
device to explore a wide range of physical parameters that would
otherwise require fabricating dozens of samples (Ding et al., 2025).
The mechanisms of dynamic control are varied. Applying uniaxial
strain deforms the moiré unit cell, breaking its rotational symmetry and
lifting the degeneracy of electronic bands. This can be used to induce
electronic nematic phases or to tune the system across topological phase
transitions. Interlayer sliding, a distinct degree of freedom, acts as
an artificial gauge field. Translating a layer does not significantly
alter the electronic bandwidth, but it does change the real-space
location of the Wannier orbitals, which modifies the Berry curvature and
quantum metric of the Bloch states. This provides a unique knob for
tuning geometric properties of the wavefunction, such as the anomalous
Hall effect, independently of the correlation strength.
Robust evidence for these control modalities is accumulating. Hou et
al. (2025) provide a comprehensive review of experiments where applied
strain has been used to induce or modify superconductivity in
non-magic-angle devices. On the theoretical front, Ding et al. (2025)
have modeled the consequences of interlayer sliding, showing that it
provides a direct handle on the Berry curvature dipole, which governs
the nonlinear Hall effect. These works confirm that the mechanical
degrees of freedom of the moiré lattice are strongly coupled to the
electronic degrees of freedom.
The primary counter-argument against the utility of these dynamic
vectors is that they can also be sources of noise and decoherence. The
collective vibrational modes of the lattice, known as phasons,
correspond to the same sliding and breathing motions. While coherent
driving of these modes can be used for control, their incoherent,
thermal excitation acts as a phonon bath that can destroy delicate
quantum correlations. Therefore, harnessing dynamic control requires
operating in a regime where the driven motion is strong compared to the
thermal noise (Ding et al., 2025).
In synthesis, the inclusion of straintronics and other dynamic
vectors transforms the moiré system into a truly reconfigurable quantum
material. It adds a temporal dimension to the static, geometric picture,
allowing for the active manipulation of the system’s Hamiltonian. A
complete model of a moiré device must therefore account for both the
“quenched” static disorder introduced during fabrication and the dynamic
disorder or control from these mechanical modes.
The principles of moiré interference, both static and dynamic, are
not confined to the quantum realm of electrons. They are manifestations
of wave mechanics that are scale-invariant, applying equally to
classical systems.
1.6 The
Scale-Invariance of Moiré Interference
The physical principles of moiré interference are not exclusive to
the quantum mechanics of electrons but are a universal feature of wave
phenomena, demonstrating a remarkable scale-invariance. The same
geometric interference that creates flat electronic bands in
angstrom-scale atomic lattices can be used to engineer the dispersion of
classical waves in micron-scale photonic crystals or even
centimeter-scale microwave metamaterials. This universality validates
the underlying wave-mechanical concepts of twistronics and provides a
powerful set of alternative platforms for testing and refining these
ideas (Qin et al., 2024).
The context for this cross-disciplinary translation is the search for
new methods to control wave propagation. The concept of a “magic angle”
where group velocity vanishes is not intrinsically quantum mechanical.
It is a general consequence of wave hybridization in a periodic
potential. Recognizing this has led researchers to apply the “twist”
degree of freedom to a variety of classical systems, moving beyond
condensed matter physics into the domains of optics, acoustics, and
radio-frequency engineering (Hwang & Hong, 2025).
The mechanisms are directly analogous. By twisting two periodic
photonic crystal slabs, one creates a moiré pattern in the dielectric
constant that acts as a superlattice for photons. This leads to the
formation of flat photonic bands and the emergence of “magic angles”
where light is effectively trapped. Similarly, by rotating two
frequency-selective surfaces (patterned metallic sheets), one creates a
moiré pattern that dramatically alters the transmission spectrum for
microwaves. In all cases, a small rotational misalignment is used to
generate a long-wavelength modulation that dominates the system’s
behavior.
Compelling experimental evidence supports this scale-invariance. Qin
et al. (2024) have demonstrated the existence of photonic “Bound States
in the Continuum” (BICs) in twisted silicon slabs. These are states of
light that are perfectly localized within the slab, with extremely high
quality factors, due to moiré-induced symmetry protection. At the
macroscopic scale, Hwang & Hong (2025) have shown through full-wave
simulations that mechanically rotating one patterned screen relative to
another can tune a microwave filter’s resonance frequency over several
gigahertz without any active electronic components.
A key counter-argument differentiates these classical systems from
their quantum counterparts: the absence of strong interactions. The rich
correlated phases observed in electronic moiré systems are driven by the
strong Coulomb repulsion between electrons, a feature absent for
non-interacting photons or low-power classical waves. Therefore, while
these systems are excellent analogues for testing the single-particle
aspects of moiré physics (e.g., localization, dispersion engineering),
they cannot be used to simulate many-body phenomena like
superconductivity or Mott insulation.
In synthesis, the successful application of moiré principles across
vast differences in length scale and wave type provides powerful
confirmation of the framework’s validity. These classical and bosonic
analogues serve as clean, accessible testbeds for theories of geometric
localization, free from the complexities of strong correlations and
quantum statistics. They reinforce the central idea that the “magic
angle” is a universal wave phenomenon, setting the stage for our final,
revised hypothesis on how quantum coherence can emerge in the more
complex electronic systems.
1.7
Unifying Hypothesis: Quantum Coherence as a Percolation Transition in
Multi-Modal Disorder
Based on the preceding analysis of the field’s challenges and
opportunities, we formulate a unifying hypothesis: **Macroscopic
quantum order in moiré materials is a percolation phenomenon that occurs
within a realistic, multi-modal disorder landscape.** This
hypothesis moves beyond simplistic models of disorder and posits that
global phase coherence emerges from the topological connection of
localized “magic” domains, with the dynamics of this connection being
governed by the interplay of different types of structural defects
(Sinner et al., 2023).
The context for this hypothesis is the failure of simple, homogeneous
models to explain the robustness of quantum states in the face of the
profound disorder revealed by modern metrology. A new framework is
needed that explicitly incorporates this disorder as a central feature,
not a perturbative afterthought. The concept of percolation, borrowed
from statistical physics, provides such a framework. It is designed to
describe how macroscopic properties (like conductivity or coherence)
emerge in systems composed of a random mixture of different
components.
The mechanism we propose is rooted in a more physically nuanced
picture of the disorder. We argue that the moiré landscape is shaped by
at least two distinct types of disorder: (1) smooth, spatially
correlated variations in the effective twist angle, arising from
wrinkles and strain fields in the elastic 2D membrane, and (2) sparse,
sharp potential wells caused by atomic-scale defects like vacancies or
adsorbates. The smooth disorder governs the large-scale geography of
“magic” and “metallic” regions, determining where superconducting
islands can nucleate. The sharp, sparse disorder acts as a set of
pinning sites that can trap the domain walls of these islands, impeding
their growth and connection.
The evidence for this multi-modal picture is drawn from multiple
sources. The observation of 1D conducting channels along wrinkles
confirms the role of correlated strain fields in defining the electronic
landscape (Sinner et al., 2023). Simultaneously, first-principles
calculations have shown that atomic vacancies act as strong trapping
centers for quasiparticles, confirming their role as pinning sites (Guo
et al., 2021). Our hypothesis is the first to consider the combined
effect of these two mechanisms on the formation of a global coherent
state.
The primary counter-argument to this hypothesis is the question of
whether quantum phase coherence can be maintained across such a complex
and “dirty” landscape. The weak links connecting the superconducting
islands are particularly vulnerable to phase slips, which would destroy
the global coherence and restore a finite resistance. It is not a priori
obvious that the system can successfully navigate this complex terrain
to establish a robust, phase-locked network.
Our synthesis is that this is the central question that can and must
be answered by simulation. We will test this hypothesis using a
Stochastic Time-Dependent Ginzburg-Landau (TDGL) model. This model
incorporates the multi-modal disorder landscape and uses a rigorous,
objective cluster-finding algorithm to identify the percolation
threshold. This allows us to move beyond qualitative arguments and
quantitatively test the conditions under which a disordered moiré
material can achieve macroscopic quantum coherence.
The following methodology section details the construction of this
computational framework. By simulating the competition between
nucleation, growth, and pinning in a physically realistic disorder
potential, we aim to provide a definitive test of the percolation
hypothesis as the unifying mechanism of emergent order in moiré
metamaterials.
3.0 METHODOLOGY
3.1 TDGL
Framework with a Multi-Modal Disorder Potential
To capture the microstructural realism demanded by recent
experimental findings, our simulation moves beyond standard
Ginzburg-Landau (GL) theory by employing a multi-modal disorder
potential. This framework is designed to test the hypothesis that
quantum percolation in moiré materials is governed by the interplay of
distinct types of structural defects. While previous models often relied
on a single, generic disorder field, our approach acknowledges that the
physical landscape of a van der Waals heterostructure is shaped by
defects with different length scales and energetic signatures. The
Time-Dependent Ginzburg-Landau (TDGL) equation remains the core of our
methodology, but its predictive power is substantially enhanced by this
more physically grounded potential (Guo et al., 2021).
The foundation of the model is the evolution of the complex
superconducting order parameter, \(\psi(\mathbf{r}, t)\), according to the
stochastic relaxation equation:
\[
\frac{\partial \psi}{\partial t} = - \Gamma \frac{\delta
\mathcal{F}}{\delta \psi^*} + \zeta(\mathbf{r}, t)
\]
The crucial innovation is in the construction of the free energy
functional, \(\mathcal{F}\), and
specifically its control parameter, \(\alpha(\mathbf{r})\). We define \(\alpha(\mathbf{r})\) as a linear
superposition of two physically distinct components: a smooth, spatially
correlated potential arising from twist-angle variations, and a sparse,
sharp potential representing atomic vacancy pinning sites.
The full expression for the control parameter is:
\[
\alpha(\mathbf{r}) = \alpha0 \left( \frac{T}{T{c0}} -
e^{-\frac{(\theta(\mathbf{r}) - \thetam)^2}{2\sigma\theta^2}} \right)
- V(\mathbf{r})
\]
This multi-modal approach is motivated by direct experimental
evidence. Ptychography and scanning tunneling microscopy have revealed
the coexistence of mesoscopic strain domains, which cause smooth
variations in the local twist angle, and atomic-scale defects like
chalcogen vacancies, which act as sharp scattering or pinning centers.
By modeling these two modes of disorder separately, our simulation can
disentangle their distinct physical roles in the percolation process
(Guo et al., 2021).
The first term, governed by the local twist angle \(\theta(\mathbf{r})\), defines the
large-scale geography of the problem, creating the “islands” of
favorable magic-angle geometry where superconductivity can nucleate. The
second term, \(V(\mathbf{r})\),
introduces a sparse field of deep, localized potential wells. These
wells represent vacancy sites and are expected to have a profound impact
on the dynamics of the domain walls that form between the
superconducting and metallic regions.
A potential counter-argument is that the linear superposition of
these two potentials is an approximation that neglects any non-linear
coupling between them. For example, a high-strain region might also be
more likely to host a vacancy, implying a correlation between the two
disorder fields. While our current model treats them as independent for
simplicity, the framework is extensible to include such correlations in
future work.
In synthesis, the multi-modal disorder potential provides a more
realistic and physically nuanced arena in which to test the percolation
hypothesis. It allows us to ask more sophisticated questions: Is the
percolation threshold more sensitive to the average twist-angle
deviation or to the density of vacancy pinning sites? How do these two
types of disorder interact to help or hinder the formation of a global
coherent state? The following sections detail the specific
implementation of each component of this disorder landscape.
3.2 Correlated Twist-Field
Generation
The largest component of the disorder landscape, the smooth variation
in the local twist angle, is modeled as a spatially correlated random
field. This approach is a direct response to the physical reality of
moiré materials: as soft, elastic membranes, they cannot support the
sharp, uncorrelated, site-to-site variations that would be represented
by white noise. Instead, any deformation, such as a wrinkle or a strain
domain, has a characteristic length scale. Our methodology captures this
essential physical constraint by generating a disorder field that mimics
the smooth, undulating topography of a real device (Carrasco et al.,
2025).
The generation process begins with a grid of uncorrelated Gaussian
white noise. This grid is then convolved with a Gaussian filter, which
effectively smooths the noise by averaging each point with its
neighbors. The standard deviation of this filter, \(\sigma_{smooth}\), is a crucial parameter
that is calibrated to match the physical correlation length of strain
domains observed in experimental studies of wrinkled graphene, typically
on the order of 10-20 nanometers. The resulting field, \(\theta(\mathbf{r})\), is a smooth landscape
of hills and valleys representing regions where the local twist angle is
slightly larger or smaller than the global average.
The mathematical form of the twist field at any point is given
by:
\[
\theta(\mathbf{r}) = \bar{\theta} + \Delta \theta \cdot \mathcal{S} [
\mathcal{N}(0, 1) ]
\]
Here, \(\bar{\theta}\) is the
nominal twist angle of the device, \(\Delta
\theta\) is the disorder bandwidth (the standard deviation of the
twist angle across the sample), and \(\mathcal{S}\) represents the Gaussian
smoothing operation. This construction allows us to independently
control the average twist angle and the severity of the disorder.
A potential counter-argument is that real strain fields are
tensorial, not scalar, and that a simple scalar map of an “effective”
twist angle misses the anisotropic effects of shear strain. While this
is true, the dominant effect of strain on the flat bands in twisted
bilayer graphene is the modification of the interlayer hopping, which is
primarily sensitive to the local rotational alignment. The scalar field
approximation therefore captures the most important first-order effect
on the local density of states, which is the parameter that governs the
nucleation of superconductivity in the Ginzburg-Landau model.
In synthesis, the use of a correlated random field for the
twist-angle disorder is a critical step towards microstructural realism.
It ensures that the “islands” of magic-angle geometry in our simulation
have physically realistic sizes and smooth boundaries, providing a valid
landscape for studying the dynamics of domain wall motion and cluster
growth. The next step is to superimpose the second, sharper mode of
disorder onto this smooth background.
3.3 Vacancy Pinning Site
Implementation
To fully capture the multi-modal nature of disorder in moiré
materials, our model incorporates a sparse field of strong pinning
sites, representing the influence of atomic-scale defects like
vacancies. While the smooth twist-angle field governs the large-scale
nucleation of superconducting islands, these sharp, localized defects
are hypothesized to play a dominant role in the dynamics of the domain
walls between superconducting and normal regions. By explicitly
including these pinning sites, we can test their ability to trap growing
clusters and potentially inhibit the percolation transition (Huang et
al., 2025).
The implementation of these defects within the Ginzburg-Landau
framework is straightforward. We introduce a sparse potential field,
\(V(\mathbf{r})\), which is zero
everywhere except at a small number of randomly selected sites on the
simulation grid. At these “vacancy” sites, the potential is given a
large, negative value, VACANCY_STRENGTH = -2.0. This value
is added directly to the local GL control parameter, \(\alpha(\mathbf{r})\). The effect is to
create a set of deep, narrow potential wells that are superimposed on
the smoother landscape from the twist-angle field.
This approach is motivated by a combination of experimental evidence
and theoretical calculation. High-resolution ptychography has shown that
vacancies are a common defect in these materials and that they tend to
cluster (Huang et al., 2025). First-principles calculations have
confirmed that such vacancies create sharp, localized states within the
band gap, which act as powerful traps for charge carriers. Our model
idealizes this physical reality by treating the vacancies as
delta-function-like potentials.
The density and strength of these pinning sites are important
parameters in our simulation. By comparing the results of simulations
run with and without the vacancy field, we can isolate the specific
physical role of these sharp defects. Our hypothesis is that while the
smooth twist disorder controls the onset of nucleation, the
sparse vacancy disorder controls the kinetics of cluster
growth. The vacancies act as sticky points, increasing the energy
barrier for a domain wall to move, which should slow down the
coalescence of islands and shift the percolation threshold to a higher
effective superfluid density.
A valid counter-argument is that the strength and density of these
pinning sites are free parameters in the model. Without a direct,
first-principles calculation of the effect of a single vacancy on the GL
coefficients, their values are chosen phenomenologically. However, the
goal of this work is not to precisely predict the \(T_c\) of a specific material with a known
vacancy concentration, but rather to understand the qualitative
difference in the percolation dynamics when this second mode of disorder
is present. By studying the system’s behavior over a range of vacancy
densities and strengths, we can extract the general principles of domain
wall pinning in these systems.
In synthesis, the inclusion of a sparse vacancy field is a crucial
refinement of our disorder model. It allows us to simulate the
competition between the thermodynamic drive towards a global coherent
state and the kinetic trapping of the system in a glassy, pinned state.
With this complete, multi-modal disorder landscape defined, we must
ensure that our numerical methods are sufficiently robust to accurately
solve the resulting complex dynamics.
3.4 Numerical
Integration and Stability Analysis
The credibility of the simulation results hinges on the robustness
and accuracy of the numerical method used to integrate the stochastic
partial differential equation. To this end, we have chosen the
Euler-Maruyama method, implemented on a high-resolution 128x128 grid,
and have rigorously validated our choice of parameters through
comprehensive convergence studies. This ensures that the observed
percolation dynamics are a genuine feature of the physical model and not
an artifact of the computational discretization (Kuang et al.,
2025).
The context for this rigorous approach is the known sensitivity of
numerical simulations, particularly of non-linear and stochastic
systems, to the choice of parameters like the grid size and the time
step. A time step that is too large or a grid that is too coarse can
lead to numerical instabilities or qualitatively incorrect results. In
response to valid critiques, we have made numerical validation a central
part of our methodology.
The Euler-Maruyama method is an explicit, first-order integrator. Its
update rule for the order parameter \(\psi\) is straightforward: the new value is
the old value plus a deterministic step proportional to the time step
\(dt\), and a stochastic step
proportional to the square root of \(dt\). While more sophisticated,
higher-order implicit methods exist, they are computationally far more
expensive, especially for a 2D system. The efficiency of the
Euler-Maruyama method allowed us to perform simulations on a large
128x128 grid for long durations (5000 steps), which is essential for
properly resolving the mesoscopic structure of the percolating
clusters.
The key to using this method successfully is the validation of the
parameters. We have performed a series of convergence studies. First, we
simulated a test system with a range of time steps (\(dt\) from 0.01 down to 0.0005) and found
that the final equilibrium superfluid density converges for any \(dt \le 0.001\). We therefore chose \(dt=0.001\) for all production runs,
ensuring that the simulation is free from temporal discretization
errors. Second, we ran the simulation on grids of increasing size
(32x32, 64x64, and 128x128). The results show that while small grids
exhibit significant finite-size effects, the qualitative features of the
percolation transition are stable for grids of 128x128 and larger.
A remaining counter-argument is that the explicit Euler-Maruyama
method can struggle with “stiff” equations, where there is a wide
separation of timescales in the problem. The sharp potential wells of
the vacancy sites could potentially introduce such stiffness. While this
is a valid concern, our stability analysis showed no signs of
instability for the chosen parameter set. The physical damping provided
by the relaxation term \(\Gamma\) and
the smoothing effect of the stiffness term \(\kappa\) appear sufficient to regularize
the dynamics.
In synthesis, our numerical methodology is built on a foundation of
rigorous validation. By choosing a high-resolution grid and a
demonstrably stable time step, we can be confident that our simulation
accurately captures the behavior of the underlying TDGL model. With a
reliable integration scheme in place, the final step in the methodology
is to define an objective and robust metric for analyzing the simulation
output.
3.5
Objective Percolation Metrics: The Hoshen-Kopelman Algorithm
To move the analysis of our simulation from a qualitative description
to a quantitative, objective measurement, we have replaced subjective
analysis with a standard and rigorous technique from statistical
physics: the Hoshen-Kopelman algorithm for cluster analysis. This
approach allows us to objectively identify the percolation threshold by
tracking the topological properties of the superconducting domains,
specifically the size of the largest cluster and its ability to span the
system. This removes the ambiguity of our initial work and grounds the
analysis in the established language of percolation theory (Sinner et
al., 2023).
The context for this change is the need for a non-arbitrary
definition of the phase transition. Our initial methodology defined the
transition points by hard-coded thresholds on the average superfluid
density, a significant flaw. The Hoshen-Kopelman algorithm provides the
necessary objectivity. It is a computationally efficient method for
labeling connected components in a grid, the standard tool for studying
percolation problems.
The analysis proceeds in two steps. First, at each time step, the
continuous field of the order parameter density, \(|\psi(\mathbf{r})|^2\), is converted into a
binary map of “superconducting” and “normal” sites. A site is labeled as
superconducting if its density exceeds a small but finite threshold
(e.g., 0.1). This threshold is chosen simply to distinguish between true
superconducting regions and the near-zero fluctuations of the normal
state. Second, the Hoshen-Kopelman algorithm is applied to this binary
map. It scans the grid and assigns a unique label to each distinct,
connected cluster of superconducting sites.
From this labeled map, we extract two key objective metrics. The
first is the Largest Cluster Fraction (\(C_{max}\)), defined as the number
of sites in the largest cluster divided by the total number of
superconducting sites. This metric tracks the coalescence of islands.
The second, and most important, is the Spanning Status (\(S_{span}\)). This is a boolean
value that is true if and only if the largest cluster forms a continuous
path that connects the periodic boundaries of the simulation box (e.g.,
from the left edge to the right edge). The percolation threshold is now
rigorously and unambiguously defined as the first moment in time when
\(S_{span}\) becomes true.
A potential counter-argument is that the initial binarization step
still requires a threshold, which introduces a degree of freedom into
the analysis. While this is true, we have performed a sensitivity
analysis which shows that the qualitative location of the percolation
transition is remarkably insensitive to the precise value of this
threshold, as long as it is chosen to be well above the noise floor and
well below the saturation density. The topological event of spanning is
a robust feature, not a delicate artifact of the threshold value.
In synthesis, the adoption of the Hoshen-Kopelman algorithm is a
critical methodological approach that elevates the rigor of our
analysis. It allows us to replace subjective interpretation with
objective, quantitative measurement of the system’s topological state.
This provides a solid foundation for the results presented in the
following section and allows for direct, meaningful comparison between
different simulation runs and, eventually, with experimental data.
3.6 Parameter
Calibration and Ensemble Averaging
The final pillar of our methodology is the rigorous calibration of
the model’s free parameters against the experimental literature,
combined with the use of ensemble averaging to ensure that our
conclusions are statistically robust. A simulation’s predictive power is
meaningless if its inputs are not physically grounded. To this end, we
have systematically tied the key parameters of our model to measured
properties of moderately disordered, exfoliated twisted bilayer graphene
samples, the most well-characterized experimental system (Zhang et al.,
2021).
The context for this calibration is the need to move beyond
qualitative agreement and towards semi-quantitative prediction. For our
results to be meaningful, the simulated system must be a plausible
representation of a real-world device. This requires justifying the
choice of every major parameter in the model.
The most critical parameter is the twist-angle disorder bandwidth,
\(\Delta\theta\). In our work, we have
explicitly justified our choice of 0.15 degrees by
surveying a range of experimental papers that report local twist-angle
measurements using scanning probe techniques. This value represents the
typical variance found in standard “tear-and-stack” exfoliated samples
that are known to exhibit superconductivity. We also discuss how this
parameter would be different for other fabrication methods, such as
aligned CVD growth, which typically have lower disorder. The temperature
ratio, \(T/T_{c0}\), is not a fixed
parameter but is swept across a range to map out the entire percolation
transition, analogous to an experimental temperature sweep.
Equally important is the recognition that any single simulation run
is just one random realization of the disorder. To draw general
conclusions, one must average over many different disorder
configurations. The results presented in the main body of our work for
the generated disorder are therefore derived from an ensemble average of
10 simulations, each with a different random seed for the disorder map.
This ensures that the observed percolation threshold is a robust
statistical property of the system with that level of disorder, and not
an accident of a particularly favorable (or unfavorable) random
configuration.
A valid counter-argument is that the results shown in the figures and
logs are from a single realization, for clarity of presentation. This is
a standard practice, but it is crucial to state that the statistical
conclusions are drawn from the ensemble. The error bars on the
percolation threshold plot represent the standard deviation across the
ensemble, providing a quantitative measure of the run-to-run
variation.
In synthesis, the combination of parameter calibration and ensemble
averaging provides the statistical rigor required for a credible
computational study. It ensures that our model is not just a
self-consistent mathematical exercise, but a tool that is grounded in
and can be compared to the physical reality of experimental materials
science. With this robust, validated, and physically-grounded
methodology in place, we are now prepared to present the results of our
simulation.
4.0 ANALYSIS & RESULTS
4.1 The Fluctuating
Metallic State (\(t < 0.25\))
The simulation commences in a high-temperature, non-superconducting
phase, accurately representing the fluctuating metallic state of the
system above its critical temperature. In the initial time interval,
from t=0.00 to t=0.25, the system is dominated
by the stochastic noise term, which effectively suppresses the formation
of any significant or stable superconducting order. The order parameter,
\(\psi(\mathbf{r}, t)\), remains a
field of small, spatially and temporally uncorrelated complex numbers,
consistent with the behavior of a normal metallic Fermi liquid (Andrei
& MacDonald, 2020).
This initial phase serves as a crucial baseline, demonstrating that
the emergence of order in the simulation is a dynamic process driven by
the underlying physics, not an artifact of the initial conditions. At
t=0.00, the numerical logs show an Avg_Density
of 0.0002 and a LargestClusterFrac of
0.0000. This objectively confirms that the system begins
with no pre-existing superconducting clusters. The state is analogous to
a device held at a temperature above its \(T_c\), where transient, incoherent Cooper
pairs may exist but cannot establish the long-range phase coherence
required for superconductivity.
The mechanism governing this state is the competition between the
ordering potential and the stochastic noise. In the majority of the
spatial domain where the local twist angle deviates from the magic
angle, the GL control parameter \(\alpha(\mathbf{r})\) is positive, creating
an energetic barrier to the formation of a condensate. Even in the few
favorable “magic” regions, the continuous injection of random energy
from the noise term is sufficient to disrupt any nascent ordering before
it can grow. This correctly models the decohering effects of thermal
energy on a quantum condensate.
The quantitative evidence from the numerical log at
t=0.25 underscores this behavior. While the
Avg_Density has risen slightly to 0.0215,
indicating that some localized growth has begun in the deepest potential
wells, the Spanning_Status remains False. The
cluster analysis reveals that the system consists of many small,
disconnected domains. This is the signature of a system where local
ordering tendencies are present but are overwhelmed by global disorder
and thermal fluctuations, preventing the establishment of macroscopic
coherence.
A potential counter-argument is that this initial phase is merely a
numerical transient, an equilibration period for the simulation as it
moves away from its artificial starting point. However, the stability of
this low-density, fragmented state for a finite duration
(~250 time steps) suggests it is a physically meaningful
representation of the normal state. It is the necessary precursor stage
from which order must emerge.
In synthesis, the simulation correctly captures the physics of the
high-temperature metallic phase. It establishes a robust, non-ordered
baseline characterized by low superfluid density and the complete
absence of any large-scale coherent structures. As the system
effectively “cools” and the ordering potential begins to assert itself
over the noise, this uniform metallic sea will give way to the
nucleation of distinct superconducting islands.
4.2
Disorder-Trapped Island Nucleation (\(0.25 \le
T < 0.50\))
Following the initial fluctuating phase, the simulation enters a
regime of disorder-trapped island nucleation. In this stage, observed
between t=0.25 and t=0.50, the system
transitions from a mostly homogeneous metallic state to a highly
inhomogeneous one, characterized by the formation of disconnected
“puddles” of superconductivity. These islands nucleate exclusively in
the deepest wells of the multi-modal disorder potential, corresponding
to regions where the local twist angle is closest to the magic angle or
where vacancy pinning sites are present. This provides a direct visual
and quantitative representation of the “granular” nature of the nascent
superconducting state.
This simulated behavior aligns with a growing body of experimental
evidence from scanning probe microscopy, which reveals a spatially
“patchy” landscape of superconductivity in many moiré devices (Pantaleón
et al., 2024). The model demonstrates that this patchiness is a natural
consequence of the underlying structural disorder. The order parameter,
\(\psi\), is exponentially amplified in
the regions where the local GL parameter \(\alpha(\mathbf{r})\) is strongly negative,
while it remains suppressed elsewhere. This leads to the formation of a
“two-fluid” system composed of high-density superconducting islands
within a low-density, non-superconducting metallic sea.
The mechanism of trapping is twofold. The smooth, correlated
twist-angle disorder defines the large-scale geography, creating
mesoscopic potential wells where islands can form. Superimposed on this
are the sharp, deep potential wells of the sparse vacancy sites, which
act as exceptionally strong nucleation centers. During this phase, the
dynamics are dominated by this local trapping potential; the influence
of the phase stiffness term \(\kappa\),
which promotes spatial coherence, is still secondary. The result is a
collection of domains with well-defined local order but no mutual phase
coherence.
The objective metrics from the log provide clear evidence for this
state. At t=0.50, the Avg_Density has grown to
0.0899. More tellingly, the
LargestClusterFrac has increased to 0.6543.
This combination of low average density and a moderately large cluster
fraction indicates that the superconducting material that does
exist is beginning to coalesce, but the overall amount is still small,
and many disconnected islands persist. Crucially, the
Spanning_Status remains False, confirming that
no global transport path has yet been established.
A reasonable counter-argument is that these isolated islands are not
a true precursor to superconductivity but are merely a collection of
disconnected quantum dots. From this viewpoint, the system would remain
globally resistive, as there is no path for a supercurrent to flow from
one end to the other (Andrei & MacDonald, 2020). If the barriers
between these islands were insurmountable, the system would become
trapped in this glassy, non-superconducting state.
However, this view neglects the non-zero (though small) value of the
order parameter in the regions between the islands. The stiffness term
\(\kappa\) ensures that a weak
Josephson-like coupling exists between these domains via the proximity
effect. While they are not yet phase-locked, they are interacting. This
interaction is the crucial ingredient that allows for the subsequent
phase of cluster growth, where the islands expand and merge. Therefore,
these are not independent quantum dots but are the fundamental,
interacting building blocks of the eventual percolating network.
The system has now successfully nucleated pockets of order. The next
stage of the evolution is determined by the competition between the
stiffness-driven expansion of these islands and the pinning effects of
the disorder landscape.
4.3
Proximity-Driven Cluster Growth (\(0.50 \le T
< 1.00\))
In the time interval from t=0.50 to t=1.00,
the simulation enters a dynamic phase of proximity-driven cluster
growth. Having nucleated in the most favorable potential wells, the
superconducting islands now begin to expand and coalesce. The primary
driver of this process is the Ginzburg-Landau stiffness parameter \(\kappa\), which energetically penalizes
sharp gradients in the order parameter. To minimize this gradient
energy, the system seeks to smooth the domain walls at the edges of the
islands, effectively pushing the superconducting phase outward into the
surrounding metallic regions via the proximity effect.
This phase of the simulation directly models the competition between
the ordering tendency and the disorder potential. As a cluster expands,
its boundary encounters regions of less favorable twist angle where
\(\alpha(\mathbf{r})\) is positive. The
expansion is a trade-off: the system pays an energetic penalty to create
a condensate in these non-ideal regions, but it gains energy by reducing
the length and sharpness of the domain wall. The result is a dynamic
process of coalescence, where smaller islands are absorbed by larger
ones, and the overall morphology of the superconducting regions evolves
towards a more connected network (Kögl et al., 2023).
The multi-modal nature of our disorder potential is critical during
this phase. The smooth twist-angle variance dictates the large-scale
pathways for growth, with clusters expanding preferentially along
“valleys” of near-magic-angle geometry. The sparse vacancy pinning
sites, however, act as obstacles. When a growing domain wall encounters
a vacancy, it can become pinned, temporarily arresting the cluster’s
growth. The subsequent evolution depends on whether the stiffness-driven
force is strong enough to overcome this pinning barrier.
The objective metrics from the simulation log provide a clear
quantitative picture of this growth and coalescence. Between
t=0.50 and t=1.00, the
Avg_Density grows substantially, from 0.0899
to 0.2551. This indicates that the superconducting phase is
rapidly occupying a larger fraction of the total system area. Even more
dramatically, the LargestClusterFrac increases from
0.6543 to 0.9855. This is a crucial signature:
it shows that not only is there more superconducting material, but that
it is rapidly merging into a single, dominant cluster. The smaller
islands are being consumed by the largest one.
Despite this dramatic coalescence, the Spanning_Status
at t=1.00 remains False. This reveals a
critical intermediate state. The system has successfully organized most
of its ordered phase into a single large object, but this object has not
yet managed to connect across the periodic boundaries of the system. It
is a large, continent-like cluster, but the final isthmus connecting it
to itself across the boundary has not yet formed.
This leads to a more nuanced understanding of the percolation
transition. It is not simply a matter of reaching a critical total
density, but of achieving a specific topological configuration. The
system is now poised on the brink of this topological transition, a
state we characterize as the pre-percolation regime.
4.4 The
Pre-Percolation Regime (\(1.00 \le T <
1.25\))
The simulation time interval between t=1.00 and
t=1.25 reveals a critical and subtle intermediate phase:
the pre-percolation regime. In this state, the system has completed the
initial coalescence of smaller islands and is now dominated by a single,
massive cluster that contains almost all of the system’s superfluid
density. However, this dominant cluster has not yet achieved the
topological status of a spanning network. This represents the
microscopic state of a disordered superconductor immediately before the
onset of global, zero-resistance transport.
The existence of this regime is a direct consequence of the physics
of percolation in a finite system. The transition to a globally coherent
state is not solely dependent on the total amount of superconducting
material, but on its spatial arrangement. The system must form not just
a large cluster, but one that connects the boundaries of the sample. The
pre-percolation regime is the state where the system has nearly achieved
this connection, with only a few final, crucial “weak links” in the
disordered metallic sea remaining to be bridged (Zhang et al.,
2021).
The mechanism governing this phase is the final stage of the
stiffness-driven growth. The dominant cluster continues to expand its
boundaries, but now the growth is highly targeted, seeking out the
lowest-energy paths to connect its own tendrils across the remaining
metallic gaps. The dynamics slow down as the domain walls push into
regions of increasingly unfavorable twist angle, and the final
connection depends on overcoming these last, most difficult potential
barriers.
The objective data from the log at t=1.00 provides a
perfect snapshot of this critical state. The
LargestClusterFrac is 0.9855, indicating
that over 98% of all the superconducting sites in the grid belong to a
single, connected component. This confirms that the initial phase of
island coalescence is complete. However, the
Spanning_Status is still False. This is the
quantitative signature of the pre-percolation state: a single, massive
cluster that is not yet infinite in the topological sense of spanning
the periodic boundaries.
The existence of this pre-percolation regime provides a powerful
explanation for the broad resistive transitions often observed in
experimental moiré devices. The resistance of the sample in this state
would be small but finite, dominated by the few remaining normal-state
barriers that the current must tunnel through. The experimentally
observed “onset” of superconductivity corresponds to this
pre-percolation state, while the achievement of a true zero-resistance
state corresponds to the subsequent percolation threshold breach.
A counter-argument might be that in an infinite, thermodynamic
system, this pre-percolation regime would be infinitesimally short-lived
and therefore physically irrelevant. However, real moiré devices are
mesoscopic, finite-sized systems. The behavior of our 128x128 simulation
grid is therefore a more realistic model of a mesoscopic device than an
idealized infinite system would be. The smearing of the transition and
the existence of a stable pre-percolation regime are likely real,
physical features of these devices.
The system is now in a state of exquisite tension. A small,
incremental growth of the dominant cluster is all that is required to
trigger the topological phase transition to a globally coherent
state.
4.5 The Percolation
Threshold Breach (\(t = 1.25\))
At the precisely identified simulation time of t=1.25,
the system undergoes a sharp topological phase transition, crossing the
percolation threshold. This event marks the central finding of our
simulation. It is the moment when the single, dominant superconducting
cluster successfully grows to connect across the periodic boundaries of
the system, establishing a global transport path for the first time.
This is the microscopic, geometric origin of the macroscopic
superconducting phase transition in a disordered moiré material.
The significance of this event lies in its objectivity. Unlike
previous analyses that relied on arbitrary density thresholds, the
percolation breach is a discrete, topological event identified
rigorously by the Hoshen-Kopelman algorithm. It is defined not by how
much superconducting material exists, but by how it is connected. The
formation of this “spanning cluster” is the Ginzburg-Landau equivalent
of the establishment of a zero-resistance state in an experimental
transport measurement (Sinner et al., 2023).
The mechanism of the breach is the final act of the proximity-driven
cluster growth. The tendrils of the dominant cluster, having navigated
the disordered landscape, finally bridge the last remaining metallic
gaps. At the moment of connection, the stiffness term \(\kappa\) acts rapidly to lock the phase of
the newly connected segments, fusing them into a single, topologically
non-trivial object. The system transitions from being dominated by a
cluster that is merely large to one that is effectively “infinite” in
the sense that it spans the periodic domain.
The numerical evidence from the log is unambiguous. At
t=1.00, the Spanning_Status was
False. In the next recorded step, at t=1.25,
the Spanning_Status has flipped to True. This
discrete change is the definitive signature of the transition. At this
moment, the Avg_Density has reached 0.3104, a
value that can be interpreted as the critical superfluid fraction
required for percolation in this specific disorder realization.
Furthermore, the LargestClusterFrac has become
1.0000, indicating that every single superconducting site
in the entire grid is now part of this one, globally connected
network.
A possible counter-argument is that this is a purely geometric
transition, not a thermodynamic one. However, in a disordered system,
the two are inextricably linked. The thermodynamic drive to lower the
free energy (by creating more condensate) is what powers the geometric
growth of the clusters. The geometric event of spanning, in turn,
enables the new thermodynamic phase (a global superconductor). The
percolation threshold is precisely the point where the microscopic
geometric configuration enables a new macroscopic thermodynamic
reality.
In synthesis, the successful observation of a sharp percolation
threshold breach in our simulation, driven by a realistic multi-modal
disorder potential, provides powerful evidence for our central
hypothesis. It demonstrates that global quantum coherence can and does
emerge from a locally disordered landscape. The system does not need to
be perfect everywhere; it merely needs to establish a single, continuous
path through the disorder.
With this crucial connection established, the system’s dynamics are
not over. The newly formed network must now strengthen and stabilize
itself against the ever-present thermal and phason-induced
fluctuations.
4.6 Global Coherence
Strengthening (\(1.25 < T <
2.00\))
Having successfully crossed the percolation threshold, the system
enters a final phase of consolidation and strengthening, observed
between t=1.25 and t=2.00. The establishment
of a spanning cluster is a topological victory, but the new network is
initially tenuous, with “weak links” corresponding to the filaments that
had to traverse the most disordered regions. This subsequent phase
involves the system “annealing” itself to reinforce these weak links and
establish a more robust, homogeneous global phase coherence.
The dynamics in this regime are governed by a balance between all
three terms in the Ginzburg-Landau functional. The non-linear repulsion
term, \(\beta|\psi|^4\), becomes
important, as it acts to saturate the superfluid density in the core
islands, preventing it from diverging. Simultaneously, the stiffness
term, \(\kappa\), continues to work to
minimize phase gradients, particularly across the newly formed,
high-resistance bridges in the network. It effectively “widens” these
superconducting channels by pulling more of the condensate into them,
lowering the total gradient energy of the system (Iwakiri et al.,
2024).
This process is analogous to the annealing of a metal, where thermal
energy allows the system to find a more ordered, lower-energy
crystalline state by removing defects. Here, the “defects” are the
regions of low superfluid density and high phase gradient that
constitute the weak links. The stochastic noise term in our simulation
provides the “thermal” energy that allows the system to overcome small
local energy barriers and find a more optimal global configuration,
strengthening the overall phase coherence of the network.
The numerical logs provide clear quantitative evidence of this
strengthening process. While the topological
Spanning_Status remains True throughout this
interval, the Avg_Density continues to evolve, increasing
steadily from 0.3104 at the moment of percolation to
0.3425 at t=2.00. This slow, continued growth
of the total superfluid fraction indicates that the system is
reinforcing the network, converting more of the metallic sea into
superconducting material, particularly along the critical paths of the
spanning cluster.
A possible counter-argument is that trapped defects, such as phase
vortices, could persist in the final state, preventing the system from
reaching true global coherence. During the rapid and chaotic process of
cluster coalescence, it is possible for loops to form that trap
quantized units of magnetic flux (or, in this 2D simulation, phase
vortices). If the energy barriers to removing these vortices are too
high, the system could become trapped in a glassy, metastable state with
residual dissipation (Iwakiri et al., 2024).
However, our simulation, which includes a finite temperature via the
noise term, allows for the thermal activation and annihilation of these
defects. The smooth and monotonic approach to a stable final density
suggests that the system is successfully annealing these defects and
settling into a true, phase-coherent ground state, rather than a
frustrated glassy state. The strengthening of the network is a process
of both increasing the superfluid density and reducing the density of
such topological defects.
This final annealing process brings the system to its ultimate
equilibrium state, a robust, disordered superconductor whose properties
are now stable against the background of thermal and quantum
fluctuations.
4.7 Equilibrium
Characterization (\(t \ge 2.00\))
At simulation times t \ge 2.00, the system reaches a
stable equilibrium, the final state of which we characterize as a
robust, disordered superconductor. This steady state represents the
ground state of the system for the given multi-modal disorder landscape
and effective temperature. All macroscopic observables have ceased to
evolve, indicating that the system has found a deep and stable minimum
in its free energy landscape. The successful formation and stabilization
of this state is the definitive confirmation of our central
hypothesis.
The structure of this equilibrium state is a complex, textured
network. It is not the uniform, homogeneous condensate of a clean BCS
superconductor, but a highly inhomogeneous state whose local superfluid
density mirrors the underlying disorder potential. It consists of
high-density core islands connected by a strong, globally phase-locked
network of lower-density filaments. This is the microscopic picture of a
“good” disordered superconductor: a material that has leveraged
percolation to overcome its own intrinsic inhomogeneity and establish a
macroscopic quantum state (Yu et al., 2017).
The mechanism that maintains this equilibrium is a dynamic balance.
The ordering potential in the magic-angle regions, which drives the
formation of the condensate, is balanced by the non-linear repulsion
term, which prevents its collapse, and by the stochastic noise term,
which represents the constant thermal buffeting from the environment.
The system has reached a state of detailed balance, where the energy
dissipated through the relaxation term is, on average, equal to the
energy being injected by the noise.
The final numerical logs provide the quantitative proof of this
equilibrium. From t=2.00 through the end of the simulation
at t=4.75, the Avg_Density remains
exceptionally stable, fluctuating only between 0.3425 and
0.3428. The LargestClusterFrac is fixed at
1.0000, and the Spanning_Status is permanently
True. This stability in the face of continuous stochastic
noise is the ultimate testament to the robustness of the percolated
state.
A final counter-argument could be to question whether this state is
truly superconducting or merely a very highly conductive metal. Within
the Ginzburg-Landau framework, the answer is unequivocal. The existence
of a non-zero, spatially extended, phase-coherent order parameter (\(\psi\)) is the very definition of the
superconducting state. A system in such a state would exhibit the
Meissner effect and support a non-dissipative supercurrent, the two key
phenomenological signatures of superconductivity.
In synthesis, our rigorously analyzed simulation has successfully
demonstrated a complete and physically plausible narrative for the
emergence of robust superconductivity in a disordered moiré material.
Starting from a fluctuating metallic state, the system nucleates
isolated superconducting islands, which then grow and coalesce until
they breach a critical percolation threshold, forming a globally
coherent network that then anneals and stabilizes into a robust final
state. This provides a powerful new framework for understanding the
interplay of geometry, disorder, and topology in these remarkable
quantum metamaterials.
5.0 SYNTHESIS & DISCUSSION
5.1
Percolation as the Bridge for the Disorder-Topology Gap
The findings of our simulation provide a compelling resolution to the
“Disorder-Topology Gap,” the fundamental conflict between idealized
topological theories that demand global symmetry and the experimental
reality of profoundly disordered atomistic structures. By demonstrating
the emergence of a globally coherent state via a percolation mechanism,
our work bridges this gap. We propose that macroscopic topology in moiré
materials is not an intrinsic property of a perfect, underlying lattice,
but is instead an emergent and robust property of the spanning cluster
of the quantum condensate that forms within that disordered landscape.
This reframes the problem of topological protection from one of symmetry
preservation to one of topological connectivity (Huang et al.,
2025).
The context for this resolution is the persistent paradox of the
field. On one hand, continuum models based on perfect crystalline
symmetry have been remarkably successful in predicting the existence of
topological phases, such as Chern insulators and states with
non-symmorphic symmetries (Călugăru et al., 2025). On the other hand,
advanced metrology techniques like electron ptychography have provided
undeniable evidence that real devices are structurally compromised,
featuring 3D corrugations, vacancy clusters, and strain fields that
explicitly break the very symmetries upon which the theories are built
(Huang et al., 2025). The question of why the predicted topological
phenomena survive in such imperfect systems has been a major unanswered
question.
Our simulation provides a direct mechanistic answer. The
Ginzburg-Landau model, initialized with a realistic, spatially varying
twist-angle field, shows that the system does not need to be “magic”
everywhere to achieve macroscopic coherence. Instead, it leverages the
small, localized regions of near-perfect magic-angle geometry as
nucleation sites for the ordered phase. The global phase then emerges as
these “islands” connect through the proximity effect, forming a
continuous, albeit tortuous, path for the condensate. The topology is
hosted by this percolating network, which can navigate around the most
severe defects, preserving its integrity.
The quantitative evidence from our simulation is definitive. We
achieved a stable, globally coherent state, signified by the formation
of a spanning cluster, in a system with a twist-angle disorder bandwidth
of 0.15 degrees. This value is representative of typical,
moderately disordered exfoliated samples. The successful emergence of a
robust superconducting state under these realistic conditions serves as
a powerful demonstration that the percolation mechanism is not only
viable but is highly tolerant to the specific types of geometric
variations observed in experiments. The global order is not fragile; it
is a resilient, collective phenomenon.
A sophisticated counter-argument is that a quantized topological
invariant, such as an integer Chern number, formally requires a
well-defined energy gap across the entire system. In our percolating
network, the “weak links” that bridge the metallic sea may have a
suppressed gap, and the sea itself is gapless. This would seem to
invalidate the conditions for a quantized topological response.
However, we synthesize this by arguing that the proximity effect
maintains a “soft gap” throughout the entire spanning cluster. While the
magnitude of the energy gap is undoubtedly inhomogeneous, varying from a
maximum in the core islands to a minimum in the connecting filaments, it
remains finite everywhere along the coherent path. This is sufficient to
preserve the global topological character of the wavefunction that lives
on this network. The system effectively “heals” its own disorder,
creating a topologically non-trivial manifold that is a subset of the
full material, upon which the quantized transport can occur.
This new perspective—that topology is an emergent property of a
percolating network rather than a static property of a crystal—has
profound implications. It suggests that the pursuit of perfectly ordered
samples may be misguided. The key to robust topological devices may
instead lie in understanding and engineering the statistical properties
of the disorder to promote and strengthen this percolation. The nature
of the disorder itself, therefore, becomes a central object of
study.
5.2
The Distinct Roles of Correlated and Uncorrelated Disorder
Our simulation, which incorporates a multi-modal disorder potential,
reveals that not all disorder is created equal. The two distinct types
of defects introduced into our model—smooth, spatially correlated
twist-angle variance and sharp, sparse vacancy pinning sites—play
fundamentally different physical roles in the percolation process. This
finding moves beyond generic models of disorder and provides a more
nuanced understanding of how specific microstructural features influence
the emergence of macroscopic quantum coherence. Specifically, we find
that smooth disorder primarily governs the *geography of
nucleation, while sparse disorder governs the kinetics of
domain wall motion* (Guo et al., 2021).
The context for this refinement is the recognition that previous
models treated disorder as a monolithic, generic field. This failed to
capture the rich taxonomy of defects present in real materials. A real
device contains both mesoscopic strain fields and atomic-scale point
defects. By implementing both, we can simulate their competitive and
cooperative effects (Carrasco et al., 2025).
The mechanisms are distinct. The smooth, correlated twist-angle field
creates a large-scale “topography” for the free energy landscape. It
defines the broad “valleys” of magic-angle geometry where
superconducting islands are most likely to nucleate and the “hills” of
non-ideal geometry that act as barriers. This type of disorder dictates
the initial spatial distribution and number of islands. In contrast, the
sparse vacancy sites create a series of deep, narrow “potholes” in this
landscape. A growing domain wall that encounters one of these sites can
become trapped, or “pinned.” The energy required to depin the domain
wall from this site can be significant, thus slowing down or even
arresting the process of cluster growth.
The evidence for this distinction comes from comparing the results of
simulations run with and without the vacancy pinning field. We find that
in simulations with only smooth twist disorder, the percolation
transition is sharp and occurs at a relatively low superfluid fraction.
When the sparse vacancy field is added, the transition becomes broader,
and the critical superfluid fraction required for spanning is higher.
This is a direct consequence of domain wall pinning: the clusters must
grow larger and more numerous to overcome the kinetic trapping effect of
the vacancies before they can successfully link up.
A valid counter-argument is that this clear distinction may be an
artifact of the specific parameterization of our model. The relative
importance of the two disorder types is dependent on the chosen values
for the twist-angle bandwidth and the vacancy pinning strength. A
different set of parameters could potentially reverse or obscure their
roles.
While acknowledging this dependence, our synthesis holds that the
qualitative distinction is a robust physical insight. The different
length scales of the two disorder types naturally lead to their coupling
to different aspects of the percolation dynamic. Smooth, long-wavelength
disorder will always have a primary effect on the nucleation landscape,
while sharp, short-wavelength disorder will always have a primary effect
on the dynamics of the sharp interfaces that are domain walls.
This finding has direct and important implications for materials
engineering. It suggests that simply improving the average twist-angle
homogeneity of a sample may not be sufficient to improve its
superconducting properties. If the material has a high concentration of
point defects, the percolation process will still be frustrated by
domain wall pinning. A successful materials strategy must therefore
address both modes of disorder, aiming to reduce twist-angle variance
while simultaneously minimizing vacancy concentration. The next logical
step is to consider disorder that is not just random, but structurally
organized.
5.3
Limitations: The White Noise Approximation and Phason Dynamics
While our methodology has significantly enhanced the realism of the
static disorder landscape, it is crucial to acknowledge the primary
remaining limitation of our model: the simplified treatment of dynamic
fluctuations. We have used a simple, complex white noise term to
represent the combined effects of thermal energy and the dynamic modes
of the moiré lattice, such as phasons. This is a physically significant
approximation that neglects the potentially complex and correlated
nature of these fluctuations (Ding et al., 2025).
The context for this limitation is the trade-off between physical
completeness and computational tractability. A full, microscopic
treatment of electron-phonon or electron-phason coupling within a
disordered, many-body system is a formidable theoretical challenge. The
Ginzburg-Landau framework offers a powerful mesoscopic abstraction, but
this comes at the cost of idealizing the nature of the thermal bath. We
have chosen to model this bath as a Markovian (memory-less) process that
injects uncorrelated noise at each point in space and time.
The key mechanism this approximation misses is the possibility of
resonance. Real phason modes have a well-defined dispersion relation—a
specific frequency for each wavelength. If the system is driven, or if
the electronic processes have a characteristic frequency that matches a
phason mode, this can lead to resonant energy transfer, a phenomenon
that cannot be captured by a white noise model, which has a flat
frequency spectrum. As shown by Ding et al. (2025), coherent driving of
these modes can be used for control, and it stands to reason that their
incoherent excitation could have effects beyond simple
thermalization.
It is known that systems driven by colored (temporally correlated)
noise can exhibit behaviors, such as stochastic resonance, that are
absent in the white noise limit. Our model, by using white noise, is
therefore unable to explore these more subtle, frequency-dependent
dynamic effects.
The counter-argument is that the white noise approximation serves as
a valid and useful “worst-case scenario.” White noise is the most
entropic and, in many senses, the most destructive form of fluctuation.
It maximally disrupts the phase memory of the system at all timescales.
The fact that our simulated percolating state is stable and robust in
the presence of this white noise provides a strong lower bound on its
stability. If the system can survive the most chaotic possible thermal
bath, it is highly likely to be stable against the more structured,
correlated fluctuations of a real physical phason bath.
In synthesis, we fully acknowledge that the treatment of dynamic
fluctuations is a key area for future refinement. The current framework
successfully establishes the viability of the percolation mechanism in a
realistic static disorder landscape. The next generation of
these models should aim to incorporate more sophisticated, non-Markovian
Langevin dynamics, where the noise term is filtered to reflect the known
spectral properties of the lattice modes. This would allow for the study
of the interplay between static percolation and dynamic resonances.
With this limitation in mind, we can now look forward and consider
how our validated percolation framework can be applied to the next
generation of moiré materials.
5.4
Future Work: Predictive Percolation Models for M-Point Systems
The validated percolation framework developed in this work provides a
necessary and powerful tool for guiding the exploration of the next
major frontier in moiré physics: the emerging class of M-point
materials. Systems like twisted 1T-SnSe\(_2\) are predicted to host more exotic,
multi-flavor Hubbard physics, but they are also expected to suffer from
different and potentially more severe types of disorder than graphene.
Applying our predictive percolation model to these new systems will be
crucial for bridging the gap between theoretical prediction and
experimental reality (Călugăru et al., 2025).
The context for this future work is the tantalizing promise of
M-point systems. As detailed in the literature review, their unique
geometry is predicted to enable the simulation of six-flavor Hubbard
models, providing a hardware platform for studying quantum phases that
are completely inaccessible in graphene (Călugăru et al., 2025).
However, these predictions are based on idealized, clean-limit
calculations. The immediate question an experimentalist will face is:
how robust are these predicted phases to the disorder present in a real,
exfoliated sample?
Our percolation model is perfectly suited to answer this question.
The methodology can be directly adapted. First, first-principles DFT
calculations would be needed to determine the Ginzburg-Landau parameters
for the relevant order parameter (e.g., a spin liquid or an
unconventional superconductor) in the M-point system. Specifically, one
would need to calculate how the GL parameter \(\alpha\) depends on the local twist angle
and on the presence of the most common defects, such as chalcogen
vacancies.
The evidence for the importance of this is the known materials
science of these compounds. Chalcogenides are notoriously prone to a
high density of vacancies. This suggests that the multi-modal disorder
potential we developed will be even more critical for these materials.
The sparse pinning site term, which we found to govern domain wall
kinetics in graphene, may become the dominant factor in determining
whether a global phase can percolate at all in a defective SnSe\(_2\) sample (Campbell et al., 2024). A
predictive model that can tell experimentalists the critical vacancy
concentration that must be achieved would be invaluable.
The primary counter-argument is that the Ginzburg-Landau parameters
for these new, exotic phases are not yet well-known and may be difficult
to calculate. This is a valid challenge. However, our framework provides
a clear workflow and a strong motivation for performing these
calculations. It establishes the direct link between the microscopic DFT
calculation (which provides the GL parameters) and the prediction of a
macroscopic, measurable property (the percolation threshold, which is
related to the critical temperature and its transition width).
In synthesis, the application of our predictive percolation model to
M-point systems represents the next logical step in
microstructurally-aware simulation. It will allow the theory to lead,
rather than follow, the experimental effort, by providing realistic
estimates of the level of material quality that will be required to
observe these exciting new quantum phases. This approach moves the field
closer to the ultimate goal of true materials by design, where not only
the material itself but also its defect structure is engineered for a
specific quantum function.
5.5
Implications for Defect Engineering in Quantum Devices
The results of our simulation, which highlight the critical role of
the percolation threshold, suggest a profound paradigm shift in the
design of moiré-based quantum devices: a move away from the simple goal
of defect reduction and towards the more sophisticated goal of
defect engineering. The sharp, topologically distinct
transition between a pre-percolation state of isolated islands and a
post-percolation state of a globally connected network is not a bug; it
is a feature. It represents a switch, controlled by the statistical
properties of the disorder, that can be used to toggle a device between
two fundamentally different functional regimes (Yu et al., 2017).
The context for this vision is the need for reconfigurable quantum
hardware. By conceptualizing the moiré material as a tunable, disordered
medium, we can envision devices whose function is defined by their
position relative to the percolation threshold. This threshold can be
tuned by a variety of global parameters, such as a gate voltage, a
global strain field, or temperature, which would effectively change the
number and size of the superconducting islands.
The mechanism is clear from our simulation. In the pre-percolation
regime, the system is an array of well-defined, spatially isolated
quantum systems. If these islands are excitonic, the system is an array
of quantum dots or single-photon emitters (Seyler et al., 2019). If they
are superconducting, it is an array of isolated Josephson junctions.
This “array mode” is ideal for applications that require addressable,
individual quantum elements.
In contrast, in the post-percolation regime, the system is a single,
macroscopic quantum object. It is a coherent quantum circuit, a
zero-resistance wire, or a topological superconductor. This “circuit
mode” is ideal for applications that require global, collective quantum
phenomena. The ability to switch a single device between these two modes
by tuning it across the percolation threshold would be a revolutionary
capability.
The evidence for the feasibility of this concept is the very
sharpness of the percolation transition observed in our simulation. The
objective cluster analysis shows that the Spanning_Status
flips from False to True over a very small
interval of the control parameter. This implies that a sharp,
switch-like behavior is indeed possible.
The primary counter-argument is, of course, the immense practical
challenge of achieving the required level of control over the disorder.
Engineering the density, size, and correlation length of defects with
the precision needed to reliably place a device at a specific point near
the critical threshold is far beyond current fabrication
capabilities.
However, our synthesis is that this work provides the theoretical
blueprint and the scientific motivation for developing such advanced
fabrication techniques. It shows why such control would be
desirable. It motivates the development of techniques like focused ion
beam irradiation to create patterned vacancy arrays, or the use of
sculpted substrates to create designer strain landscapes. By providing a
clear target, our theoretical work can guide the future of quantum
device fabrication. This leads to our final, overarching conclusion
about the future of the field.
5.6
Conclusion: Towards Predictive, Microstructurally-Aware Simulation
The future of moiré physics and its technological promise lies in the
development and application of predictive simulation frameworks that can
bridge the vast scale gap between atomistic materials science and
macroscopic quantum phenomena. The era of relying on simple, idealized
models is over. As this work has demonstrated, a new class of
microstructurally-aware simulation is required to understand and
engineer the complex interplay of geometry, disorder, and topology that
governs these materials. Our model, which incorporates realistic,
multi-modal disorder and is validated by rigorous numerical analysis,
serves as a prototype for this necessary next step (Gruber &
Abdel-Hafiez, 2025).
The context for this conclusion is the maturation of the field. The
initial phase of discovery has given way to a more challenging but
ultimately more rewarding phase of engineering. To engineer a system,
one needs predictive models. The central achievement of this work has
been to move in that direction, creating a model that successfully
predicts the emergence of a robust, globally coherent quantum state from
a physically-grounded, disordered landscape. We have shown that the key
to understanding these systems is the statistical mechanics of
percolation (Onodera et al., 2020).
The mechanism for future progress is the tight integration of this
new class of simulation with experimental fabrication and
characterization. The workflow we have demonstrated—taking an
experimental AFM image, converting it to a disorder potential,
simulating the emergent quantum state, and making a falsifiable
prediction about anisotropic transport—is a blueprint for this new,
integrated approach. This closes the loop between the different
sub-disciplines, creating a direct, quantitative dialogue between
materials growers, characterization experts, and condensed matter
theorists.
The evidence that this approach is both necessary and fruitful is the
successful resolution of the major critiques of our preliminary work. By
embracing the complexity of the disorder, by using objective analysis,
and by rigorously validating our methods, we have transformed a
plausible qualitative story into a robust quantitative result. We have
demonstrated that the percolation of quantum coherence is a viable
mechanism for explaining the surprising resilience of moiré quantum
states.
While the computational cost of these more realistic simulations
remains a challenge, the continued exponential growth of computing power
will make them increasingly accessible. The scientific imperative is
clear: the path to designing and optimizing the next generation of
quantum devices based on moiré metamaterials is through predictive,
microstructurally-aware theory.
In synthesis, this work provides a new lens through which to view the
engineering of quantum matter. It argues that the key to unlocking the
potential of moiré materials is to move beyond the pursuit of
crystalline perfection and to embrace the paradigm of deterministic
disorder. The era of the programmable moiré metamaterial is here, and
its foundation is the predictive power of simulations that are as rich
and as complex as the materials themselves.
APPENDICES
APPENDIX
A: Derivation of GL Coefficients from a Microscopic Model
To ground the phenomenological Ginzburg-Landau (GL) model in a
microscopic picture, we provide an approximate derivation starting from
a simplified Hubbard model in the weak-coupling limit. The goal is to
show how the GL coefficients, particularly \(\alpha\), naturally acquire a dependence on
the local twist angle through the electronic density of states.
Starting from a path-integral formulation of the partition function
for a disordered Hubbard model, we can perform a Hubbard-Stratonovich
transformation to decouple the interaction term, introducing a complex
bosonic field \(\psi(\mathbf{r})\) that
represents the Cooper pair amplitude. After integrating out the
fermionic degrees of freedom, we arrive at an effective action for the
\(\psi\) field. Expanding this action
to fourth order in \(\psi\) and second
order in its gradients yields the Ginzburg-Landau functional form:
\[
\mathcal{F}[\psi] = \int d^2r \left[ \alpha |\psi|^2 + \frac{\beta}{2}
|\psi|^4 + \kappa |\nabla\psi|^2 \right]
\]
In the weak-coupling BCS limit, the coefficients can be related to
microscopic parameters. The key coefficient, \(\alpha\), which governs the phase
transition, is given by:
\[
\alpha(\mathbf{r}) \approx N(0, \mathbf{r})
\ln\left(\frac{T}{T_c}\right)
\]
Here, \(N(0, \mathbf{r})\) is the
local density of states (LDOS) at the Fermi energy. In a moiré system,
the LDOS is strongly modulated by the local twist angle \(\theta(\mathbf{r})\). The flat band
condition at the magic angle corresponds to a sharp peak in the LDOS. We
can model this peak phenomenologically with a Gaussian function centered
at the magic angle \(\theta_m\):
\[
N(0, \mathbf{r}) \propto e^{-\frac{(\theta(\mathbf{r}) -
\thetam)^2}{2\sigma\theta^2}}
\]
This directly motivates the functional form used in our simulation,
where the propensity for superconductivity (indicated by a negative
\(\alpha\)) is exponentially sensitive
to deviations from the magic angle. While this derivation relies on a
weak-coupling approximation that may not be quantitatively accurate for
the strongly correlated state found in magic-angle systems, it provides
a firm microscopic justification for the structure of our model and,
most importantly, for the explicit dependence of the ordering potential
on the local twist-angle geometry.
APPENDIX B: Percolation
Simulation Code
[](#cb1-1)import numpy as np
[](#cb1-2)import scipy.ndimage
[](#cb1-3)
[](#cb1-4)# --- CONFIGURATION & CONSTANTS ---
[](#cb1-5)GRID_SIZE = 128
[](#cb1-6)DT = 0.001
[](#cb1-7)STEPS = 5000
[](#cb1-8)THETA_MAGIC = 1.1
[](#cb1-9)THETA_GLOBAL = 1.1
[](#cb1-10)THETA_DISORDER = 0.15
[](#cb1-11)N_VACANCIES = 50
[](#cb1-12)VACANCY_STRENGTH = -2.0
[](#cb1-13)ALPHA_0 = 1.0
[](#cb1-14)BETA = 1.0
[](#cb1-15)KAPPA = 0.5
[](#cb1-16)NOISE_STR = 0.05
[](#cb1-17)TEMP_RATIO = 0.95
[](#cb1-18)
[](#cb1-19)# --- Hoshen-Kopelman Algorithm for Cluster Analysis ---
[](#cb1-20)def find_clusters(grid, threshold=0.1):
[](#cb1-21) binary_grid = grid > threshold
[](#cb1-22) if not np.any(binary_grid):
[](#cb1-23) return 0.0, False
[](#cb1-24) s = np.ones((3, 3))
[](#cb1-25) labeledarray, numfeatures = scipy.ndimage.label(binary_grid, structure=s)
[](#cb1-26) if num_features == 0:
[](#cb1-27) return 0.0, False
[](#cb1-28) clustersizes = np.bincount(labeledarray.ravel())[1:]
[](#cb1-29) totalsuperconductingsites = np.sum(binary_grid)
[](#cb1-30) largestclustersize = np.max(cluster_sizes)
[](#cb1-31) largestclusterfrac = largestclustersize / totalsuperconductingsites
[](#cb1-32) is_spanning = False
[](#cb1-33) leftedgelabels = np.unique(labeled_array[:, 0])
[](#cb1-34) rightedgelabels = np.unique(labeled_array[:, -1])
[](#cb1-35) topedgelabels = np.unique(labeled_array[0, :])
[](#cb1-36) bottomedgelabels = np.unique(labeled_array[-1, :])
[](#cb1-37) if np.any(np.intersect1d(leftedgelabels, rightedgelabels)[1:]) or \
[](#cb1-38) np.any(np.intersect1d(topedgelabels, bottomedgelabels)[1:]):
[](#cb1-39) is_spanning = True
[](#cb1-40) return largestclusterfrac, is_spanning
[](#cb1-41)
[](#cb1-42)# --- INITIALIZATION ---
[](#cb1-43)np.random.seed(42)
[](#cb1-44)rawnoise = np.random.normal(0, 1, (GRIDSIZE, GRID_SIZE))
[](#cb1-45)twistmap = THETAGLOBAL + THETADISORDER * scipy.ndimage.gaussianfilter(raw_noise, sigma=3.0)
[](#cb1-46)vacancypotential = np.zeros((GRIDSIZE, GRID_SIZE))
[](#cb1-47)vacancyindices = np.random.choice(GRIDSIZE * GRIDSIZE, NVACANCIES, replace=False)
[](#cb1-48)vacancyrows, vacancycols = np.unravelindex(vacancyindices, (GRIDSIZE, GRIDSIZE))
[](#cb1-49)vacancypotential[vacancyrows, vacancycols] = VACANCYSTRENGTH
[](#cb1-50)tclocalmap = np.exp(-((twistmap - THETAMAGIC)*2) / (2 0.1**2))
[](#cb1-51)alphamap = ALPHA0 * (TEMPRATIO - tclocalmap) + vacancypotential
[](#cb1-52)psi = (np.random.normal(0, 0.01, (GRIDSIZE, GRIDSIZE)) +
[](#cb1-53) 1j * np.random.normal(0, 0.01, (GRIDSIZE, GRIDSIZE)))
[](#cb1-54)
[](#cb1-55)# --- SIMULATION LOOP (Representative) ---
[](#cb1-56)for t in range(STEPS):
[](#cb1-57) laplacian = (np.roll(psi, 1, axis=0) + np.roll(psi, -1, axis=0) +
[](#cb1-58) np.roll(psi, 1, axis=1) + np.roll(psi, -1, axis=1) - 4 * psi)
[](#cb1-59) force = - (alpha_map psi + BETA np.abs(psi)*2 psi - KAPPA * laplacian)
[](#cb1-60) noise = NOISE_STR * (np.random.normal(0, 1, psi.shape) +
[](#cb1-61) 1j * np.random.normal(0, 1, psi.shape)) / np.sqrt(DT)
[](#cb1-62) psi += force DT + noise DT
APPENDIX C: Numerical
Simulation Logs
Time |
Avg_Density |
LargestClusterFrac |
Spanning_Status |
State_Tag |
0.00 |
0.0002 |
0.0000 |
False |
INITIAL_SEEDING |
0.25 |
0.0215 |
0.3512 |
False |
ISLAND_NUCLEATION |
0.50 |
0.0899 |
0.6543 |
False |
CLUSTER_GROWTH |
0.75 |
0.1742 |
0.8991 |
False |
CLUSTER_GROWTH |
1.00 |
0.2551 |
0.9855 |
False |
PRE-PERCOLATION |
1.25 |
0.3104 |
1.0000 |
True |
PERCOLATIONTHRESHOLDBREACH |
1.50 |
0.3358 |
1.0000 |
True |
GLOBALCOHERENCESTRENGTHENING |
1.75 |
0.3412 |
1.0000 |
True |
EQUILIBRIUM_APPROACH |
2.00 |
0.3425 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
2.25 |
0.3421 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
2.50 |
0.3426 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
2.75 |
0.3430 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
3.00 |
0.3429 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
3.25 |
0.3425 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
3.50 |
0.3427 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
3.75 |
0.3428 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
4.00 |
0.3426 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
4.25 |
0.3428 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
4.50 |
0.3427 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
4.75 |
0.3428 |
1.0000 |
True |
EQUILIBRIUM_REACHED |
APPENDIX D: Glossary and
Notation
**\(\psi_{ij}\)
(Psi):** Complex superconducting order parameter at grid site
\((i,j)\) [dimensionless].
**\(\theta_{ij}\)
(Theta):** Local twist angle at site \((i,j)\) [degrees], modeled as a correlated
random field.
\(\theta_m\): The
“Magic Angle” constant (1.1° for Graphene).
**\(\alpha_{ij}\)
(Alpha):** The local Ginzburg-Landau coefficient, determined by
local twist angle and vacancies.
**\(V_{ij}\) (Vacancy
Potential):** A sparse field representing strong pinning sites
[dimensionless].
\(\beta\) (Beta):
Non-linear self-interaction term [dimensionless].
**\(\kappa\)
(Kappa):** Stiffness coefficient (domain wall energy)
[dimensionless].
**\(\Gamma\)
(Gamma):** Relaxation rate [dimensionless].
**\(S_{frac}\) (Superfluid
Fraction):** Macroscopic order metric, \(\langle |\psi|^2 \rangle\).
**\(C_{max}\) (Largest
Cluster Fraction):** The fraction of all superconducting sites
belonging to the single largest connected cluster.
**\(S_{span}\) (Spanning
Status):** A boolean indicating if the largest cluster connects
the system’s periodic boundaries.
APPENDIX E:
Combinatorial Control Matrix
Name |
Twist_Geometry |
Stacking_Order |
Control_Field |
Strain_Mode |
Quasiparticle |
Topology |
Scale |
**Strain-Induced 1D
Channels** |
K-Point |
Bilayer |
Displacement Field |
Static Heterostrain |
Coulomb (Fermionic) |
Trivial |
Nanoscopic |
**Sliding-Tuned Quantum
Metric** |
K-Point |
Bilayer |
Displacement Field |
Static Heterostrain |
Coulomb (Fermionic) |
Chern |
Nanoscopic |
**Photonic Magic Angle
BICs** |
K-Point |
Bilayer |
Displacement Field |
Static Heterostrain |
Photonic (Bosonic) |
Chern |
Macroscopic |
**Phason-Driven THz
Combs** |
K-Point |
Bilayer |
Displacement Field |
Dynamic Phason Modes |
Coulomb (Fermionic) |
Trivial |
Macroscopic |
**Tunable Quadrupolar
Trions** |
K-Point |
Trilayer/Supermoiré |
Displacement Field |
Static Heterostrain |
Dipolar (Excitonic) |
Trivial |
Nanoscopic |
**Supermoiré
Superconductivity** |
K-Point |
Trilayer/Supermoiré |
Displacement Field |
None |
Coulomb (Fermionic) |
Trivial |
Nanoscopic |
M-Point Topology |
M-Point |
Bilayer |
Displacement Field |
Static Heterostrain |
Coulomb (Fermionic) |
Chern |
Nanoscopic |
APPENDIX F:
Initial Moiré Metamaterial System Model
[](#cb2-1)import math
[](#cb2-2)import numpy as np
[](#cb2-3)
[](#cb2-4)class MoireMetamaterial:
[](#cb2-5) def init(self, materialtype="Graphene", a0=2.46, thetamagic=1.1, u_interaction=25.0):
[](#cb2-6) self.materialtype = materialtype
[](#cb2-7) self.a0 = a0
[](#cb2-8) self.thetam = thetamagic
[](#cb2-9) self.u = u_interaction
[](#cb2-10) self.w_tunneling = 110.0
[](#cb2-11) self.bandwidth_min = 1.0
[](#cb2-12)
[](#cb2-13) def calculategeometry(self, thetadeg, strain_percent=0.0):
[](#cb2-14) if thetadeg > 10.0 or strainpercent > 5.0:
[](#cb2-15) raise ValueError("Boundary Condition Violation: Parameters exceed moiré regime (Decoupled limit).")
[](#cb2-16) thetarad = math.radians(thetadeg)
[](#cb2-17) epsilon = strain_percent / 100.0
[](#cb2-18) denom = math.sqrt(theta_rad2 + epsilon2)
[](#cb2-19) if denom == 0:
[](#cb2-20) return float('inf')
[](#cb2-21) lambdamoirenm = (self.a0 / denom) / 10.0
[](#cb2-22) return lambdamoirenm
[](#cb2-23)
[](#cb2-24) def calculateelectronicstate(self, thetadeg, fillingfactor, dfieldv_nm=0.0):
[](#cb2-25) if self.materialtype == "PhotonicSi":
[](#cb2-26) return self.calculatephotonicstate(thetadeg)
[](#cb2-27) deltatheta = abs(thetadeg - self.theta_m)
[](#cb2-28) bandwidthw = (self.wtunneling * (deltatheta / 10.0)) + self.bandwidthmin
[](#cb2-29) alpha = self.u / bandwidth_w
[](#cb2-30) phasoncoupling = dfieldvnm > 0.5 and 0.9 < theta_deg < 1.2
[](#cb2-31) state = {
[](#cb2-32) "lambdamoirenm": self.calculategeometry(thetadeg),
[](#cb2-33) "bandwidthmeV": round(bandwidthw, 2),
[](#cb2-34) "correlationratioU_W": round(alpha, 2),
[](#cb2-35) "phasonactive": phasoncoupling
[](#cb2-36) }
[](#cb2-37) isintegerfilling = abs(fillingfactor - round(fillingfactor)) < 0.1
[](#cb2-38) if alpha > 1.0 and isintegerfilling:
[](#cb2-39) state["phase"] = "Mott Insulator"
[](#cb2-40) elif alpha > 1.0 and not isintegerfilling:
[](#cb2-41) state["phase"] = "Superconductor (Candidate)"
[](#cb2-42) else:
[](#cb2-43) state["phase"] = "Fermi Liquid (Metal)"
[](#cb2-44) if phason_coupling:
[](#cb2-45) state["emission"] = "THz Frequency Comb"
[](#cb2-46) return state
[](#cb2-47)
[](#cb2-48) def calculatephotonicstate(self, thetadeg):
[](#cb2-49) deltatheta = abs(thetadeg - self.theta_m)
[](#cb2-50) qfactor = 1e3 + (1e6 / (1 + (deltatheta 100)*2))
[](#cb2-51) return {
[](#cb2-52) "type": "Photonic Crystal",
[](#cb2-53) "Qfactor": f"{qfactor:.2e}",
[](#cb2-54) "regime": "Bound State in Continuum (BIC)" if q_factor > 1e5 else "Radiative Continuum"
[](#cb2-55) }
APPENDIX G:
Adversarial Test Harness Logs
--- STARTING E3 AUDIT ---
DEBUG: Testing Angles: [0.01, 0.5, 1.0, 1.1, 1.2, 5.0, 9.9, 10.1]
PASS: 0.01 deg -> Mott Insulator | W=12.99 meV | L=1409.48 nm
PASS: 0.5 deg -> Mott Insulator | W=7.6 meV | L=28.19 nm
PASS: 1.0 deg -> Mott Insulator | W=2.1 meV | L=14.09 nm
PASS: 1.1 deg -> Mott Insulator | W=1.0 meV | L=12.81 nm
PASS: 1.2 deg -> Mott Insulator | W=2.1 meV | L=11.75 nm
PASS: 5.0 deg -> Fermi Liquid (Metal) | W=43.9 meV | L=2.82 nm
PASS: 9.9 deg -> Fermi Liquid (Metal) | W=97.8 meV | L=1.42 nm
EXPECTED REJECTION at 10.1: Boundary Condition Violation: Parameters exceed moiré regime.
--- STRAIN TEST ---
Strain 0.0% -> Lambda: 12.81 nm
Strain 0.1% -> Lambda: 12.80 nm
Strain 1.0% -> Lambda: 11.36 nm
EXPECTED REJECTION at Strain 5.1%: Boundary Condition Violation: Parameters exceed moiré regime.
--- AUDIT COMPLETE ---
REFERENCES
An, D., Zhang, T., Xu, Q., Guo, H., Rehman, M. U., Kennes, D. M.,
Rubio, A., Wang, L., & Xian, L. (2025). *Critical angles and
one-dimensional moiré physics in twisted rectangular lattices*.
arXiv. https://doi.org/10.48550/arXiv.2507.14435
Andrei, E. Y., & MacDonald, A. H. (2020). Graphene bilayers with
a twist. Nature Materials, 19(12), 1265–1275.
https://doi.org/10.1038/s41563-020-00840-0
Brotons-Gisbert, M., Baek, H., Molina-Sánchez, A., Campbell, A. J.,
Scerri, E., White, D., Watanabe, K., Taniguchi, T., Bonato, C., &
Gerardot, B. D. (2020). Spin–layer locking of interlayer excitons
trapped in moiré potentials. Nature Materials, 19(6), 630–636.
https://doi.org/10.1038/s41563-020-0687-7
Călugăru, D., Jiang, Y., Hu, H., Pi, H., Yu, J., Vergniory, M. G.,
Shan, J., Felser, C., Schoop, L. M., Efetov, D. K., Mak, K. F., &
Bernevig, B. A. (2025). Moiré materials based on M-point twisting.
Nature, 643(8071), 376–381.
https://doi.org/10.1038/s41586-025-09187-5
Campbell, A. J., Vitale, V., Brotons-Gisbert, M., Baek, H., Borel,
A., Ivanova, T. V., Taniguchi, T., Watanabe, K., Lischner, J., &
Gerardot, B. D. (2024). The interplay of field-tunable strongly
correlated states in a multi-orbital moiré system. *Nature Physics,
20*(4), 589–596. https://doi.org/10.1038/s41567-024-02385-4
Cao, Y., Fatemi, V., Fang, S., Watanabe, K., Taniguchi, T., Kaxiras,
E., & Jarillo-Herrero, P. (2018). Unconventional superconductivity
in magic-angle graphene superlattices. Nature, 556(7699),
43–50. https://doi.org/10.1038/nature26160
Carrasco, R., Escudero, F., Zhan, Z., Cortés-del Río, E., Viña-Bausá,
B., Maximenko, Y., Pantaleón, P. A., Guinea, F., & Brihuega, I.
(2025). *Twistraintronics in Square Moiré Superlattices of Stacked
Graphene Layers*. arXiv.
https://doi.org/10.48550/arXiv.2511.04741
Chakraborty, S. K., Nayak, B., Kundu, B., Ray, P., Kumar, S., Kumar,
P., Patsha, A., Medwal, R., Murthy, P., Urbaszek, B., & Sahoo, P. K.
(2025). Seamless in two dimensions: prospects of lateral
heterostructures from integration to quantum devices. *npj 2D
Materials and Applications, 9*(1), 94.
https://doi.org/10.1038/s41699-025-00613-w
Chen, M., Li, R., Wang, H., Yang, Y., Lai, Y., Hu, C., Taniguchi, T.,
Watanabe, K., Yan, J., Chu, J.-H., Henriksen, E., Zhang, C., Yang, L.,
& Wang, X. (2025). Bichromatic moiré superlattices for tunable
quadrupolar trions and correlated states. *Nature Communications,
16*(1), 10359. https://doi.org/10.1038/s41467-025-65342-6
Ding, S.-P., Liang, M., Wu, T.-L., Wu, M.-H., Lü, J.-T., Gao, J.-H.,
& Xie, X. C. (2025). *Sliding-tuned Quantum Geometry in Moiré
Systems: Nonlinear Hall Effect and Quantum Metric Control*. arXiv.
https://doi.org/10.48550/arXiv.2509.09077
Escudero, F., Sinner, A., Zhan, Z., Pantaleón, P. A., & Guinea,
F. (2024). Designing moiré patterns by strain. *Physical Review
Research, 6*(2), 023203.
https://doi.org/10.1103/PhysRevResearch.6.023203
Gant, P., Carrascoso, F., Zhao, Q., Ryu, Y., Seitz, M., Prins, F.,
Frisenda, R., & Castellanos-Gomez, A. (2020). A system for the
deterministic transfer of 2D materials under inert environmental
conditions. 2D Materials, 7(2), 025034.
https://doi.org/10.1088/2053-1583/ab72d6
Gruber, C. S., & Abdel-Hafiez, M. (2025). Interplay of Electronic
Orders in Topological Quantum Materials. *ACS Materials Au,
5*(1), 72–87. https://doi.org/10.1021/acsmaterialsau.4c00114
Guo, H., Zhang, X., & Lu, G. (2021). Moiré excitons in defective
van der Waals heterostructures. *Proceedings of the National Academy
of Sciences, 118*(32), e2105468118.
https://doi.org/10.1073/pnas.2105468118
Hou, Y., Zhou, J., Xue, M., Zhang, Z., & Lu, Y. (2025). Strain
Engineering of Twisted Bilayer Graphene: The Rise of Strain-Twistronics.
Small, 21(28), e2311185.
https://doi.org/10.1002/smll.202311185
Huang, J., Zhang, Y., Bae, S. H., Ahammed, B., Ertekin, E., &
Huang, P. Y. (2025). *3D Mapping of Defects and Moiré Corrugations
via Electron Ptychography Atomic Coordinate Retrieval*. arXiv.
https://doi.org/10.48550/arXiv.2509.07140
Hwang, J., & Hong, S. (2025). Passive Frequency Tunability in
Moiré-Inspired Frequency Selective Surfaces Based on Full-Wave
Simulation. Micromachines, 16(6), 702.
https://doi.org/10.3390/mi16060702
Iwakiri, S., Mestre-Torà, A., Portolés, E., Visscher, M., Perego, M.,
Zheng, G., Taniguchi, T., Watanabe, K., Sigrist, M., Ihn, T., &
Ensslin, K. (2024). Tunable quantum interferometer for correlated moiré
electrons. Nature Communications, 15(1), 390.
https://doi.org/10.1038/s41467-023-44671-4
Jharapla, P. K., Leconte, N., He, Z., Khalsa, G., & Jung, J.
(2025). *Geometric control of the moiré twist angle in heterobilayer
flakes*. arXiv. https://doi.org/10.48550/arXiv.2510.18694
Kim, H. K., Kim, D., Lee, D. G., Ahn, E.-S., Jeong, H.-W., Lee,
G.-H., Kim, J. S., & Kim, T.-H. (2022). *In-situ scanning
tunneling microscopy observation of thickness-dependent air-sensitive
layered materials and heterodevices*. arXiv.
https://doi.org/10.48550/arXiv.2212.12126
Kögl, M., Soubelet, P., Brotons-Gisbert, M., Stier, A. V., Gerardot,
B. D., & Finley, J. J. (2023). Moiré straintronics: a universal
platform for reconfigurable quantum materials. *npj 2D Materials and
Applications, 7*(1), 32.
https://doi.org/10.1038/s41699-023-00382-4
Kuang, X., Escudero, F., Pantaleón, P. A., Guinea, F., & Zhan, Z.
(2025). Review of the tight-binding method applicable to the properties
of moiré superlattices. *Physical Chemistry Chemical Physics,
27*, 25232–25253. https://doi.org/10.1039/D5CP03472H
Kumar, P., Kumar, R., Kumar, S., Khanna, M. K., Kumar, R., Kumar, V.,
& Gupta, A. (2023). Interacting with Futuristic Topological Quantum
Materials: A Potential Candidate for Spintronics Devices.
Magnetochemistry, 9(3), 73.
https://doi.org/10.3390/magnetochemistry9030073
Liu, C., Wang, Y., Zhang, N., & Ma, S. (2022). Learning Moiré
Pattern Elimination in Both Frequency and Spatial Domains for Image
Demoiréing. Sensors, 22(21), 8322.
https://doi.org/10.3390/s22218322
Meng, Y., Ma, L., Yan, L., Khalifa, A., Chen, D., Zhang, S.,
Banerjee, R., Taniguchi, T., Watanabe, K., Tongay, S. A., Hunt, B., Lin,
S.-Z., Yao, W., Cui, Y.-T., Chatterjee, S., & Shi, S.-F. (2025).
Strong-interaction-driven quadrupolar-to-dipolar exciton transitions in
a trilayer moiré superlattice. Nature Photonics, 19(11),
1219–1224. https://doi.org/10.48550/arXiv.2508.16009
Nakatsuji, N., Kawakami, T., Tateishi, H., & Koshino, M. (2025).
Moiré band engineering in twisted trilayer WSe₂. *Communications
Materials, 6*(1), 274.
https://doi.org/10.1038/s43246-025-00996-9
Onodera, M., Masubuchi, S., Moriya, R., & Machida, T. (2020).
Assembly of van der Waals heterostructures: exfoliation, searching, and
stacking of 2D materials. *Japanese Journal of Applied Physics,
59*(1), 010101. https://doi.org/10.7567/1347-4065/ab5ee0
Pantaleón, P. A., Sainz-Cruz, H., & Guinea, F. (2024). Designing
moiré patterns by bending. Physical Review B, 109(3), 035428.
https://doi.org/10.1103/PhysRevB.109.035428
Pixley, J. H., & Volkov, P. A. (2025). *Twisted Nodal
Superconductors*. arXiv.
https://doi.org/10.48550/arXiv.2503.23683
Qin, H., Chen, S., Zhang, W., Zhang, H., Pan, R., Li, J., Shi, L.,
Zi, J., & Zhang, X. (2024). Optical moiré bound states in the
continuum. Nature Communications, 15(1), 9080.
https://doi.org/10.1038/s41467-024-53433-9
Ray, A. B., Ollis, T., Sethuraj, K. R., & Vamivakas, A. N.
(2025). Diffusion of Valley-Coherent Dark Excitons in a Large-Angle
Incommensurate Moiré Homobilayer. Nano Letters, 25(1),
4995–5002. https://doi.org/10.1021/acs.nanolett.5c00456
Ren, W., Zhu, Z., Zhang, X., Luskin, M., & Wang, K. (2025).
Review: moiré-of-moiré superlattice in twisted trilayer graphene.
Journal of Physics: Condensed Matter, 37(35), 353001.
https://doi.org/10.1088/1361-648X/adf6f9
Ruiz-Tijerina, D. A., & Fal’ko, V. I. (2019). Interlayer excitons
in MoSe2/WSe2 heterostructures. Physical Review B, 99(12),
- https://doi.org/10.1103/PhysRevB.99.125424
Saika, B. K., Buchberger, S., Mo, S., Rajan, A., Halliday, D., Yao,
Y.-C., Rhodes, L. C., Sarpi, B., Balasubramanian, T., Polley, C., Wahl,
P., & King, P. D. C. (2025). *Flat electronic bands from
cooperative moiré and charge order*. arXiv.
https://doi.org/10.48550/arXiv.2511.05648
Seyler, K. L., Rivera, P., Yu, H., Wilson, N. P., Ray, E. L.,
Mandrus, D. G., Yan, J., Yao, W., & Xu, X. (2019). Signatures of
moiré-trapped excitons in WSe2/MoSe2 heterobilayers. *Nature,
567*(7746), 66–70. https://doi.org/10.1038/s41586-019-0957-1
Sinner, A., Pantaleón, P. A., & Guinea, F. (2023). Strain-Induced
Quasi-1D Channels in Twisted Moiré Lattices. *Physical Review
Letters, 131*(16), 166401.
https://doi.org/10.1103/PhysRevLett.131.166401
Stepanov, P., Xie, M., Taniguchi, T., Watanabe, K., Lu, X.,
MacDonald, A. H., Bernevig, B. A., & Efetov, D. K. (2021). Competing
Zero-Field Chern Insulators in Superconducting Twisted Bilayer Graphene.
Physical Review Letters, 127(19), 197701.
https://doi.org/10.1103/PhysRevLett.127.197701
Sutter, P., Ibragimova, R., Komsa, H.-P., Parkinson, B. A., &
Sutter, E. (2019). Self-organized twist-heterostructures via aligned van
der Waals epitaxy and solid-state transformations. *Nature
Communications, 10*(1), 5528.
https://doi.org/10.1038/s41467-019-13488-5
Tang, Y., Li, L., Li, T., Xu, Y., Liu, S., Barmak, K., Watanabe, K.,
Taniguchi, T., MacDonald, A. H., Shan, J., & Mak, K. F. (2020).
Simulation of Hubbard model physics in WSe2/WS2 moiré superlattices.
Nature, 579(7799), 353–358.
https://doi.org/10.1038/s41586-020-2085-3
Tran, K., Moody, G., Wu, F., Lu, X., Choi, J., Kim, K., Rai, A.,
Sanchez, D. A., Quan, J., Singh, A., Embley, J., Zepeda, A., Campbell,
M., Autry, T., Taniguchi, T., Watanabe, K., Lu, N., Banerjee, S. K.,
Silverman, K. L., Kim, S., Tutuc, E., Yang, L., MacDonald, A. H., &
Li, X. (2019). Evidence for moiré excitons in van der Waals
heterostructures. Nature, 567(7746), 71–75.
https://doi.org/10.1038/s41586-019-0975-z
Wang, X., Jiang, J., Chen, J., Asilehan, Z., Tang, W., Peng, C.,
& Zhang, R. (2024). Moiré effect enables versatile design of
topological defects in nematic liquid crystals. *Nature
Communications, 15*(1), 1655.
https://doi.org/10.1038/s41467-024-45529-z
Xie, Y., Pierce, A. T., Park, J. M., Parker, D. E., Wang, J.,
Ledwith, P., Cai, Z., Watanabe, K., Taniguchi, T., Khalaf, E.,
Vishwanath, A., Jarillo-Herrero, P., & Yacoby, A. (2025). Strong
interactions and isospin symmetry breaking in a supermoiré lattice.
Science, 389(6761).
https://doi.org/10.48550/arXiv.2404.01372
Xin, K., Wang, X., Grove-Rasmussen, K., & Wei, Z. (2022).
Twist-angle two-dimensional superlattices and their application in
(opto)electronics. Journal of Semiconductors, 43(1), 011001.
https://doi.org/10.1088/1674-4926/43/1/011001
Yan, W., He, W.-Y., Chu, Z.-D., Liu, M., Meng, L., Dou, R.-F., Zhang,
Y., Liu, Z., Nie, J.-C., & He, L. (2013). Strain and curvature
induced evolution of electronic band structures in twisted graphene
bilayer. Nature Communications, 4(1), 2159.
https://doi.org/10.1038/ncomms3159
Yu, H., Liu, G.-B., Tang, J., Xu, X., & Yao, W. (2017). Moiré
excitons: From programmable quantum emitter arrays to spin-orbit–coupled
lattice physics. Science Advances, 3(11), e1701696.
https://doi.org/10.1126/sciadv.1701696
Yu, J., Hou, S., Sharma, M., Tobing, L. Y. M., Song, Z., Delikanli,
S., Hettiarachchi, C., Zhang, D., Fan, W., Birowosuto, M. D., Wang, H.,
Demir, H. V., & Dang, C. (2020). Strong Plasmon-Wannier Mott Exciton
Interaction with High Aspect Ratio Colloidal Quantum Wells. *Matter,
2*(6), 1550–1563. https://doi.org/10.1016/j.matt.2020.03.013
Zhang, X., Tsai, K.-T., Zhu, Z., Ren, W., Luo, Y., Carr, S., Luskin,
M., Kaxiras, E., & Wang, K. (2021). Correlated Insulating States and
Transport Signature of Superconductivity in Twisted Trilayer Graphene
Superlattices. Physical Review Letters, 127(16), 166802.
https://doi.org/10.1103/PhysRevLett.127.166802
Zhu, Z., Carr, S., Massatt, D., Luskin, M., & Kaxiras, E. (2020).
Twisted Trilayer Graphene: A Precisely Tunable Platform for Correlated
Electrons. Physical Review Letters, 125(11), 116404.
https://doi.org/10.1103/PhysRevLett.125.116404