Duality between atomic conﬁgurations and Bloch states in twistronic materials

The relative orientation (twist) of successive layers of stacked two-dimensional (2D) materials creates variations in the interlayer atomic registry. The variations often form a superlattice, called a moiré pattern, which can alter electronic properties. In this work we introduce a classiﬁcation of the single-particle electronic structures that can occur in twisted stacks of 2D layers by characterizing them as “moiré molecules” or “moiré crystals.” The molecules generate localized electronic states and moiré ﬂat bands, while the crystals are sometimes unconventional and produce electronic banding in the conﬁguration basis. The underpinning of this classiﬁcation is the duality between interlayer conﬁguration and monolayer Bloch momentum in moiré Hamiltonians. We apply this understanding to diagrams of local electron density in untwisted geometries to produce intuitive and quantitative predictions of twistronic properties. We provide a conceptual introduction to this framework through a one-dimensional model, and then apply it to 2D twisted bilayers of the semimetal graphene, and of MoS 2 , a representative material of the transition metal dichalcogenide family of semiconductors. This level of thorough understanding of twistronic phenomena is vital in the search for new material platforms for localized moiré electrons. DOI: 10


I. INTRODUCTION
The number of two-dimensional (2D) materials that have been experimentally isolated as single layers of atomic-scale thickness has been growing at a rapid pace [1][2][3][4][5][6].Alongside these advances, methods to combine single layers into multilayer heterostructures [7][8][9][10] have developed a remarkable level of control over the relative orientation (twist) of successive layers, on the order of 0.1 • .The resulting heterostrucures represent a new type of composite materials with properties spanning a vast range that includes insulators, semiconductors, metals, and superconductors [11].The different arrangements of layers of 2D materials provide an intriguing platform for exploring new physics and potential applications based on their electronic, optical, magnetic, and thermal properties.The relative orientation of successive layers, often characterized by a twist angle between the ideal in-plane lattices, represents an additional "knob" for adjusting the system properties.A well-studied example is twisted bilayer graphene (tBLG), where the twist angle between successive layers of graphene causes controllable interlayer electronic hybridization.At a magic twist angle (≈1 • ) a symmetric hybridization between the two Dirac cones from the individual layers introduces flat bands in the bilayer band structure, resulting in localized Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
In the study of twisted bilayers, the use of a supercell approximation [21][22][23][24][25] or a continuum model [13,[26][27][28] imposes a periodic length scale for the system to provide interpretable band structures.We show that instead of relying on the bands of Bloch theory in momentum space, the local density of states (LDOS) in either the space of atomic configurations or Bloch states provides a surprising amount of clarity in the study of twisted electronic structure.The experimentally important notions of flat bands and moiré band gaps are still obtainable in this context.In addition, the LDOS can form sharply defined features in configuration space and energy, which we call "configuration banding."These two regimes of interesting twistronic features, hosting localized modes or banding, are analogous to the electronic structure observed in conventional molecules (localized states) or crystals (extended states or bands).This pattern arises because of a duality between the momentum and position operators that can occur in specific scenarios for moiré Hamiltonians.
To illustrate the utility of the duality interpretation in a simple context, we first examine a one-dimensional (1D) incommensurate chain model.Moiré flat bands occur when only a finite number of Bloch states or local configurations are needed to capture how two parabolic band extrema are coupled through the interlayer coupling; this corresponds to the "moiré molecule."Configuration banding occurs when an infinite number of Bloch states are needed to accurately capture the interlayer interaction between two band structures; this corresponds to the "momentum crystal," which often exists alongside conventional Bloch bands in the moiré systems.Expanding on the simple 1D example, we then investigate similar features in the electronic structure of realistic materials, such as graphene and MoS 2 , a representative of the transition metal dichalcogenide (TMDC) family.The moiré molecule interpretation can be modeled by a harmonic oscillator and allows for accurate prediction of twist-induced flat bands.The momentum crystal interpretation can also provide highly localized electronic states, and provides an alternate explanation of the localization caused by strong incommensurate potentials.
This paper is organized as follows: Section II contains a discussion of the methods used for modeling electronic structure of 1D and 2D moiré systems, Sec.III categorizes twistronic spectral features for 1D systems, Sec.IV applies this categorization to realistic 2D moiré bilayers, and Sec.V summarizes the results.The appendix examines how configuration banding manifests in the real space LDOS.

II. METHODS
It is a daunting computational task to calculate accurately the properties of layered structures that involve a small twist angle.This is because these atomic arrangements are either incommensurate or have periodicity over length scales several orders of magnitude larger than the primitive unit cells of the individual layers.In tBLG, the Dirac cone yields a simple low-energy Hamiltonian from which momentum-scattering interactions between layers can be computed [13,26,27,[29][30][31].For transition-metal dichalcogenides (TMDC) similar models for the parabolic band extrema have been used [28,32,33].Here, we employ accurate ab initio tight-binding Hamiltonians [34,35] derived from density functional theory (DFT) and maximally localized Wannier functions (MLWF) [36] within a framework specifically developed for twisted 2D systems [37,38].
The tight-binding modeling makes it feasible to investigate the electronic features in twisted 2D materials in two different ways.First, by calculating the band structure through diagonalization of twisted supercell Hamiltonians.Second, by calculating the local density of states (LDOS) of aperiodic systems by computing local spectral properties at the center of large finite regions with the kernel polynomial method (KPM) [38,39].For the latter method, we used circular disks of diameter of 80 nm for TMDCs and of 200 nm for graphene, with both geometries containing millions of atomic orbitals.
In Fig. 1 we introduce an important concept for exploring the physics, namely, displaying the LDOS for atomic structures chosen along lines in configuration space which connect the high-symmetry stackings of the bilayer.This is analogous to displaying band structures in reciprocal (momentum) space along lines connecting high-symmetry points in the Brillouin zone (BZ).The Wigner-Sietz cell of the moiré supercell is outlined in Fig. 1(a), with insets showing the local atomic configuration at selected points.The LDOS results displayed in Figs.1(b)-1(g) are not generated by looking at different atoms in the same structure.Each configuration is calculated independently of all the others by uniformly shifting the top layer relative to the bottom one.
For both tBLG and twisted bilayer MoS 2 (tBLMoS 2 ) near 0 • , there are three high-symmetry stacking configurations: AA, AB, and BA.The sublattices of the honeycomb structure, A and B, are used to describe the alignment between the two layers.AA has both the A and B atoms aligned vertically, while AB and BA have only one pair of opposite sublattice atoms aligned vertically between layers.For graphene, both sublattices are identical carbon atoms, while the hexagonal TMDC hosts a metal on the A site and a dimer of chaclogen atoms on the B site.Note that the AA and AB notation is often used in the TMDCs to instead distinguish between the two distinct phases of an aligned bilayer: the 0 • ± 120 • alignment and 60 • ± 120 • .Although near 0 • twist the TMDC bilayers have the same symmetry as tBLG (with AB and BA stacking related), due to the nonidentical sublattices that symmetry is broken for TMDC bilayers near 60 • twist.
The smooth variations in LDOS at θ = 0 • can be attributed solely to the local variation in interlayer coupling over configuration space [40].For θ = 0 • , localized modes appear in both tBLG and tBLMoS 2 [Fig.1(g)].Their appearance in tBLG at θ = 1.0 • corresponds to the localized wave function caused by moiré flat bands at the so-called magic angle [13].The tBLMoS 2 LDoS of Fig. 1(f) is more complicated.The curved and straight features in the LDOS smoothly expand or contract as the energy changes, caused by surfaces of high electron density in the configuration-energy space.This is an example of electronic banding in the configuration basis.

III. ONE-DIMENSIONAL MODEL
To investigate the nature of twistronic phenomena, we consider a simple 1D model consisting of two chains of single-atom unit cells with a starting lattice parameter of L = 1.The mismatch is introduced by changing the lattice spacing of one layer to (1 − )L, for a small value as shown in Fig. 2(a).encodes a lattice mismatch between the chains, and is the 1D equivalent to the twist angle.Although we focus on twisted structures in the case of 2D bilayers, moiré superlattices generated by a small lattice mismatch show similar phenomena.
For the single-layer electronic structure, a nearest-neighbor hopping of strength T 0 = 1 defines the characteristic energy scale for the system.With two atoms in different layers separated by an in-plane distance d, the interlayer coupling is ( The local density of states (LDOS) for this model for various choices of , ξ , and T 1 is shown in Fig. 2.This LDOS is calculated by the conventional method: diagonalizing a periodic supercell with 10 k-point samples in the 1D Brillouin zone.
Focusing on the = 0 cases, there are a number of common features.Near the extrema of the electronic states, bright lines correspond to parabolic extrema in the one-dimensional band structure.The energy of the band extremum depends on the interlayer stacking d because of the stacking-dependent interlayer coupling.In this simple system, this is caused by a splitting of the two monolayers' identical spectra by the interlayer coupling strength, creating copies of the monolayer band structure that are either shifted down or up in energy (corresponding to bonding or antibonding combinations).
When the atoms of both layers are on top of each other (d = 0 or L) the interlayer coupling is strongest, causing the largest splitting of the two copies.In the two ξ = 0.3 cases, the splitting is minimal as the layers have almost no interlayer coupling at d = L/2.Although a larger value of ξ (longer-range coupling) will generally increase the absolute strength of the interlayer coupling for a given stacking, it reduces the variation in the interlayer coupling as it begins to average over adjacent sites.This is visible when comparing the d variation in the positive LDOS extrema between ξ = 0.3 and 0.6 calculations.Increasing the variation in the interlayer coupling, and not its absolute strength, is the key to inducing moiré phenomena, so generally a small ξ is preferred.
The LDOS of the = 0 systems fall into three general categories.Near band extrema isolated bright spots can occur, corresponding to the emergence of moiré flat bands.Away from band extrema two scenarios are possible.The LDOS can either be featureless in both energy and d, or it can have sharply defined lines crisscrossing one another.The former is the signature of conventional Bloch waves, while the latter is a unique feature of the moiré systems, the phenomenon we defined as configuration banding.In the following subsections we will investigate each of these categories with the 1D model, explaining how they arise and illustrating the duality between the configuration and Bloch bases in twistronic materials.

A. Flat bands from a harmonic approximation
The moiré flat bands can be explained in this model by performing a perturbative expansion in both the monolayer eigenvalues E (k) and the interlayer coupling T (d).As we are concerned only with the band extremum in this case, the monolayer Hamiltonians can be modeled as a uniform electron gas with effective mass m * , where σ = ±1 depends on if we are interested in a valence (hole) band maximum or conduction (electron) band minimum.The interlayer coupling can be approximated as a perturbing potential in space T (r).When the band extremum is located at the point (k i 0 = 0), the full Hamiltonian can be written entirely in the position basis Defining both a bonding and an antibonding pairing between the layers, , separates this into two singlevariable Hamiltonians: If the interlayer coupling can be locally expanded at its extrema in a quadratic form T (r) = σ 2 m * ω 2 r 2 , one obtains the harmonic oscillator equation This concept is illustrated in Figs.3(a) and 3(b) for a holelike band (σ = −1), and the effective parameters m * and ω are easily extracted from the simple 1D model.The band structure for a 1D chain with nearest-neighbor coupling T 0 and lattice parameter L is given by 2T 0 cos kL Assuming is small, we can take m * = (2T 0 L 2 ) −1 for both chains.Here, we have kept terms of T 0 and L (both equal to 1) to ensure a general result.The value of ω depends on both the strength of interlayer coupling T 1 and the moiré length λ = 1/ .Making the substitution d = r/λ in the definition of T (d ) gives The -independent part of the above expression, ω , is a useful parameter when analyzing moiré flat bands that arise from parabolic extrema.
Comparing the expected energy levels of the harmonic oscillator ( 12 + n)ω to the calculated DOS in Fig. 3(d) shows good agreement for small ( 0.1).For larger values, the electronic density fans away from the ideal oscillator modes as the width of the associated moiré bands grows.For ideal "flatness" of the nth oscillator state, we find the condition to be (n + 2.5)ω < T 1 , or at least three harmonic oscillator states must fit in the interlayer potential well, as shown in Fig. 3(e).

B. Moiré molecules
Beyond its use to predict moiré flat bands, the harmonic oscillator equation also serves as a clear example of momentumconfiguration duality.The two terms in the harmonic equation, quadratic in either the momentum operator or the position operator, are interchangeable provided their prefactors are comparable.For the momentum term, the prefactor was simply the curvature of the electronic bands (i.e., 1/m * ): For the position term, it was the curvature of the interlayer coupling potential Putting these together, the frequency of the harmonic oscillator is This illustrates the dual nature of momentum and configuration spaces in moiré problems.Indeed, we can re-examine the LDoS calculation to see the duality more clearly.In configuration space, we computed the LDOS as a function of configuration, but implicitly summed over all the momentum degrees of freedom for the moiré supercell: Reformulating this equation under the duality, we can instead calculate the LDOS in momentum space by summing over all configurations for a specific moiré wave number: In Fig. 4, this is presented for the positive band extremum of the 1D model.Clearly defined states are isolated in energy in both the configuration and momentum bases.If the curvature of E (k) is equal to that of T (d), the two sets of figures would be indistinguishable near the band extremum.In both configuration or momentum space, a strong dependence of the spectrum on the local variable (d or k) causes pockets of isolated eigenvalues in the energy spectrum.These couple to one another, forming an effective molecule and creating electronic structure reminiscent of molecular orbitals in both configuration and momentum space.The higher-energy harmonic states have more nodal points, corresponding to higher "momentum" in 1D.This concept generalizes to 2D crystals, where the higher-energy moiré molecular states have larger angular momentum (s wave, p wave, etc.)

C. Moiré crystals
The question arises, when and how does the previous expansion for ω and the moiré molecule fail?Such a failure is visible in the low-energy regions of Fig. 4 where the smooth density of a banded Bloch wave appears.This occurs exactly at the energy when the d-dependent extremum in the LDOS reaches its minimum and changes direction.There is no longer any spatial confinement of the spectrum, but more importantly T (d) can be treated as a constant in this energy range.From the Bloch-state point of view at this energy, T (d) is constant (and thus can be neglected), reverting to the band structure of a typical 1D crystal.This is illustrated by the schematic picture and LDOS calculations in Figs.5(a)-5(c).The standard result is obtained for the LDOS in configuration space (a continuum) and momentum space (a band).
Connecting this to the derivation of ω , we see that the harmonic approximation fails because ∂ 2 T ∂d 2 has vanished.Taking T 1 = 3T 0 and ξ = 0.3 [Figs.5(g)-5(i)] gives an interesting situation: both the conventional and unconventional crystals occur in the same spectrum.Isolated spots of density, associated with moiré flat bands, occur in both spaces.The two crystals occur, showing a continuous spectra in the original space and a clear band structure in its dual.Comparing this to the = 0 calculation of the same parameters in the lower-left panel of Fig. 2(c), the conventional crystal regime is predicted by a d-independent energy region of the LDOS.For the unconventional crystal, there is an equivalent k-independent energy region in the momentum-projected LDOS at = 0.

IV. TWO-DIMENSIONAL BILAYERS
We now turn our attention to moiré crystals in 2D.First, we will show the ubiquitous ability for the configurationdependent LDOS in the untwisted geometry to predict localized electronic structures in twisted geometries by looking at examples of tBLG and tBLMoS 2 .Second, we will show how conventional band structure calculations of bilayer MoS 2 shed light on how the moiré molecule interpretation can describe the appearance and properties of moiré flat bands in generic twisted semiconductors.

A. Predicting twistronic features from LDOS
The three distinct regimes of electronic structure in the 1D model, the moiré molecule, and the real and momentumspace moiré crystal, can be directly predicted by careful interpretation of the untwisted LDOS.The moiré molecule, and its associated flat bands, occur when wells appear in the moiré potential.For a 2D crystal, this manifests in the LDOS as extrema or bright lines that have strong configuration dependence.The conventional moiré crystal, hosting traditional Bloch wave states, occurs in regions of the LDOS that have little configuration dependence.The unconventional moiré crystal and its configuration banding require a region of LDOS that has little k dependence after averaging over all configurations.This manifests in the configuration LDOS as regions with steep configuration dependence.
Following these rules, the configuration dependence of the LDOS in untwisted bilayer graphene [Figs.6(a the density of states.This singular feature (of the monolayer) is split by the effective interlayer coupling in the bilayer, and thus has strong configuration dependence.The largest splitting occurs at AA, with slightly weaker splitting at AB and BA, consistent with the microscopic origins of the interlayer coupling.
Looking at the 1 • twisted bilayer, a singular feature appears at 0 eV and AA stacking: this is precisely the well-known magic angle flat band.It can be interpreted in the same fashion as solutions of the harmonic continuum model by replacing the quadratic monolayer Hamiltonians with the conical Dirac spinor [13], which still provides analytic solutions [41].At 2 • , there are two bright features at low energy, and these correspond to the Van Hove singularities of the twisted bilayer.These features arise from the conical moiré potentials visible in the untwisted LDOS.
For the two other energy regions, twisting causes isolated levels to occur in each of the T (d) wells of the untwisted LDOS.In particular, between AB and BA a rather steep harmonic moiré potential well (relative to the monolayer band curvature) generates many stable harmonic oscillator levels, creating checkerboard patterns in the LDOS.Doubling the twist angle doubles the energy spacing between the states, showing that the harmonic approximation is still in good agreement with these features.Comparing the relative spacing between the states centered at AA and the midpoint (between AB and BA) we see that the curvature of T (d) in the untwisted case proportionally affects the spacing between levels.This will be more thoroughly validated in the following subsection.
Similar behavior can be observed in our calculations of twisted bilayer MoS 2 (tBLMoS 2 ).In Fig. 7 we show the LDOS intensity of the p z orbitals of sulfur as a function of energy and local atomic structure.These orbitals extend the farthest out of the plane of each layer and have the strongest interlayer coupling, so we have ignored other orbital contributions to the LDOS here.The other orbitals do couple weakly with the p z states in certain energy regions, and thus show faint signatures of twistronic features, but they primarily consist of a smooth LDOS background due to uninteresting Bloch waves (the conventional moiré crystal regime).
The first energy region is the valence band maximum [Figs.7(a has singular features in the LDOS.This region arises from a Lifshitz transition in a p z -majority band, and is functionally identical in shape to the feature in the graphene model.This implies that this particular twistronic structure is a signature of p z bands on a triangular lattice.As it occurs in both graphene and MoS 2 , the structure does not seem to require sublattice symmetry on the honeycomb lattice.Similar to graphene, moiré molecular states are visible at each of the band extrema.In contrast, weak configuration banding is visible near the center of the feature (−2.0 eV), which was missing in tBLG, due to stronger interlayer coupling relative to the monolayer band curvature in the TMDCs.This agrees with the requirement of a sufficiently strong interlayer moiré potential in the 1D model [Fig.2(c)].
The last energy region [Figs.7(g)-7(i)] again shows LDOS features consistent with a 2D band extremum, however, the configuration dependence is much stronger than at the valence band maximum.Accordingly, clear moiré molecular states exist near −4.4 eV, and more robust configuration banding occurs at the center (−4.0 eV).Two independent sets of AAcentered harmonic oscillator states appear more clearly in this last case, owing to the two AA-centered T (d) wells at −4.4 eV with different curvature.The LDOS corresponding to the n = 0 and 1 harmonic energy levels of the shallower moiré potential are shown in Figs.7(l) and 7(m) displaying s and p x ± ip y symmetries, respectively.In Fig. 7(k) a bright center inside a ring is visible, a signature of d orbital momentum from the n = 2 level of the sharper well.At slightly higher energy (−4.3 eV) a T (d) well centered in between AB and BA also hosts moiré molecular states, although they are not as well confined.The coexistence of multiple unique moiré flat band modes at the same energy could give rise to unique quantum phases.
These results introduce an additional constraint necessary in the classification of possible twistronic phenomena of 2D materials: the connectivity of the constituent monolayer band structures.The moiré molecule and its associated flat bands requires a small number of connected Bloch states of similar energy and, therefore, it naturally occurs in parabolic band extrema and Dirac points.The momentum-space crystal can occur at 2D Lifshitz transitions in the monolayer band structure.A previous work [42] showed that these transitions can create an infinite lattice of coupled Bloch states, preventing convergence of the truncated continuum models for twistronic systems.This describes the formation of a momentum-space crystal, albeit only on a subset of the Brillouin zone, and the failure of the continuum models is directly related to the breakdown of the dual nature of the moiré molecule.This means saddle points in the 2D band structure are good candidates for the realization of configuration banding.An alternate approach would be to find a system where the interlayer coupling T is larger than the entire width of a monolayer band.In van der Waals materials, the interlayer distance between orbitals is larger than the in-plane distances, so this is a challenging task.
The Lifshitz regions also host robust moiré flat bands at the interlayer coupling extrema, which justifies experimental study of twisted bilayers made of metallic materials.Unlike flat bands arising from Dirac cones or band extrema, it seems necessary for these to exist within a background of metallic states.The LDOS calculations are particularly useful for identifying these modes, as band structure calculations are difficult to interpret within a metallic background and continuum methods are unlikely to converge [42].

B. Flat bands in MoS 2
The harmonic approximation of band extrema for the incommensurate 1D chain is a general interpretation of any moiré flat bands arising from interlayer hybridization.This approximation can also give reliable predictions of band flattening at the valence band maximum of MoS 2 .These results are in agreement with the LDOS calculations of the previous section, but by performing band structure calculations the energy resolution limitation is removed, allowing for careful comparison to the harmonic oscillator.This band extremum (maximum) is at and is maximized at AA stacking, as shown in Fig. 8(a).Accordingly, we take this maximum value as E = 0. Following the 1D case, we then extract the effective mass m * from a monolayer calculation and the configuration dependence of the band extremum energy, which provides the effective moiré potential T (d).Both are obtained by fitting up to the quadratic term to extract the effective harmonic frequency parameter ω θ = −7.24meV per degree.For a given twist angle, the harmonic frequency is given by ω = ω θ θ .
The band structure calculations show excellent agreement to the energy levels of the effective harmonic oscillator ω(n + 1).Unlike the 1D case, the nth 2D harmonic oscillator level is (n + 1)-fold degenerate.Looking at the band structure of the θ = 2.28 • supercell (λ = 8 nm) in Fig. 8(b), we find a single flat band with nearly zero bandwidth at the expected value of ω θ θ = −16.5 meV.Lower in energy, a pair of nearly degenerate flat bands have a value just shy of the next energy level 2ω θ θ = −33 meV.Below this, a set of three bands have separated from the remainder of the valence band manifold, although they still have a substantial bandwidth.This analysis is performed for each supercell and we report the band average energy and total bandwidth along the high-symmetry lines in Fig. 8(c).Near 2 • , the first flat band is in perfect agreement with the QHO model, but as θ increases the band becomes less flat and its average begins to drift from the QHO model.This same behavior can be seen for the n = 1 and 2 levels, and so in the θ → 0 • limit we expect every n will emerge as a flat band in this simple tight-binding model.Following a similar constraint of the band flattening in the 1D model, we find the nth level will become sufficiently flat (b w < 5 meV) when the (n + 2)th level is contained within the moiré potential well.The moiré potential for the MoS 2 model used here is approximately 100 meV deep.That means the critical twist angle for the n = 0 flat band should satisfy |(2 + 1)ω θ θ (0)  c | < 100 meV, yielding θ (0) c ≈ 4.6 • .For the second and third levels of the oscillator, this similarly provides critical values of θ (1)  c ≈ 3.5 • and θ (2)  c ≈ 2.8 • .These are all in excellent agreement with the results of the tight-binding calculations, namely, that when θ < θ (n)  c the bandwidth for the nth level is below 5 meV.Comparing this result to full DFT calculations [23] we find good agreement to the unrelaxed case: the flat bands are associated with AA stacking spots at the valence band maximum.However, the inclusion of atomic relaxations changes the valence flat-band characteristics.This is primarily due to large variations in the interlayer separation as a function of atomic configuration [23,25], which significantly modifies the interlayer moiré potential.The predicted moiré potential also depends sensitively on the choice of van der Waals DFT functional [43].For this reason, the results presented in Fig. 8 are meant to serve as an introduction to the band flattening phenomena in twisted semiconductors, not as a fully selfconsistent prediction of flat bands.See Refs.[23,25] for fully ab initio predictions of band flattening in semiconductors, including important relaxation effects.
Although the derivation of an explicit harmonic oscillator equation relies on the maximum being at , it is not a necessary condition for the harmonic oscillator interpretation of the flat bands.Any band extremum in any material can host moiré flat bands, with its own characteristic twist angle depending on the interlayer coupling strength and effective mass.From this perspective, the Dirac cones of graphene are an exception to the rule.The two-sided nature of the Dirac equation causes the flat-band condition to occur only under a fine tuning of the twist angle to a magic value.For generic insulators or semiconductors, this fine tuning is not necessary, as reducing the twist angle always causes a shallower moiré potential and leads to more confined harmonic states and thus flatter bands [44].
In experiment, one wants not only to minimize the bandwidth b w , but also to keep each set of flat bands separated sufficiently in energy.In the unrelaxed model studied here, it is clear that although one can make the bandwidth arbitrarily small as θ → 0 • , the band gaps (given by ω θ θ ) will vanish.
A tradeoff between bandwidth and band gap should be considered depending on the goal of the proposed experiment.For angles well above the first critical angle θ (1)  c ≈ 4.6 • , no bands are guaranteed to be flat or gapped from the rest of the valence manifold.If more than one harmonic oscillator level is required, or a certain minimum bandwidth is desired, this will set an upper bound for the desired twist angle.In both cases, any constraint on the minimum allowed band gap will set a lower bound on the twist angle.Unlike twisted bilayer graphene, variations in the twist angle [45] are not likely to fundamentally alter the correlated behavior of electrons in the material, as the band flatness is generally monotonic with respect to the twist angle.
Additional corrections should be carefully considered in the future.The anharmonicity in the effective mass and interlayer coupling will be necessary for accurate comparison to experimental results, but the overall intuition for the moiré flat bands is unchanged.In addition, in-plane relaxations at small twist angles form domain-wall structures [46][47][48][49] and will cause the moiré potential to become independent of the twist angle [50].This prevents the tuning of ω to arbitrarily small values, and will define a minimum possible bandwidth for the flat bands.As the bands of MoS 2 are primarily p z character and thus have negligible spin splitting, spin-orbit coupling was not needed here.If the moiré bands instead arise from the K valley in a TMDC material other than MoS 2 , spin-orbit coupling can be important [28,33,44].

V. CONCLUSION
We have provided a classification of twist-induced electronic structure in 1D and 2D moiré systems through the local density of states.The self-dual nature of the moiré Hamiltonians puts the momentum and configuration bases on equal footing, providing insight into twistronic features.In particular, moiré flat bands can be robustly predicted from either the LDOS or band structures of untwisted systems.The LDOS approach also allows for prediction of flat bands even in metallic backgrounds, and explains the formation of effective momentum crystals, problems not yet addressed by band structure methods.The study of moiré heterostructures made of arbitrary 2D materials is made much easier within this framework.Apart from the strengthening of correlated effects, many applications that take advantage of the electronic localization in these structures can be envisaged, such as arrays of quantum dots or networks of 1D electron channels.1D model [Fig.5(g)], but one may ask what these electronic states look like in real space.Calculations of the LDOS for a system with a momentum-space crystal are presented in Fig. 9 for various choices of .It begins with a commensurate cell, 0 = 0.1.
To understand the role of incommensurate systems in the context of configuration banding we introduce incommensuration through the parameter , setting = 0 × (1 + ).must be irrational to generate a truly incommensurate system.However, in calculating the d-dependent LDOS via the finite-sized method (for both 1D and 2D bilayers), we find that small values of have no qualitative effect when 0 is small.The d-dependent LDOS is a stable quantity under small variations in the twist angle.For this reason, in the following discussion and in Fig. 9 we always use the configuration LDOS calculation for = 0 , and have only affect the sampling of d to generate the real-space LDOS.
The interlayer configuration at the origin (r = 0) is controlled by the parameter d 0 .For commensurate systems, this choice of d 0 is only unique up to 1/λ, as units of the inverse moiré length correspond to a translation of the origin.Modifying d 0 will show the relative stability of the real-space LDOS under small changes in interlayer alignment in commensurate systems.
For all choices of , the high-and low-energy regions show isolated electronic states that are robust under changes in the initial configuration d 0 and incommensuration .These are the harmonic oscillator states of the moiré molecule, and their stability with respect to is explained by the insensitivity of the previous harmonic expansion to the microscopic details of atomic geometry.
For commensurate geometries, a slight relative movement of the layers can drastically alter the electronic structure in the momentum-space crystal regime.Near zero energy in Figs.9(b) and 9(c) the energies of the localized modes change significantly when the configuration is moved by half of .This is because a commensurate system only samples a finite number of configurations, and so changing the initial configuration d 0 causes different sections of the configuration bands to manifest.A visual aid for inferring the real-space LDOS from the configuration space LDOS is a "moiré comb" which has a spacing between each tooth, as shown in Fig. 9(a).The two red combs correspond to two values of d 0 which generate geometries with vastly different LDOS near zero energy.Changing the initial configuration of the commensurate system corresponds to moving the comb left or right, while changing the energy of interest moves it up or down.The electronic structure of these commensurate geometries in the momentum crystal energy range looks markedly similar to the stable isolated states of the moiré molecule, but they are highly sensitive to small changes in alignment between the layers.Care must be taken in moiré supercell calculations to avoid interpreting these features as stable flat bands.
Incommensurate geometries, or geometries which are not periodic on the moiré length, with = 0 can be interpreted as a sequence of moiré supercells of length 1/ 0 but with the variable d 0 varying in space.The isolated electronic states that result band with a wavelength determined by a second-order moiré length, which is the length required for d 0 to go once through the range [0, L].This electronic structure hosts highly localized wave functions, which are not directly provable from the LDOS calculations, but have been studied extensively for incommensurate potentials acting on 1D tight-binding models [51][52][53].In particular, self-dual models of this form [51,54] bear strong resemblance to the momentum-space crystal structure studied here.Making a direct connection between the moiré system and these incommensurate potential models via explicit comparison of the tight-binding models is difficult, owing to the inability to pair atoms of the top and bottom layers together in a consistent manner across the surface of a moiré pattern [55].

FIG. 1 .
FIG. 1.(a) Geometry of a twisted bilayer formed by two identical honeycomb lattices.The bottom (top) layer is red (blue).The moiré supercell is shown in gray and the moiré Wigner-Seitz cell is outlined in black, with colored boxes highlighting the high-symmetry stacking locations shown in the insets, corresponding to AB, AA, and BA stackings.(b) Configuration-dependent local density of states (LDOS) for p z orbitals in untwisted bilayer MoS 2 in an energy region 4 eV below the valence band maximum.(c) Configuration-dependent LDOS for a 4 • twisted bilayer of MoS 2 .(d), (e) The same as (b) and (c), but for twisted bilayer graphene at 0.0 • and 1.0 • .(f), (g) Top-down views of the real-space variation in the LDOS at the energy values identified by the dashed colored lines.The LDOS in (g) appears in both tBLG and tBLMoS 2 , and includes the moiré Wigner-Seitz cell outlined in white.

FIG. 2 .
FIG. 2. (a) The 1D two-chain model, with the top (bottom) chain in red (blue) with its lattice parameter in the same color.The relative interlayer configuration between the top and bottom layer (d) and its interlayer coupling [T (d)] are labeled for a pair of atoms.(b) The interlayer coupling function T , with parameters ξ = 0.3 (0.6) and T 1 = 3 (1) in black (red).(c) Configuration-dependent LDOS of the top atom in the two-chain model.Each of the four panels consists of an "untwisted" geometry with = 0, and two "twisted" geometries with = 0.05 and 0.1, for two different values of ξ , ξ = 0.3 (left set), and ξ = 0.6 (right set), and two different values of T 1 , T 1 = 1 (top set), and T 1 = 3 (bottom set).

FIG. 3 .
FIG. 3. (a) Band structure of a single 1D chain, with a harmonic approximation of a band extremum in red.(b) Configuration dependence of the interlayer coupling function for the bilayer 1D chain model, with a harmonic approximation of its maximum in red.(c) Total density of states (DOS) for commensurate bilayer chains for in [0, 0.4].A fractal butterfly spectrum is evident, with a region highlighted in the top left corner.(d) Expanded view of the region highlighted in (c), that shows a number of isolated -dependent states in the DOS.The dashed red lines correspond to the energy levels ω(n + 1/2) of a real-space harmonic oscillator based on the red fits in (a) and (b).The red numbers indicate the n for that level.(e) Schematic of the fitted real-space harmonic well, with the first three energy levels of the oscillator indicated by the red lines and numbers.The frequency ω = ω is displayed between the first two energy levels.

FIG. 4 .
FIG. 4. Electronic structure calculations for bilayer chains with T 1 = 3 and ξ = 0.3, with varying .Top: LDOS for atoms of the top layer with specific configurations d.The harmonic approximation of T (d) is highlighted in red.Bottom: LDOS for Bloch waves of the top layer with specific momenta k.The harmonic approximation of the bands E (k) is highlighted in red.

FIG. 5 .
FIG. 5. (a) Diagram of the bilayer chain Hamiltonian when the monolayer couplings (H i ) are much larger than the interlayer coupling (T ).(b) Configuration-dependent LDOS for the atoms of the top layer in this limiting case with = 0. (c) Momentum-dependent LDOS for Bloch waves of the top layer in this limiting case with = 0. (d)-(f) Same as (a)-(c), but the interlayer coupling is now much larger than the monolayer couplings.(g)-(i) Configuration and momentum LDOS for the bilayer chain model ( = 0.1, T 1 = 3, ξ = 0.3) with the central panel describing the twistronic regime at that energy.

FIG. 6 .
FIG. 6. (a)-(i) Configuration-dependent LDOS for twisted bilayer graphene (tBLG) for selected energy regions and various twist angles.The configurations correspond to a diagonal traversal of the moiré supercell.The results were calculated for finite flakes of diameter 200 nm with a kernel polynomial method.(j) Band structure for a simple nearest-neighbor tight-binding model of graphene.The energy regions associated with the selected LDOS energy regions of the twisted bilayer are highlighted in red.(k) Top-down view of (d).The Bloch states that are within the interlayer coupling strength (0.3 eV) of each selected energy are highlighted in red.
FIG. 7. Configuration-dependent LDOS for the p z orbitals of bilayer MoS 2 in selected energy regions of the valence bands.The configurations go across the diagonal of the moiré cell, starting and ending at aligned AA-type stacking regions.(j)-(m) Show a top view of the configuration-dependent LDOS at four energy values for the 4 • twisted bilayer, interpolated from a 20 × 20 sampling of configuration space.All results were calculated for finite flakes of radius 40 nm with the kernel polynomial method.The energy spectrum resolution is limited by roughly 10 meV (p = 4000).
)-7(c)].The LDOS shows a step-function behavior in the vertical (energy) direction, consistent with the density of states of two-dimensional parabolic bands.At 4 • , signatures of moiré flat bands are visible, but hard to resolve as the KPM broadens the eigenvalues by a Gaussian of width σ = 5 meV in these calculations.The next energy region [Figs.7(d)-7(f)]

FIG. 8 .
FIG. 8. Interpretation of the moiré flat bands in unrelaxed twisted bilayers of MoS 2 as a quantum harmonic oscillator model.(a) Band structure comparison between a monolayer (orange) and AA-stacked bilayer (black) of MoS 2 .The effective interlayer coupling can be easily interpreted at any band and momentum by looking at the strength of the band splitting.The dashed red line indicates the bilayer's valence band maximum, which is the zero-energy reference for the following panels.(b) Band structure of the twisted bilayer at θ = 2.28 • , showing the characteristic splitting and flattening of the bands, as captured by the QHO model.(c) Summary of band structures for twisted bilayer MoS 2 near the valence band maximum.The average band energy for each band is plotted as a function of the twist angle θ, with the bandwidth b w given in color.The red lines indicate the expected energy levels of a harmonic oscillator defined by the monolayer's effective mass m * , and the shape of the moiré potential T (r), both expanded around .The dashed lines are a guide to the eye for θ dependence of the first (singly degenerate) level and the second (double degenerate) level.

FIG. 9 .
FIG. 9. (a) LDOS for the 1D bilayer model as a function of local configuration d at = 0.1.For a commensurate geometry, only a finite number of configurations are sampled.The sampled configurations for two different choices of d 0 are indicated by the red "moiré combs."(b), (c) Real-space LDOS for the two commensurate geometries (C) introduced in (a).The configuration banding is lost due to the finite sampling of configurations.The moiré length λ is given at the bottom of (c).(d)-(f) Real-space LDOS for three selected incommensurate geometries (I), parametrized by = 0.1 × (1 + ).As these geometries are not periodic on the moiré length, the configuration bands are visible in real space.