Short Versus Long Range Exchange Interactions in Twisted Bilayer Graphene

This study discusses the effect of long‐range interactions within the self‐consistent Hartree‐Fock (HF) approximation in comparison to short‐range atomic Hubbard interactions on the band structure of twisted bilayer graphene (TBG) at charge neutrality for various twist angles. Starting from atomistic calculations, it determines the quasi‐particle band structure of TBG with Hubbard interactions for three magnetic orderings: modulated anti‐ferromagnetic (MAFM), (NAFM) and hexagonal anti‐ferromagnetic (HAFM). Then, it develops an approach to incorporate these magnetic orderings along with the HF potential in the continuum approximation. Away from the magic angle, it observes a drastic effect of the magnetic order on the band structure of TBG compared to the influence of the HF potential. Near the magic angle, the HF potential plays a major role in the band structure, with HAFM and MAFM being secondary effects, but NAFM appears to still significantly distort the electronic structure at the magic angle. These findings suggest that the spin‐valley degenerate broken symmetry state often found in HF calculations of charge neutral TBG near the magic angle should favor magnetic order, since the atomistic Hubbard interaction will break this symmetry in favor of spin polarization.

A significant limiting factor of self-consistent atomistic approaches for broken symmetry phases is their computational cost [69,72,73].Some examples exist in the literature [69][70][71]74], but a full phase diagram -over twist angle and doping level, amongst other experimental variables -has not yet been achieved.For example, Stauber and González [70,71] developed a theory based on Green's functions, where only some of the states were retained.They were able to study long-ranged interactions, and also the interplay between long and short ranged interactions, but only at 1.16 • and either at charge neutrality or −2 electrons per moiré unit cell.Moreover, Vahedi et al. [69] investigated several twist angles (1.08 • , 1.30 • and 1.47 • ) but only focused on charge neutrality.Usually, this has been overcome from re-scaling the TB parameters [69,[75][76][77] or applying hydrostatic pressure [78][79][80], such that flat bands can be created with unit cells only containing a few hundred to a few thousand atoms.
Here we develop an approach which can include shortranged interactions, such as the on-site Hubbard interaction of the p z orbitals, in the continuum model.Starting from the RPA spin susceptibility calculations of Klebl and Honerkamp [63], we perform self-consistent atomistic Hubbard calculations to obtain the mean-field magnetic order parameters for different ordering tendencies (at a large twist angle and charge neutrality).We develop analytical forms for the magnetic ordering in real space, which we are able to include in the continuum model as an effective a scalar sublattice potential.In order to elucidate the angle dependence of the interplay between the magnetic orderings and the Hartree-Fock potential, we perform self-consistent Hartree-Fock calculations at charge neutrality, to which we later add the effective magnetic potential at different twist angles.Overall, it is found that the long range contribution dominates at the magic angle, but away from the magic-angle, these magnetic orders become more significant.We discuss the competition between these long and short range exchange interactions in detail, and finish with a discussion of future directions.

A. Atomistic Calculations
We study commensurate moiré unit cells of TBG [36], starting from AA stacked bilayers and rotating the top layer anticlockwise about an axis perpendicular to the graphene sheets that passes through a carbon atom in each layer.The moiré lattice vectors of the commensurate structures are R 1 = na 1 +ma 2 and R 2 = −ma 1 +(n+ m)a 2 , where n and m are integers which define the commensurate TBG structure, and a 1 = ( √ 3/2, −1/2)a 0 and a 2 = ( √ 3/2, 1/2)a 0 are the lattice vectors of graphene, where a 0 = 2.46 Å is the lattice constant of graphene.
To investigate the electronic structure of TBG, we use the Hamiltonian where ε iσ and ĉ † iσ (ĉ iσ ) denote the on-site energy of atom i with spin σ and the electron creation (annihilation) operator associated with atom i and spin σ, respectively.The hopping parameters between atoms i and j, t(r i − r j ), are calculated using the Slater-Koster rules [89] where We take the pre-factor for the ppσ-hopping and ppπ-hopping to be V 0 ppσ = 0.48 eV and V 0 ppπ = −2.7 eV, respectively.The carbon-carbon bond length is given by a = a 0 / √ 3 and the interlayer separation is taken to be d AB = 3.35 Å.We take the decay parameters of the Slater-Koster rules to be q σ = d AB /0.184a 0 and q π = 1/0.184√ 3 [36,37].Hoppings between carbon atoms separated by more than R c = 10 Å are neglected [90].
To include the effects of short-range Hubbard interactions, the on-site energy is determined by the mean-field Hubbard interaction where U is the Hubbard parameter of the carbon p z orbital, and n iσ is the mean-field electron density on atom i with the spin σ being the opposite to σ.
The electron density can be determined from the Bloch eigenstates ψ nkσ (r) (with subscripts n and k denoting a band index and the crystal momentum, respectively) according to where f nkσ = Θ(ε F − ε nkσ ) is the occupancy of state ψ nkσ with eigenvalue ε nkσ (where ε F is the Fermi energy), χ j (r) = R φ 2 z (r − t j − R) (with R denoting the moiré lattice vectors) and n jσ is the total number of electrons in the j-th orbital with spin σ.To characterizes the magnetic ordering, we calculate the spin polarisation To start the mean-field calculations, instead of guessing random magnetic configurations, we perform RPA spin-susceptibility calculations, following the methods outlined in Refs.63, 64.The eigenvalues of these calculations provide the critical interaction strength of an instability (U c ) and the eignevector is the form of the magnetic order.By using these eigenvectors as an initial on-site interaction, we can induce spin polarisation, which can then be used to perform self-consistent calculations.
To obtain a self-consistent solution of the atomistic Hubbard equation, we use a simple mixing scheme with a mixing parameter of 0.1 typically (0.1 of the new electron spin density is mixed into the same spin density).When determining the Fermi energy, the total electron number is again forced to be N + ν, but this does not restrict the spin densities to be the same.We mix the up and down spin density by the same amount, instead of choosing to work with the total electron density and magnetic order parameter, as we find it is sometimes more stable.
Only the leading instability was able to be stabilised with this method [91].Therefore, we also perform constrained calculations.The constrained calculations require an analytical form for the magnetic order to be specified (which is outlined in Section III A), say ζ (j) , and we self-consistently determine the magnitude of this magnetic order through projecting onto the magnetic order parameter, ζ.The on-site energy term is then given by where the sign depends on the spin.A self-consistent solution is obtained through a linear mixing of the spin densities.

B. Continuum Model Calculations
The mini-Brillouin Zone (mBZ) of the continuum model is spanned by the two reciprocal lattice vectors given by, is the moiré period and a 0 is the lattice constant of graphene.These vectors form the basis to define any reciprocal lattice vector, The non-interacting continuum model is 4-fold degenerate, since it accounts for valley and spin quantum numbers, with the Hamiltonian at crystal momentum k being written as where Ĥl,ξ is the continuum single layer graphene Hamiltonian of valley ξ, given by with v F = ( √ 3V 0 ppπ a)/(2 ) denoting the Fermi velocity, K l is the position of the Dirac point of layer l, τ θ,l = e iξτzθ/2 τ x , ξτ y e −iξτzθ/2 , with τ i being the Pauli matrices acting on the sublattice degree of freedom.The matrix T is a periodic function in the moiré unit cell that hybridises layers.For small angles, the main contribution comes from the first three reciprocal lattice vectors, where u 1 = 0.0797 eV and u 2 = 0.0975 eV [92] are, respectively, the hopping amplitudes between AB/BA and AA stacking, which takes into account the atomic relaxation in the continuum model.
To account for electron-electron interactions, we include the mean-field Hartree-Fock terms to the Hamiltonian.The Hartree contribution to the Hamiltonian is given by where i ∈ [1,4] labels the combined sublattice and layer degree of freedom, σ accounts for the spin, Ω is the area of the mBZ, and the local Hartree potential is given by Here δρ (r) ≡ ρ (r) − ρ CN (r) denotes the fluctuation in charge density, with ρ (r) = ξ,σ ψ † ξ,σ (r) ψ ξ,σ (r) corresponding to the charge density and ρ CN (r) is the average density corresponding to the non-interacting TBG at charge neutrality point.We assume that the Coulomb interaction is screened by a double-metallic gate [45] where d = 40 nm is the distance to the metallic gates and ε = 10 is the dielectric constant [49,58,67].
The Fock contribution to the Hamiltonian is given by where i, j run over the sublattice and layer indexes.The non-local Fock potential is described by, As we want to express the matrix elements of ĤF , defined in Eq. ( 13), in the reciprocal space for which we must transform the non-local Fock potential into the Fourier space.By this procedure we compute the Fock matrix elements as where the index n runs over the occupied bands at a given Fermi energy.For the Hartree-Fock calculations we work with a continuum model of TBG expanded up to the third star.We use a density of points between 2 − 6 × 10 5 Å2 in the mBZ, depending on the twisting angle.The convergence of the Hartree-Fock potential is normally reached after 5 − 6 self-consistency steps and it has been proved for increasing density of points.
To include the α-magnetic potential (α = M, N, H) in the continuum model we use a scalar sublattice and spin dependent potential expressed through its harmonic decomposition in the first star reciprocal lattice vectors The full form of U α δ (G i − G j ) is discussed in the Section III C. The final Hamiltonian that combines both the effective magnetic potential derived from atomistic calculations and the Hartree-Fock potential at half-filling is given by, Note that this final Hamiltonian is not treated in a selfconsistent way since the effective magnetic Hamiltonian is just added to the self-consistent Hartree-Fock Hamiltonian, which would be equivalent to a first-order approximation of the magnetic orderings in perturbation theory.

A. Short Range Atomistic Hubbard Interactions
From the RPA (q = 0) spin-susceptibility calculations [63,64], a number of leading antiferromagnetic instabilities are found: modulated antiferrogmagnetic order (MAFM), nodal antiferromagnetic order (NAFM) and hexagonal antiferromagnetic order (HAFM).To find which instability is the ground state at a given twist angle and doping level, (mean-field) atomistic Hubbard calculations must be performed.As these atomistic calculations are extremely computationally expensive at the magic angle, we focus on 1.54 • at charge neutrality.At this twist angle and doping level, TBGs leading instability was found to be MAFM (with a critical Hubbard interaction of U c ≈ 5.1 eV), with NAFM and HAFM having slightly larger critical interaction strengths (of U c ≈ 5.4 eV).For these latter instabilities, we perform unconstrained and constrained atomistic Hubbard calculations, as outlined in Section II A.
The MAFM instability, as seen in Fig. 1(b) and (c), is characterised by a sub-lattice oscillation in the magnetic order parameter ζ = (n ↑ − n ↓ )/(n ↑ + n ↓ ) that is modulated throughout the moiré supercell [63].In Fig. 1(b), which plots the magnetic structure along the diagonal of the moiré supperlattice, sublattice A is shown in black and sublattice B is shown in grey for the top graphene layer.To show its real-space structure more clearly, Fig. 1(c) plots the magnetic order parameter on sublattice B of the top layer over the moiré superlattice.The constant contribution of the MAFM order is larger than the moiré-scale variation, which means the sub-lattice polarisation is the same throughout the moiré unit cell.The moiré scale variation of the magnetic order peaks in the AA regions, as might be expected from the LDOS of TBG peaking peaked in the AA regions [36].The MAFM order can be approximated with the following analytical form where Therefore, the magnetic structures from these atomistic calculations are involving many more electrons than just the 4 flat band electrons, which is in agreement with other works [69][70][71].
In Fig. 1(a), we show the corresponding self-consistent quasi-particle band structure.The different valleys, K and K', have been identified by applying the valley operator to the states (see Refs. [76,77,93] for details of this calculation), and shown in solid black and dotted grey, respectively.Since the MAFM order breaks C 2 symmetry, it causes a gap to open at the Dirac cones at the K/K' points of the moiré Brillouin zone of TBG.This instability was not able to be stabilised at U = 5.1 eV, but for U 's larger than 5.4 eV, the constant contribution dominates and it becomes graphene-like (ζ s ζ s , such that there is only a constant sub-lattice polarisation in each graphene sheet), with the gap at the K/K' points becoming very large (100's of meV) [91].
Similarly, NAFM also has a moiré-scale peak in the magnetic order in the AA regions, but it does not possess a constant contribution to ζ, as seen in the self-consistent  (c,f,i) Corresponding plots in real space for a single layer and sublattice, where only sublattice B of the top layer is shown.Note, the U 's were chosen to be slightly above the critical values.For MAFM, if U = 5.94 eV is used, the gap at the Dirac point is extremely large [91].
values plotted in Fig. 1(e).The corresponding real-space structure is shown in Fig. 1(f), where nodes in the magnetic order around the AA region separate the regions of opposite signs of ζ spin polarisation.This magnetic order is referred to as nodal anti-ferromagnetic order because ζ goes through 0 between the AA and AB/BA regions, causing the sign of ζ to change on each sub-lattice between these types of stacking [63].Therefore, it can be described by This instability was not the leading instability [63], and we found the unconstrained calculations could never sta-bilise this order, as it would always eventually revert to MAFM order.Therefore, we performed constrained mean-field Hubbard calculations to find the "excited" magnetic order, as explained in Section II A. The quasiparticle band structure is shown in Fig. 1(d), and again it is a Mott insulator.A large gap is obtained in the Dirac cone because of the Hubbard interaction parameter being taken as U = 5.94 eV.This instability could be stabilised with U = 5.67 eV and also larger U 's in the constrained method.
Finally, the HAFM instability is similar to NAFM, but where the magnetic order parameter peaks on the AB/BA regions instead of the AA regions.The self-consistent HAFM magnetic structure is shown in Fig. 1(h) and (i).The sign of ζ changes between the AB and BA regions of the moiré unit cell.As the peaks of ζ occur on the AB/BA regions, which form a hexagonal lattice on the moiré scale, this ordering is referred to as hexagonal anti-ferromagnetic order.This magnetic order is shown in Fig. 1(h) and (i), and can be described by the analytical form This instability never appears to be a leading instability, but its critical Hubbard interaction was ∼ 5.4 eV, which is only slightly higher than the leading instability.Therefore, we again found that constrained calculations were required to obtain mean-field values of its order parameter, as explained in Section II A. In Fig. 1(g) we show the mean-field quasi-particle band structure for U = 5.94 eV (this was the smallest U which could stabilise the order with constrained calculations), where the different valleys have been coloured solid black and dotted grey.We find that this magnetic order does not create a gap at the K/K' point, despite it having a sublattice oscillation.This is because of the moiré-scale sine nature of its variation, which means it does not break C 2 on the moiré scale.It does cause the valleys to split at the K/K' points, however, resulting in one valley being pushed higher in energy and the other to lower energies.This is analogous to the effect of a perpendicular electric field has on the electronic structure.
These calculations give some insight into the magnetic order from Hubbard interactions in TBG.However, we performed these well away from the magic angle and for Hubbard interaction parameters that are large (U ≈ 5.5 eV for these calculations, but the value is thought to be ∼ 4 eV [75]).Therefore, to understand the role of these instabilities close to the magic angle, we aim to include these magnetic states in the continuum model.

B. Long Range Interactions in the Continuum Model
In Fig. 2 we show the Hartree-Fock quasi-particle band structure (solid blue line, for all subplots) for a number of twist angles at charge neutrality, in addition to the non-interacting band structure (solid red line, for all subplots).At 1.54 • , Figs. 2(a), (d) and (g), we find that the Hartree-Fock potential slightly modifies the non-interacting band structure, indicating that this twist angle is too large for the formation of an insulating state.At a twist angle of 1.25 • , Figs. 2(b), (e) and (h), the non-interacting bands are flat enough for the onset of a small gap due to the Hartree-Fock potential, significantly smaller than the bandwidth, at the Dirac cones K and K' points in the moiré Brillouin zone.Right at the magic angle of 1.05 • , Figs. 2(c), (f) and (i), the Hartree-Fock interactions induce a large gap at the K and K' points, on a similar scale to the bandwidth.At half filling (charge neutrality) these insulating states are characterised by a broken sublattice symmetry and a preserved spin and valley symmetry with respect to the non-interacting picture.This implies that the mean-field ground state associated to the Hartree-Fock band structure is actually a linear combination of 4 states with different spin and valley indexes.These calculations are in good agreement with a large body of literature which investigates these long ranged interactions ...

C. Short Range Interactions in the Continuum Model
In Section III A, we described the leading antiferromagnetic instabilities obtained from RPA calculations of the atomistic model, performed mean-field Hubbard calculations of these, and found analytical forms which approximate the spin polarisation well, with the magnetic potential they create simply being U ζ/2. Taking inspiration from how the Hartree potential is calculated in the atomistic model relative to the continuum model, a general expression for the potential induced by these anti-ferromagnetic instabilities would be where δ 1 corresponds to the constant sublattice polarization strength, δ 2 (G i ) is the moiré modulated part of the potential, i.e. the weight associated to the expansion of the potential in the i-th BZ reciprocal lattice vector, and τ z is the Pauli matrix.
For the MAFM instability both δ 1 and δ 2 are real valued.In the case of NAFM order there is not a constant contribution to ζ N so δ 1 = 0, and δ 2 is a real number.Similarly for the HAFM instability δ 1 = 0 but δ 2 is a purely imaginary number so the modulation is a sine function.Both MAFM and NAFM orderings are degenerate in the valley index but the HAFM is not and the modulated contribution to the potential must be complex conjugated when exchanging valleys.The value of the parameters δ 1 and δ 2 are obtained numerically within the self-consistent atomistic Hubbard calculation at a twist angle of 1.54 • and at charge neutrality, as described in These parameters should be linear in the interaction strength and the spin-polarised electron density, meaning an overall non-linear dependence on U , but we shall not seek self-consistent solutions to this potential in the continuum model here.Instead, these potentials are included in the continuum model through perturbation theory on the converged self-consistent Hartree-Fock calculations, as explained in Section II B. Note that the sublattice is polarised in the Hartree-Fock calculations, which means the overall sign of the potential can be chosen in two different ways, but we always choose the sign such that the sublattice polarisation matches, as this should increase any gaps, and further lower the energy.In the Appendix, we show the band structures in the continuum model (without the Hartree-Fock contribution) at 1.54 • with these perturbed potentials, and find good agreement with the atomistic calculations.
For twist angles smaller than 1.54 • , the Hubbard interaction in the continuum model should scale as 2 , [41].However, extending the parameters to the magic angle in this way does not yield very physical results, as the gaps in the flat bands are much larger than is ever found in experi-ments, as shown in Fig. A2(c), (f) and (i).This is because to obtain a mean-field solution of the magnetic order at 1.54 • a U = 5.4 − 6 eV was required.However, a more physical value is U = 4 eV, or smaller [75].Moreover, the large interaction strength gives rise to large polarised spin densities, as seen in Section III A, where the ζ values indicated that not just the flat band electrons were being polarised.Thus, we rescale these parameters to obtain more suitable estimates for these for the changes to the electronic structure.In Fig. 2 we show the results with a scaling factor of 3 in the main text, with other scaling factors being shown in the Appendix.As these perturbed calculations are not self-consistent, we cannot say what is the ground state.However, we can interpret the changes to the electronic structure, which provides information about where these interactions might be significant.First we shall discuss MAFM order, as shown in Figs.2(a), (b) and (c).At a twist angle of 1.54 • , Fig. 2(a), we find a small gap at the K/K' points, which is more significant than the Hartree-Fock distortions.For a smaller twist angle of 1.25 • , Fig. 2(b), a similar situation is found, the magnetic order opens a gap at the K/K' points while the Hartree-Fock potential is responsible for minor adjustments in the band structure.At the magic angle of 1.05 • , Fig. 2(c), the situation totally changes, now the magnetic potential only slightly modifies the gap at the K/K' points, while the Hartree-Fock contribution dominates the deformations to the electronic structure.
Next, we move on to describing the NAFM ordering, as seen in Fig. 2(d), (e) and (f).At the largest twist angle of 1.54 • , Fig. 2(d), this magnetic order creates a significant gap at the K/K' points.This could be attributed to the non-self-consistent nature of the calculations, and form using large values of the parameters.This large band deformations persists at the smaller twist angles of 1.25 • , Fig. 2(e), and 1.05 • , Fig. 2(f).Even if smaller values of the parameters are utilised, this NAFM order induces large band deformations, lowering the energy of the occupied valence band.Therefore, it appears that this magnetic order can compete with the Hartree-Fock contribution.
Finally, we describe the effect of HAFM.As can be seen in Figs.2(g), (h) and (i), the HAFM order does not create a gap at the K/K' points.Instead it causes the Dirac cones at K and K' to shift up and down, respectively, for the single spin and valley channel.At 1.54 • , Fig. 2(g), this effect is almost imperceptible but stronger than the one induced by the Hartree-Fock interactions.For the smaller twist angle of 1.25 • , Fig. 2(h) the situation is similar but the energy gap between K and K' points increases.While at the magic-angle, Fig. 2(i), it slightly contributes to reshaping the band structure, which in contrast is heavily affected by the Hartree-Fock contribution, so we can safely say that the HAFM is a secondary effect in this case.

IV. DISCUSSION
Overall, it appears that these magnetic potentials are more significant away from the magic angle, but at angles close enough that there could still be broken symmetry phases [26,64].For example, 1.25 • seems to be the most significantly affected by these Hubbard potentials relative to the Hartree-Fock contribution.At the magic angle, the Hartree-Fock contribution dominates, and at large twist angles the effects are small relative to the bandwidth, suggesting that these magnetic orders are not significant at these twist angles.This twist angle dependence could suggest why the predictions of Klebl et al. [64], in terms of the twist angle and doping dependence of magnetic states, agreed well with subsequent experiments [26].This could be because these Hubbard interactions are important close to the onset of broken symmetry phases, but close to the magic angle these Hubbard interactions are dominated by long-ranged Hartree-Fock interactions [41].
The NAFM order appears to effect the electronic structure most significantly, significantly lowering the eigenvalues of the occupied valence band, and therefore, it is a possible candidate for magnetic order in TBG.In contrast, the MAFM order effects the electronic structure more weakly.Finally, the HAFM appears to only effect the electronic structure slightly.This magnetic order should, however, couple to perpendicular electric fields [75,91], which could make this ordering tendency more important.These perturbative calculations are interesting to be performed because they are a very natural explanation for the correlated insulating states in TBG [4].From the Hartree-Fock calculations, a spinvalley degenerate insulating state is obtained.The atomistic Hubbard interaction, however, should break this symmetry and cause the onset of magnetic order.
We have focused on TBG here, but many more moiré materials comprised of graphene exist [94,95].Perhaps the most promising ones are where there is a ±θ twist between each adjacent graphene layer [93,[96][97][98][99][100][101][102].These moiré graphene multilayers have been shown to host highly tunable superconducting phases [103][104][105][106][107], and as the number of layers increases, the superconducting phase occurs over wider and wider doping ranges [108,109].Fischer et al. [100] has shown similar types of magnetic order occur in these systems, which means it will be possible to use the developed method.Another example of moiré structures is graphene twisted on a graphene multilayer, such as twisted monobilayer graphene [110,111].The magnetic structure of these systems was shown to be more complex by Goodwin et al. [112], which suggests the approach described here could be difficult to utilise.Finally, another class of moiré graphene multilayers is twisted bilayers composed of graphene multilayers, such as twisted double bilayer graphene [113][114][115][116][117][118].Further investigation of this system would be of interest.

V. CONCLUSION
In summary, starting from atomistic methods, we studied several leading magnetic instabilities of charge neutral TBG at a large twist.These calculations permitted analytical forms for the magnetic ordering.These Hubbard potentials were investigated perturbatively from selfconsistent Hartree-Fock calculations in the continuum model, allowing a comparison of long and short range exchange interactions.From these calculations, our takehome conclusions are: 1.These atomistic Hubbard interactions break the spin-valley degeneracy of the insulating state at charge neutrality of TBG obtained from selfconsistent Hubbard calculations.Therefore, these insulating states are likely to have some Mottness.
2. These magnetic orders are most significant for intermediate twist angles between the magic-angle and angles where non-interacting physics is sufficient.At the magic-angle, the Hartree-Fock contribution dominates, and at large angles the bandwidth dominates.
3. Out of the studied magnetic orders, nodal antiferromagnetic order appears to be the most significant for changes to the electronic structure.
It is hoped that these results further motivate inclusion of atomistic effects in the continuum model.Moreover, performing self-consistent magnetic calculations should also be possible, and investigating such ordering tendencies in other moiré graphene multilayers is possible now.
i are reciprocal lattice vectors, and the sublattice oscillation ζ s is the constant sublattice polarisation and ζ s describes how this sublattice polarisation changes on the moré scale.For sublattice A l (B l ), where l = 1, 2 is the layer index, the sign of the polarisation is − (+), or vice versa.Note this equation assumes the AA region is located at the origin of the moiré unit cell.In Fig. 1(b) and (c) the unconstrained, self-consistent ζ for 1.54 • at charge neutrality and U = 5.4 eV is shown.The peaks of ζ ≈ 0.1 reside in the AA regions.If one electron is delocalized across all atoms in the moiré unit cell, which contains 5548 atoms at 1.54 • , then ζ ≈ 10 −4 .

Figure 1 .
Figure 1.(a,d,g) Self-consistent quasi-particle band structures for the studied magnetic orders of TBG at 1.54 • and charge neutrality along the high symmetry path.For MAFM, U = 5.4 eV and unconstrained calculations were performed; whereas, for NAFM and HAFM, U = 5.94 eV and constrained calculations were performed.(b,e,h) Corresponding self-consistent magnetic order parameter plotted along the diagonal of the moiré supperlattice, where R1 and R2 are the moiré lattice vectors.Sublattice A is shown in black and sublattice B is shown in grey for the top graphene layer (bottom layer not shown).(c,f,i) Corresponding plots in real space for a single layer and sublattice, where only sublattice B of the top layer is shown.Note, the U 's were chosen to be slightly above the critical values.For MAFM, if U = 5.94 eV is used, the gap at the Dirac point is extremely large[91].
Parameters |δ1| and |δ2| extracted form the atomistic Hubbard calculations for the different orderings considered in this work.All values are expressed in meV.
VI. ACKNOWLEDGMENTS A.J.P., P.A.P. and F.G. acknowledge support from the Severo Ochoa programme for centres of excellence in R&D (Grant No. SEV-2016-0686, Ministerio de Ciencia e Innovación, Spain); from the European Commission, within the Graphene Flagship, Core 3, grant number 881603 and from grants NMAT2D (Comunidad de Madrid, Spain) and SprQuMat (Ministerio de Ciencia e Innovación, Spain).ZG was supported through a studentship in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (EP/L015579/1).We acknowledge funding from EPSRC grant EP/S025324/1 and the Thomas Young Centre under grant number TYC-101.We acknowledge the Imperial College London Research Computing Service (DOI:10.14469/hpc/2232) for the computational resources used in carrying out this work.This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie grant agreement No. 101067977.The Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) is acknowledged for support through RTG 1995, within the Priority Program SPP 2244 "2DMP" and under Germany's Excellence Strategy-Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC2004/1 -390534769.We acknowledge support from the Max Planck-New York City Center for Non-Equilibrium Quantum Phenomena.Spin susceptibility calculations were performed with computing resources granted by RWTH Aachen University under projects rwth0496 and rwth0589.