1 Introduction

Recently it has been showed that twisted bilayer graphene systems support superconducting phases at certain special twist angles ^{[1, 2]}, where electronic correlations are maximized due to the existence of flat-bands ^[3]. Moreover, Cao et al. ^[1] found that such TBG systems have a Mott insulating phase that appears in the middle of unconventional superconducting phases, similar to the phase diagram found in cuprates and other high temperature superconductors ^[4]. The appearance of a Mott correlated insulator and unconventional superconducting phases in the flat band of magic-angle TBG at a small carrier density cannot be explained by weak-coupling BCS theory ^{[5- 14]}. In fact, although T _c is very small, of the order of 1.7K, the electron-electron coupling turns out to be very high. In Fig. 1 we show the variation of In(T _c / T _F ), where T _F is the Fermi temperature and T _c is the critical temperature of the superconducting phase as a function of the charge doping n’. This plot intends to compare the electron-electron coupling once the density of carriers and T _F are taken into account. The magic-angle TBG T _c / T _F ratio is above the trend lines on which most heavy fermions, cuprates and organic superconductors lie ^[1]. Thereafter, it has been experimentally found that trilayer graphene rotated by magic angles turns out to be the highest ever found electron-electron correlated material ^[15]. Up to now, the origin of superconductivity in TBG remains under debate.

Figure 1 Variation of In(T _c / T _F ) (logarithmic scale) as function of n’ (charge doping) in scaled units of 10¹¹ cm^-2 for magic-angle TBG (blue points). Here T _c is the critical temperature for a superconductor state and T _F is the Fermi temperature. The horizontal lines indicate the approximate T _c / T _F values of the corresponding family of materials. Adapted from Ref. ^[1].  

An interesting strong coupling theory of superconductivity based on skyrmions could explain the mechanism of superconductivity in TBG ^[11]. For this purpose is important to understand the quantum geometry of flat bands in magic-angle TBG which is purely the mathematical structure necessary to measure distances at the quantum regime. Also is important to analyze the interplay between interaction range and Berry curvature inhomogeneity. A deeper theoretical study by Ledwith et al. ^[16] based in vortexability and the importance of chiral symmetry to induce vortexable bands as a generalization of the LLL (Lowest Landau Level) for hosting a short-range interacting ground state (SRI-GS) is of great relevance. Vortexable bands are a class of bands to which one can attach vortices while remaining within the subspace defined by the bands. This treatment of vortexable bands allows the construction of exact many-body FQHE (Fractional Quantum Hall Effect) ground states in the limit of short-range interactions, and possibly vortexable bands of equal and opposite Chern numbers related to superconductors based on skyrmions. On the other hand, magic-angle TBG has flat Chern bands at zero magnetic field. Therefore, TBG promises a route towards stabilizing zero-fields Fractional Chern Insulators (FCIs) ^[17].

In fact, the study of the electronic properties of twisted bilayer graphene started before the discovery of superconductivity at magic angles. In works of 2007 by J. Santos ^[18] and 2011 by A. Macdonald ^[19], an effective low energy continuous Hamiltonian model was derived.

In such studies the idea was to generate moiré patterns as a function of the twist angle between graphene layers. In Fig. 2 we show the unit cell for the moiré pattern of TBG. The new unit cell vectors in the moiré lattice are scaled as the inverse of the twist angle between layers, therefore, for small twist angles, it is expected to have a bigger moiré unit cell. Form there is possible to define a moiré Brillouin zone (mBZ). For small twist angles, the mBZ is also small due to the moiré modulation vector. Flat-band were found at certain angles and from there superconductivity was predicted ^[19].

Figure 2 Picture of twisted bilayer graphene (TBG) system. The graphene monolayer 1 (green) and monolayer 2 (purple) have a relative twist angle that produce a moiré pattern (cyan) with a unit cell much bigger than the original unit cell of each monolayer graphene. The unit cell of the moiré pattern is indicated (red).  

The reason for the appearance of flat bands and its quantized nature is not understood. Therefore, lots of research have been conducted in this direction. Notice that flat bands can also occur in graphene over substrates even without twists yet TBG systems are paradigmatic. To make the TBG model more realistic, G. Tarnopolsky et al. ^[20] took into account the structural relaxation due to carbon-carbon repulsion between layers and a chiral continuum model was produced. Perhaps, so far, it is the simplest and more realistic model that best captures the nature of magic angles; at these special angles the dispersion energy at the lowest bands becomes flat and has a recurrence behavior. Also at these magic angles the Fermi velocity goes to zero. Due to its chiral symmetry, the Hamiltonian of this model produces an intra-valley inversion symmetry ^[21] so the energy dispersion is inversion symmetric at all twist angles and thus symmetry protected at any twist angle. The zero-mode have some resemblance to the ground state of a quantum Hall effect wave function on a torus ^{[20, 22]}, and, therefore, the solution is of the quantum harmonic oscillator type, where Landau levels arise ^{[23, 24, 25]}. Another interesting mathematical characteristic is that flat-band modes are constructed from a complex analytic functions that are ratios of the theta Jacobi functions ^{[20, 22]}.

Furthermore, due to the quantized nature of the magic angles it appears that there is an adiabatic change in integers of this topological invariant. However, there are many open question concerning the possible relationship between the quantum hall effect (QHE) and the TBG, the quantized nature of flat-bands and the nature of wave functions.

In a series of previous papers, the authors showed that by making a supersymmetric transformation akin to take the square of the original chiral TBG hamiltonian, the physics behind the magic angles becomes clear ^{[26, 27, 28]}. The aim of this work is to review the topic and clarify the physical behavior of this model as well as to develop an effective equation for such squared Hamiltonian in order to study all other non flat band states. Although one of the authors presented before in this journal some details of the squared Hamiltonian model ^[29], here we go further by showing that the effective system is related to pseudo-magnetic fields, bringing the possibility of writing the squared chiral TBG Hamiltonian in a non-Abelian fashion where effective magnetic fields appear, making the connection with the QHE transparent. Then we discuss the physical picture that arises from this minimal model including an analytical calculation of the first magic angle to show how this arises as a balance between kinetic, confinement energies and interlayer currents.

Let us finally add that the original 4 × 4 Hamiltonian is an effective model that contains hundreds of atoms in its unitary cell which is of the size of nanometers (see Fig. 2). However, it is made from a basis with bi-spinors and thus can be thought as representing the dynamics of two coupled quasiparticles. Therefore, one can still interpret the original 4 × 4 as two quasiparticles in a non Abelian field. However, this is much more difficult to visualize than in the approach presented here, where one of the quasiparticles is renormalized by squaring the Hamiltonian.

2 Chiral TBG system

The Bistritzer-MacDonald Hamiltonian (BMH) model arises by a simple stacking of two rotated graphene layers with an interaction between them ^[19]. One of the authors gave a simple derivation of this model in this journal before and thus we refer the reader to it ^[29]. Although the BMH contains the magic angle physics, in real TBG stacking points where carbon atoms are on top of another carbon (known as AA stacking) are pulled apart out of the intralayer plane due to Coulomb repulsion. As a consequence, interlayer electron tunnelling in AA sectors is greatly reduced and can be safely neglected. This result in a chiral Hamiltonian (CH) ^{[13, 20]}. Let us write this model. Consider as basis the wave vectors Φ(r)=ψ1(r),ψ2(r),χ1(r),χ2(r)T where the index 1,2 represents each graphene layer and ψj(r) and χj(r) are the Wannier orbitals on each site of the graphene’s unit cell, i.e., in sites A and B of the graphene’s bipartite lattice. Using this base, the CH is ^{[20, 30, 31]},

where we defined the zero mode operator as,

and its rotated and conjugated version,

We also defined the antiholonomic operator,

and the holonomic operator,

Both operators eventually play an important role in the theory, as the flat-band wavefunctions are in part determined by the fact that the antiholonomic derivative of any analytic function is zero ^[28]. Although these operators appear in graphene, the lack of a position dependent potential in graphene allows to solve the problem in a very simple way without the need to look for the analytical properties of solutions. The interlayer potential is ^[20],

Here ϕ=2π/3 and q1=kθ(0,-1), q2=kθ(3/2,1/2), q3=kθ(-3/2,1/2).

The distance between each layer graphene Dirac cone is kθ=2kDsinθ/2 with kD=4π/3a0 the Dirac wave vector and a0 the lattice constant of graphene. The physics of this model is captured via the parameter α, defined as,

where w1 is the interlayer coupling of stacking AB/BA with value w1=110 meV and v0 is the Fermi velocity with value v0=19.81  eV/2kD.

Here we will use units where v0=1 and re-scaled vectors such that qμ'=qμ/kθ and r'=kθr, therefore, r⋅qμ=r'⋅qμ' is independent of kθ thus the antiholonomic and holonomic operators are dimensionless. In other words, we can consider the system as if it were with a fixed geometry, i.e., setting q1':=q1=(0,-1), q2':=q2=(3/2,1/2), q3':=q3=(-3/2,1/2) while the twist angle enters only in the coupling parameter α. The vectors b1,2=q2,3-q1 turn out to be the moiré reciprocal vectors and for its utility, we define a third vector b3=q3-q2.

To finish the model, we found useful to define a set of unitary vectors q^μ⊥ given by q^1⊥=(1,0), q^2⊥=(-1/2, 3/2), q^3⊥=(-1/2,-3/2).

3 Squared TBG Hamiltonian and Pseudo-magnetic fields

In a previous work ^[26] we showed how, by taking the square of H, which is akin to consider a supersymmetric model, it is possible to write the Hamiltonian as a 2 × 2 matrix. This transformation takes into account the particle-hole symmetry and folds the band around E = 0. Physically, it removes one of the bipartite lattices on each layer, leaving two triangular lattices on each layer. In a series of previous works we investigated several properties of this Hamiltonian, among the most important is the mapping of the flat-band into a ground state separated by a gap form the rest of the states. This state has an antibonding nature in a triangular lattice and then has frustration, associated with a massive degeneration ^[32]. The connection with supersymmetry is important as there are many recent articles exploring several properties such as topology and supersymmetry on lattice systems and which are of great relevance, thus we refer the reader to them ^{[33, 34, 35, 36]}. In this context is specially intriguing the relationship found before between the squared TBG Hamiltonian and a boson Hamiltonian describing phonons in a flexible system^[26].

For the present work we decided to change the notation used in such works and write H ² in terms of a more physically suggestive picture by using pseudo-magnetic fields,

where U± denotes the intralayer confinement potentials ^[26],

and,

which turns out to be a two-dimensional pseudo-magnetic vector potential in the Coulomb gauge as is easy to prove that ∇⋅A±=0. We also define a pseudo-magnetic scalar potential,

which has the property,

where the cross product between two-dimensional vectors given by A=(Ax,Ay) and B=(Bx,By) is A×B=AxBy-AyBx. Observe that the related pseudo-magnetic flux over the real-space unit cell vanishes,

Such cancelling can be graphically understood in Figs. 3a)-b) where we plot the pseudo-magnetic vector potential in real space. Notice that B has an imaginary part which at first sight seems to be an odd fact. However, after the first experimental observation for such effective fields in the evolution of quantum coherence of a spin coupled to an Ising-type spin bath ^[37], there are lots of works dealing with the subject ^{[38, 39]}. Here it can be interpreted as in Ref. ^[37], i.e., the quantum space evolution introduces a complex phase factor inside the field. On the other hand, in Figs. 3c)-d) is shown the pseudo-magnetic scalar field. In both cases, we have a Kagome geometry for the pseudo-magnetic fields of the squared hamiltonian.

Figure 3 In panels a) and b), we show the real and imaginary parts of the pseudo-magnetic vector potential A∓(r), respectively. In panels c) and d), we show the real and imaginary parts of the pseudo-magnetic scalar potential B∓(r), respectively. The arrows denote the direction of the vectors while the color is used for the magnitude. In all panels we indicate the stacking points AA (red), BA (yellow) and AB (green) over the real-space unit cell (black hexagon) and the cyan hexagon represents the magnetic unit cell. The magenta dots are related to the centers of charge for α→∞ in the electronic density, and are also related to the QHE ^[28]. 

4 Non-Abelian model for TBG

Let us now explore the nature of Eq. (8). In analogy to an electron in a pseudo-magnetic vector potential, we define the canonical momentum,

where Π~ is the canonical SU(2) momentum, e the electron charge and A^ is a non-Abelian SU(2) (see below) pseudo-magnetic vector potential,

with A^x=A1,xτ^1+A2,xτ^2 and A^y=A1,yτ^1+A2,yτ^2, where we used the set of Pauli matrices τ^j (with j=1,2,3) in the pseudo-spin of the layer degree, and the 2 × 2 identity matrix τ^0. Explicitly, the components of A^ are,

Note that A^ is non-Abelian as follows from the fact that A^μ(r),A^ν(r')≠0, i.e., μ≠ν.

Using the canonical momentum with natural units ℏ,e,m,c=1, the squared Hamiltonian can be written as,

Eq. (17) can be further simplified to,

where F^μν=∂μA^ν-∂νA^μ-iαA^μ,A^ν is the Zeeman coupling. Therefore, Eq. (18) proves that indeed the squared chiral Hamiltonian describes an electron in a non-Abelian SU(2)pseudo-magnetic field as was suggested before ^[40].

Now we investigate the Zeeman term. We start by observing that,

from where,

Substituting Eq. (20) in the commutator Eq. (19), it follows that for r=r' ,

where we specified that the commutator must be taken at the same point. We also used that 2iτ^3=τ^1,τ^2. From Eq.(21) it follows that,

multiplying by α2 Eq. (22) and substituting in Eq. (21),

where we defined the potential of the non-Abelian term h3(r)=3α2∑μ‍(-1)μ-1sin(bμ⋅r)=α2/2(U--U+).

In Fig. 4 we show the contour plot in real space for the potential h3(r). As in Fig. 3, it can be seen that the pseudo-magnetic fields of the squared Hamiltonian behave like a Kagome lattice. Moreover, note that these pseudo-magnetic fields are described in a larger lattice due to the renormalization process of the graphene lattice into a triangular lattice. Therefore, the original real-space unit cell is now three times larger where the pseudo-magnetic fields operate.

Figure 4 Representation in real space of a) A2(r) and b) h3(r). Where is indicated the stacking points AA (red), BA (yellow) and AB (green) over the real-space unit cell (black hexagon) and the cyan hexagon represents the magnetic unit cell. The magenta dots are related to the centers of charge for α→∞ in the electronic density, and are also related to the QHE ^[28]. 

The other term of the Zeeman coupling is related to the pseudo-magnetic scalar potentials as follows,

but qμ×q^μ⊥=q^μ⊥yqμx-q^μ⊥xqμy=1, therefore Eq. (24),

where B±(r) is the pseudo-magnetic interlayer potential and the non-Abelian pseudo-magnetic vector potential is written in matrix form as,

In the same way, the square of the non-Abelian pseudo-magnetic vector potential A^ is related to a new potential V(r)=α2A2 where,

It follows that,

and therefore α2A^2=V(r)τ^0. Figure 4 presents A ² which in fact is a confinement potential akin to a quantum dot.

Using all these results, it is possible to write the square Hamiltonian Eq. (8) in a way that manifest more clearly its non-Abelian nature. We start by writing H ² purely in terms of A,

where Ax,Ay denotes the commutator of operators Ax and Ay. Also we can write H ² as follows,

where we defined h3(r)=-α2iAx(r),Ay(r) and V(r)=α2A2. Also we define the off-diagonal components as h1(r)=-∑μ‍sin(qμ⋅r)+2icos(q⋅r)q^μ⊥⋅∇ and h2(r)=∑μ‍cos(qμ⋅r)-2isin(q⋅r)q^μ⊥⋅∇. This last expression is very important as it manifest in a clear fashion the topological nature of the hamiltonian as it includes all three Pauli matrices ^[26].

As proved in Ref. ^[29], in the limit α→∞ we have that ∂μA^ν-∂νA^μ≫-iαA^μ,A^ν. Thus the Zeeman coupling is F^μν≈∂μA^ν-∂νA^μ. Therefore h3(r) is neglected and the system has a pure Abelian field for higher magic angles ^[28]. The pure Abelian field means that a constant effective magnetic can be found and eventually, it leads to the proof that in such limit, the squared TBG hamiltonian at magic angles and the Quantum Hall effect hamiltonian are the same ^[28]. The other limit α→0 can be explored using perturbation theory, as done in the following section.

5 Strong non-Abelian limit: first magic angle perturbative solution at Γ-point

In this section we will explore the strong non-Abelian limit of the hamiltonian and its relationship with the first magic angle. Before doing so, and as explained elsewhere ^[29], the Γ point can be used to reveal the magic angles as E2(Γ) is always the highest energy for the central band. This is seen in Fig. 5, where we plot E2(Γ) as obtained from the numerical calculation detailed in the Appendix. For a flat band to exist, the energy E2(Γ) must be zero.

Figure 5 

The expected values for a given layer j = 1,2 are, the kinetic energy,

the confinement energy,

and the energy contribution from the off-diagonal operators,

here A^j and B^j are the projections of the non-Abelian pseudo-magnetic vector potential and scalar potential, respectively, on each layer j of the pseudo-spin degree. Therefore, the total energy for the squared Hamiltonian is,

where the index j = 1,2 for take into account the two layers.

For the Γ point it is very illustrative to use perturbation theory in the limit α→0. The corresponding wave function is,

and ψΓ,2(r)=iμαψΓ,1(-r), where Rγr is a counterclockwise rotation on angle γ with tan(γ)=3/5 and μα=±1, the minus sign is used for odd magic angles. In the Γ point at a given layer we obtain up to second order in α that,

It follows that,

In Fig. 6 we compare these perturbative expected values with those obtained from the numerical results (see the appendix) showing a very good agreement for α→0. Our perturbative analysis allows to get a glimpse of the first magic angle position as using Eq. (39) the condition for having a flat band turns out to be,

Figure 6 Expected values in the Γ point vs α. The numerical results are indicated with dashed lines with markers. The kinetic Energy 〈Γ|T^|Γ〉 (purple), confinement energy 〈Γ|A^2|Γ〉 (cyan), 〈Γ|A^⋅∇|Γ〉 (magenta) and 〈Γ|B^|Γ〉 (brown). The solid lines are the perturbative solutions. The vertical line indicates the first magic angle value α1=0.586. Notice how 〈Γ|T^|Γ〉=〈Γ|A^2|Γ〉 at α1. 

Therefore, α1≈1/2, a value which is close to the numerically found first magic angle which is at α1=0.586. Higher order terms in the expansion are needed to increase the accuracy, but yet the main principle behind a magic angle is already present in this simple approach.

Then we conclude that in going from α=0 to α1, the confinement potential starts to contribute and reaches the kinetic energy at the magic angle. The off-diagonal operators always diminish the energy. As expected, the first magic angle is thus produced when the sum of the kinetic plus confinement energies are equal in magnitude to the expected values of the off-diagonal operators.

6 Weak Non-Abelian limit and higher order magic angles: equivalence with the Quantum Hall effect

As we discuss below, for higher-order magic angles the non-Abelian effects are weak since the system can be treated as if the magnetic field is constant due to a strong confinement. In fact, we recently proved that the H ² hamiltonian converges into the Quantum Hall effect hamiltonian, i.e., it is equivalent to a two dimensional quantum oscillator ^[28]. We refer the reader to this reference for more details. Here we only make some remarks concerning the strong confinement. To show how it arises, consider the zero-mode equation D(r)ψ1(r),ψ2(r)T=0. Although all k points are suitable to be studied, it is easier to obtain the Γ point solution as due to symmetry, ψ2(r)=iμαψ1(-r) where μα=±1 depending on the magic angle parity ^[20]. The zero-mode equation is thus translated into the following set of equations,

In the limit α→∞, boundary layer theory of differential equations ^[42] applied to Eqs. (41)-(42) indicates that if gradients are small, the derivatives can be neglected and thus the solution must satisfy ψj(r)→0. Only in the boundary layer, i.e., when ∂-ψj(r) is of order αU(±r)ψj(r) a solution different from zero is found. We conclude that zero-modes solutions are strongly peaked around certain regions of space indicating a strong confinement.

From there it is possible to prove that flat-band states converge into coherent Landau levels ^[28]. These are the smallest wave-packets that can be constructed from states within the lowest Landau level ^[28]. The resultant electron density is made from a linear combination of Gaussians ^[43] with mean deviation given by σ≈1/3α. As σ→0 when α→∞, the non-Abelian nature of the field is effectively decreased as the electron will not probe all regions of space. As the problem is thus similar to a constant magnetic field, from there it is obtained that the squared Hamiltonian becomes the Quantum Hall effect Hamiltonian ^[28].

7 Conclusion

The chiral TBG Hamiltonian is equivalent to an electron inside of a non-Abelian pseudo-magnetic field. Here we gave its corresponding general expression. Then we explored the first magic angle using perturbation theory to show how the kinetic, confinement energy and interlayer currents produce a magic angle. For the first magic angle results that the non-Abelian properties are more important than other higher magic angles because the commutator becomes less important. Therefore, in the limit α→∞, we recover an effective Abelian behavior that explains in part the effective equivalence between the Quantum Hall effect hamiltonian and the TBG hamiltonian ^[28]. The emergent Kagome structure observed in the pseudo-magnetic fields and for small angles α→∞ ^{[28, 44]}, opens the question of whether there is any mechanism related to the Kagome lattice of the chiral magic angle TBG like spin-liquids or magnetic frustration, and what role it plays in the mechanism of superconductivity, FQHE and the physics of such strong correlations in TBG. This allows us to take advantage of this ideal scenario to build highly entangled ground states and opens many exciting technological possibilities in moiré materials.