Nematic stripe density wave order and Mott insulating ground states in small-angle twisted bilayer graphene

In this work, we determine states of electronic order of small-angle twisted bilayer graphene. Ground states are determined for weak- and strong-couplings being representatives for varying distance of the twist-angle from its magic value. In the weak-coupling regime, charge density waves emerge which break translational and $C_{3}$-rotational symmetry. In the strong coupling-regime, we find rotational and translational symmetry breaking Mott insulating states for all commensurate moir\'e band fillings. Depending on the local occupation of superlattice sites hosting up to four electrons, global spin-(ferromagnetic) and valley symmetries are also broken giving rise to characteristic Landau level degeneracies in detailed agreement with observations for commensurate band fillings of $\nu=0,\pm1/2,\pm3/4,$. The formation of those particular electron orders is traced back to the important role of characteristic non-local interactions which connect states sharing one hexagon of AB- and BA-stacked regions of the superlattice.

Quantum oscillations reveal that the insulating states differ in their Landau level degeneracy [1][2][3][4]. This indicates the presence (or absence) of (global) symmetries which generate Kramer-like degeneracies of singleparticle states. In our description we assume the presence of spin-rotation and valley symmetry which is justified in the limit of vanishing spin-orbit coupling [42] and small twist-angles [15]. This implies for the insulator state at charge neutrality (ν = 0), which is found to be 4-fold degenerate, the presence of both spin-and valley-symmetry, for half electron-or hole-filling (ν = ±1/2) being twofold degenerate that either spin-or valley-symmetry is broken, and for band filling ν = ±3/4 being single degenerate that both spin-and valley-symmetry are absent [2][3][4]. Those finding put strong requirements on possible * markus.klug@kit.edu Figure 1. By tuning the twist-angle, twisted-bilayer graphene undergoes a transition from a weak into a strong-coupling regime. The moiré bandwidth (red) represents the kinetic energy scale, the amplitudes of interaction matrix elements of Eq. 4 (green) the potential energy scale. Their ratio (black) grows towards the magic-angle regime signaling the crossover. ground states of insulating states.
In addition to transport experiments [1][2][3][4][5], which were performed on samples identified with the so-called magicangle regime [6] hosting among others Mott-like insulator states and superconductivity, STS and STM experiments revealed correlated electron states for large doping ranges around charge neutrality (CN) [7][8][9][10]. Here interaction effects manifest themself as significant redistribution of single-particle spectral weight setting in at certain electron-or hole-fillings of moiré bands. Additionally, it is unanimously reported that the correlated states break C 3 -rotational symmetry as seen in spatially resolved charge distribution measurements and that the effect is largest at CN. Those results were obtained for samples with moiré bandwidths significantly larger than 10meV perhaps placing them in a close-to-magic-angle regime, where correlation effects are expected to be still present yet with interaction energies smaller than in the magic-angle regime. Notice, the role of interaction effects is in general not determined by the twist-angle alone. As the magic angle regime is determined by the dimensionless quantity α = 3 8π at ⊥ vF sin θ/2 , which takes on the value α m ≈ 1/ √ 3 in the magic angle regime [6], the bandwidth of the single-particle spectrum and thus the role of interaction effects depend on the twist-angle θ, the interlayer coupling t ⊥ , the graphene lattice constant a and the bare Fermi velocity v F . A graphical depiction of the crossover from weak-to strong-coupling by tuning the twist-angle (or equally the interlayer coupling) is given in Fig. 1.
In this work, we investigate possible electronic ground states in the presence of finite-range Coulomb electronelectron interaction as function of doping. We will distinguish between two parameter regimes, a strong-coupling regime α = α m representative for magic-angles where interaction effects dominate, and a weak-coupling regime |α − α m | > 0 where interaction effects are weak. As outlined above, both regimes are likely to be realized in experiments.
In a first step, we set up a single-particle continuum theory for TBG as discussed in Refs. [6,43], with superlattice translational T , crystalline point group D 6 , spin SU (2), and valley conservation U v (1) approximate symmetries [16]. Interaction effects are considered by employing a two-orbital model introduced in Refs. [12,13], where maximally localized Wannier orbitals centered at AB-and BA-stacked regions of the superlattice are constructed from the four moiré bands. Subsequently, interaction matrix elements between the Wannier states in the direct, exchange and pair-hopping channel are computed. Interaction processes are found to be highly nonlocal connecting, in leading order, all orbitals sharing one hexagon of AB-and BA-stacked regions of the superlattice.
In the weak-coupling regime, our results, which are ob- Figure 2. In the weak-coupling regime, our ground state analysis reveals a formation of stripe charge density wave orders which break translational T and C3-rotational symmetries.
In the strong-coupling regime, we find three types of Mott insulating ground states (a, b, c) for all commensurate moiré band fillings ν = 0, ±1/4, ±1/2, ±3/4, which resemble the stripe-type orders in the weak-coupling limit (a) but differ, in part, by the absence of spin SU (2) (b) as well as valley Uv(1) symmetry (c).
tained in an unrestricted mean-field analysis, reveal a formation of stripe charge density wave order with commensurate ordering vectors of half of the reciprocal lattice vectors of the moiré superlattice. The charge inhomogeneities form when the interaction strength and the moiré band filling reaches a certain threshold. However, we find no indication for a link to the van-Hove singularities of the single-particle spectrum.
In the strong-coupling regime, we perform an infinite coupling limit where we drop the kinetic part of the theory and consider only interaction processes which can be expressed in terms of density-density interactions. This approach would be futile in the limit of local Hubbard interactions. However, for the present model, the important role of non-locality combined with ferromagnetic and ferrovalley exchange interactions distinguishing between spin and valley number allows us to determine the nature of ordered states even without kinetic energy contributions of electrons. For all commensurate band fillings ν = 0, ±1/4, ±1/2, ±3/4, we find Mott insulating ground states which break translational T and rotational C 3 symmetry (ν = 0), as well as spin SU (2) symmetry (|ν| = 1/4, 1/2, 3/4) and valley U v (1) symmetry (ν = ±1/4, ±3/4), and being of stripe-type order for |ν| < 1/2, see Fig. 2.
Microscopic model : We obtain the weakly dispersing moiré bands by considering a continuum model [6,43] where states near the slightly twisted Dirac cones of two graphene layers hybridize due to a finite inter-layer coupling. Because of the large separation in momentum space, states near non-equivalent cones, in the following labeled by the quantum number ξ = ±, are effectively decoupled generating an emergent U v (1) valley symmetry. The effective Hamiltonian is therefore given by where band λ, spin σ, valley ξ and crystal momentum k, element of the tiny moiré Brillouin zone, label the superlattice Bloch states ψ λkσξ . Those four moiré bands are separated by significant band gaps from a continuum of bands at higher and lower energies, which is phenomenologically modeled by taking lattice relaxation effects into account [13].
Wannier basis: We construct the Wannier basis by employing a two-orbital model [13][14][15] where localized Wannier functions, though centered at the AB-and BAstacked regions, possess highest weight at the AA-stacked regions of the underlying TBG superlattice. Despite the occurrence of Wannier obstructions, which render certain exact symmetries non-local [15,16] and are resolved by incorporating auxiliary bands [17,18] adding a vast number of degrees of freedom, we assume that the relevant physics of our work is captured by the two-orbital model. This assumption is further supported by the fact that topological properties are irrelevant when we discuss the limit of strong couplings.
In particular, we follow the approach outlined in Ref. 2 sin θ/2 for θ = 1.05°. All interactions connecting sites which share one hexagon are found to be relevant. Longer-distance interactions are numerically small and are therefore negligible. Since the twist-angle dependence of interaction matrix element is rather weak and approximately determined by the superlattice length scale U abcd ∝ L −1 M ∝ sin(θ/2) as depicted in Fig.  1 and being complementary to Ref. [19], the numerical values are representative for both the strong-and weak-coupling regime.
[13] and construct orbitals centered at the AB-and BAstacked regions having definite valley number. More specifically, the Wannier function with valley ξ and spin σ number located in superlattice unit cell i centered at the α = AB,BA-stacked region is given as superposition of Bloch wave functions by where the superlattice vector A i , and the unitary matrix (U (α) λk ) is chosen such that its spread functional is minimized [44,45]. Details about our construction of the maximally localized wannier basis is found in Appendix A. Eventually, we obtain orthogonal and exponentially localized Wannier orbitals, which are labeled by the single-particle quantum numbers a = (i, α, ξ) and σ giving rise to 8 states per superlattice unit cell. The effective tight-binding model is given by where the transition amplitudes {t ab } are, by construction, diagonal in valley number and obtained by projecting the Bloch states into the Wannier basis.
Interaction matrix elements: We compute interaction matrix elements between Wannier states using an unscreened Coulomb kernel, with electron charge e and relative permittivity ǫ ≈ 7 for hexagonal boron nitride. We distinguish between density (a = d and b = c introducing U = U abba ), exchange (a = c = b = d introducing J = U abab ), and pair-hopping (a = b = c = d introducing X = U aabb ). We also checked for charge-bond interaction matrix elements (a = d = b = c), but find them at least one order of magnitude smaller than comparable processes and are therefore safely neglected. We find that the amplitude of interaction matrix elements drops with distance but remains significant for processes connecting all orbitals which share one hexagon formed by AB-and BA-stacked regions of the superlattice, as depicted in Fig. 3. In general, interactions are dominated by direct interaction processes. Since this type of interaction is insensitive to local valley or spin configurations, we expect that possible charge modulations are determined by U . Though at least one order smaller in amplitude but being sensitive to valley and spin number, intra-and inter-valley exchange J and J IV , as well as intra-and inter-valley pair-hopping X and X IV processes are expected to be relevant. Due to rapid phase fluctuations, inter-valley processes are much smaller. However, because of the coupling of otherwise decoupled valley sectors, they are considered relevant to determine the exact ground state. For both the intraand inter-valley case, all matrix elements are found positive J, J IV > 0 causing neighboring spin-and orbitaldegrees of freedom to align. Concluding, the hierarchy of amplitudes of interaction matrix elements, is characteristic for small-angle TBG. Ground state analysis: The effective two-orbital model Eq. (3) established in the previous sections constitutes the basis for the subsequent ground state analysis, which is two-fold: We first investigate the weak-coupling regime (I) for |α − α m | > 0, and second, conduct a strong coupling analysis (II) for α = α m .
Weak-coupling regime (I): For small coupling strength, we conduct a mean-field analysis in which interaction terms are locally decoupled in all possible channels, This ensures that our analysis is susceptible to various kinds of electron instabilities. The mean-fields are considered as variational parameters which are optimized to minimize the ground state energy functional. Spatial correlations c † aσ c bσ ′ are here determined numerically using Lanczos' method [46]. Mutually independent mean-fields are introduced for a sublattice of 6 × 6 superlattice unit cells with imposed periodic boundary conditions which is embedded in a lattice of 30 × 30 unit cells to capture the rather slowly decaying transition amplitudes {t ab }.
Results (I): For a large amount of electron-or holedoping, we find no electronic order (recall, we do not probe for superconductivity) and the characteristic density of states (DOS) at van-Hove singularities as well as the linear dispersion relation near charge neutrality are observed. In contrast, for moderate electron-or holedoping, we observe the formation of a stripe charge density wave order which breaks translational and rotational C 3 -symmetry. The corresponding order parameter is given by with particle number operatorn aσ = c † aσ c aσ and lattice site vector R a . Possible ordering vectors are Q ∈ {G 1 /2, G 2 /2, (G 1 + G 2 )/2} with G i reciprocal lattice vectors. Its onset is determined by a critical concentration of electrons or holes signaled by a significant distortion of the single-particle spectrum, as depicted Fig. 4. When tuning towards the magic-angle the lower (upper) critical band filling decreases (increases). We however do not find any indication for a link between the observed electron instability and the van-Hove singularities.
A possible real space realization of the stripe density charge order is depicted in the inlet of Fig. 4. Note that the corresponding real-space charge distribution, accessible e.g. in STM measurements, differs qualitatively because of the highly non-local shape of Wannier orbitals with highest localization probability at the hexagon's centers. The local charge distribution would rather resemble a distorted version of the disordered state breaking rotational-C 3 symmetry.
Strong-coupling regime (II): In the magic-angle regime, kinetic energy contributions of electrons are expected to be much smaller than contribution from interaction processes (see Fig. 1). We therefore consider density-density interaction processes only. This approximation renders the local occupation number a good quantum number and the theory integrable. Later, kinetic contribution may be incorporated perturbatively in orders of ∼ t/U , which is however not part of this work. In this limit of infinite coupling, the Hamiltonian contains only contributions from direct and exchange interaction processes, wheren ασ represents the local occupation number operator. This model is particle-hole symmetric, i.e. invariant undern aσ → 1 −n aσ , which allows us to study either hole or electron doping. Since electronic single-particle states are either occupied or empty, n ασ = 0, 1, the ground state is effectively determined using the conventional Monte Carlo based simulated-annealing-algorithm [47].

Results (II):
For all commensurate moiré band fillings ν = 0, ±1/4, ±1/2, ±3/4 we find Mott-insulating ground states. When adding or removing electrons we expect that the insulator turns into a conductor where single particles move in a landscape of potential barriers generated by electrons and holes constituting the nearest Mott insulating state. Our results of possible ground states of the hole-doped side are depicted in Fig. 5.
For moiré band fillings |ν| < 1/2, we observe stripetype orders which resemble our findings of the the weakcoupling approach with charge inhomogeneities are described by ordering vectors as given in Eq. (7). It indicates, that this particular density configuration minimizes the potential energy costs generated by the dominant direct interaction processes irrespective of kinetic energy contributions. It supports our earlier presumption that the formation of density inhomogeneities in the weak coupling-regime is not linked to features of the single-particle spectrum (e.g. nesting). For |ν| ≥ 1/2, we find charge configurations which maximize the distance between charges similar to the principal of Wigner crystallization.
We conclude that the charge distribution is decisively determined by direct interaction processes U which are characterized by a significant coupling of all sites sharing one hexagon. It is determined solely by the ratio of direct interaction terms, which was determined to (U on-site : U NN :U NNN :U NNNN )/U on-site = (1:0.79:0.63:0.58), where the exact numerical values matter (i.e. simple ratios of type (1 : 2 3 : 2 3 : 1 3 ) being suggestive when counting the direct overlap of wannier orbitals lead to qualitatively different results). Furthermore, since the local singleparticle states are either empty or occupied, the particular ground state is required, unless occupied lattice sites are always fully occupied, to additionally break the spin-and / or valley-symmetry. Since direct interaction processes do not discriminate between spin and charge degrees-of-freedom, the energetically most favorable con- figuration is here determined by exchange interactions. Their matrix elements, for both the intra-and intervalley channel, are always found to be positive and therefore favor an alignment of spins non-locally (because of intra-valley exchange) and locally (because of inter-valley exchange). This results for |ν| = 1/4, 3/4 in a condensation of local degrees of freedom of partially occupied sites in one particular spin and valley sector, whereas for |ν| = 1/2 in one particular spin sector because on-site intra-valley exchange interaction is absent, as depicted in Fig. 5. These findings are in agreement with the experimental observations outlined above.
Conclusion: In this work, we employed a continuum model of TBG to obtain the characteristic moiré band structure. Subsequently, we constructed Wannier orbitals to determine the relevant interactions between localized states on the effective moiré superlattice. We found a hierarchy of interaction processes, Eq. (5), where direct interaction processes between states sharing one hexagon dominate, followed by intra-and inter-valley exchange interaction processes which are at least one order in magnitude smaller but always positive. Combined with a matrix elements' dependence on distance, those characteristic were found decisive for the analysis of possible electronic ground states. Our findings have to be contrasted to a similar strong coupling approach presented in Ref. [39] where the authors assumed an averaged interaction strength for all processes connecting states sharing one hexagon and also included processes beyond the density channel. As our results depend decisively on the distance dependence of interaction elements, our ground states for commensurate fillings |ν| = 1/2, 3/4 differ from those determined in Ref. [39].
The most robust finding of our analysis, that occurs at weak and strong coupling, is the emergence of a nematic state that breaks the three-fold rotation symmetry of the moiré lattice. The details of the related translational symmetry breaking and of additional broken symmetries depend then on the strength of the interaction and the filling fractions. While critical fluctuations, not included in our formalism, may render charge density waves, spin, or valley order finite ranged, the discrete nematic symmetry breaking should give rise to a sharp second-order phase transition with a finite transition temperature. Even for moderate symmetry-breaking substrate-induced strain we expect a well defined crossover temperature. In addition to the nematic state, we find an onset of spinand valley-polarized orders at strong coupling. This effect is due to the Hund's coupling of the non-local and positive intra-valley and inter-valley exchange couplings, which maximizes first the spin and subsequently the valley number when filling up superlattice sites with electrons. This yields the observed degeneracy pattern of the observed insulating states.
If the nematic order exist away from incommensurate filling the reduced symmetry at T c exclude more complex superconducting order parameters, such as chiral d+id or nematic d + id states. On the other hand, the abundance of nematic order in twisted bilayer suggests that nematic fluctuations may be important in inducing or amplifying superconductivity in these materials [29,48,49]. is minimized. Here, r iα represents the coordinates of the Wannier orbital's center located at the center of the ABor BA-stacked regions of the ith superlattice unit cell. Exemplarily, the obtained Wannier orbitals for a twist-angle of θ = 1.05°, which are checked to be exponentially localized, are depicted in Fig. 6. Eventually, the transition amplitudes of the effective superlattice tight-binding model, introduced in Eq. λk . (A5) The single-particle spectrum using transition amplitudes {t iα,jβ } is depicted in Fig. 7.