Magnetism of magic-angle twisted bilayer graphene

We investigate magnetic instabilities in charge-neutral twisted bilayer graphene close to so-called"magic angles"using a combination of real-space Hartree-Fock and dynamical mean-field theories. In view of the large size of the unit cell close to magic angles, we examine a previously proposed rescaling that permits to mimic the same underlying flat minibands at larger twist angles. We find that localized magnetic states emerge for values of the Coulomb interaction $U$ that are significantly smaller than what would be required to render an isolated layer antiferromagnetic. However, this effect is overestimated in the rescaled system, hinting at a complex interplay of flatness of the minibands close to the Fermi level and the spatial extent of the corresponding localized states. Our findings shed new light on perspectives for experimental realization of magnetic states in charge-neutral twisted bilayer graphene.

], (c) total density of states (DOS) corresponding to panel (b). The almost flat minibands at zero energy and corresponding large DOS peaks exhibit good agreement between the rescaled and non-scaled systems.

Geometry of twisted bilayer graphene (TBG)
Let us start by explaining the geometry of TBG in more detail. A single layer of graphene consists of carbon atoms arranged in a honeycomb lattice such that the unit cell includes two sites. We then construct a periodic commensurate bilayer structure parameterized by two integers m, n using the method of Refs. [27,28,31,38,39]. m and n are coordinates with respect to the lattice vectors of a single graphene layer a 1,2 = a( √ 3, ±1)/2. The rotation angle for such a commensurate structure (moiré pattern) is then given by cos θ = n 2 + m 2 + 4mn 2(n 2 + m 2 + mn) , (1) and the fundamental vectors of the TBG superlattice are t 1 = na 1 + ma 2 and t 2 = −ma 1 + (m + n)a 2 . The number of atoms in the moiré cell is given by N c = 4(n 2 + m 2 + mn) .

Model Hamiltonian
We start from the tight-binding model for the p z orbitals of the carbon atoms in chargeneutral TBG:Ĥ =Ĥ 0 +Ĥ int , whereĤ 0 is the single-electron Hamiltonian andĤ int is the electron-electron interaction. This leads to the one-band Hubbard model whered † iσ andd iσ are the creation and annihilation operators of an electron with spin projection σ = {↑, ↓} at site i andn i = σd † iσd iσ is the total electron density at site i. The hopping parameters t(r i ; r j ) between two p z orbitals located at r i and r j are given in Refs. [28,31]. The second term in Eq. (3) describes the on-site Coulomb repulsion. The resulting non-interacting band structure (U = 0) at the first magic angle θ = 1.08 • is shown by the full blue line in Fig. 1(b). This case corresponds to (n, m) = (30,31) and thus to a moiré cell with N c = 11164 sites. Dealing with such big unit cells will be challenging even for a one-band model and even within mean-field theory (MFT) and thus we will explore an idea of Ref. [40] to reduce the numerical effort.
The precise non-interacting band structure depends not only on the geometry, but evidently also on the hopping parameters t(r i ; r j ), and in particular the ratio between intraand interlayer hopping. Let θ and θ be the angles corresponding to two commensurate moiré structures and Then the rescaling t 0 = Λ t 0 of the nearest-neighbor intralayer hopping while keeping the interlayer hopping unchanged maps the low-energy band structure from the unprimed to the primed geometry [40]. The panels (b) and (c) of Fig. 1 illustrate this mapping for the first magic angle from θ = 3.89 • to θ eff ≡ θ = 1.08 • . Indeed, the dashed red line reproduces both the low-energy band structure and the density of states well at the expense of reducing the nearest-neighbor intralayer hopping from the physical value t 0 = 2.7 eV [5,41] to t 0 ≈ 0.75 eV, i.e., modifying the high-energy physics. With different rescaling factors, i.e., t 0 ≈ 0.90 eV and 1.02 eV, we can also model the angles θ = 1.30 • and 1.47 • in the (n, m) = (25,26) and (22,23) systems, respectively by the same effective (n, m) = (8, 9) system. Ref. [40] suggested that the on-site Coulomb interaction should scale in the same way as the intralayer hopping parameters, U = Λ U although this is less evident than the rescaling of the hopping parameters, as we will also see in the results to be presented below.
In the following section 4 we will first explore this rescaling trick in order to perform a detailed study using the case (n, m) = (8,9) (N c = 868). In section 5 we will then check for some representative cases to what extent the conclusions do indeed apply to the unscaled system, including the first magic angle, i.e., (n, m) = (30, 31) (N c = 11164).

Rescaled system
In this section, we investigate the Hubbard model (3) for twisted bilayer graphene (TBG) using rescaled interlayer hopping parameters, as outlined in the previous section. We will start with a systematic study using static MFT and then use a more sophisticated dynamical mean-field theory (DMFT) to argue that the findings of the simple MFT are qualitatively correct even if there is a quantitative renormalization of the values of the on-site Coulomb interaction U .

Static mean-field theory (MFT)
Static MFT is a well-established method to investigate the magnetism in graphene (see, e.g., chapter 3.1 of [5] and Refs. [35,36,40,42,43]) such that here we summarize only the essential features. It amounts to the Hartree-Fock approximation of the interaction term in Eq. (3), where n iσ is the average electron occupation number with spin σ at site i. Note that the approximation (5) decouples the operators for the two spin sectors and thus gives rises to a quadratic Hamiltonian in each of them where the other spin sector enters only via the site-dependent mean fields n iσ that have to be determined self-consistently. We focus on charge-neutral TBG that has exactly one electron per site, i.e., we work with the half-filled Hubbard model. A self-consistent solution is found iteratively, where in each step N c × N c matrices need to be diagonalized and an integral over the moiré Brillouin zone has to be calculated, that we approximate by a uniform grid of k points. We iterate this procedure until the maximum change of a density is below 10 −6 . Given the necessity to diagonalize a large number of moderately-sized matrices, even this elementary MFT approach becomes CPU-time intensive in the present situation. Some checks indicate that a k-grid of at least 9 × 9 points is required to eliminate artifacts of this discretization while more points do not change the conclusions. We therefore show results below that have been obtained for 9 × 9 k points.  The RPA analysis that we present in appendix A reveals different competing magnetic instabilities at different values of q for the present model. There is a periodic solution with an antiferromagnetic internal structure. The dominant instabilities are actually found at q = 0, i.e., they should have a larger unit cell than the twisted bilayer lattice, and they have a ferromagnetic structure inside a moiré cell. Motivated by the fact that the Hubbard model on a single honeycomb layer becomes antiferromagnetic at large U [44], we focus here on the antiferromagnetic mean-field solution. The RPA analysis of appendix A predicts a critical value U c ≈ 0.23 t 0 for the antiferromagnetic state of the twisted bilayer system with θ eff = 1.08 • , an order of magnitude below the critical value of a single layer, that for nearestneighbor hopping is known to be U MFT c /t ≈ 2.23 [44].
respectively. We first focus on the first magic angle θ eff = 1.08 • (red data in Fig. 2). Here, we find a small albeit finite magnetization for values of U/t as low as We note that convergence is delicate close to U c1,MFT and sensitive to the chosen k grid.
The result (8) should thus be considered as an upper bound. Thus, we conclude that this value is consistent with the prediction of the RPA analysis of appendix A. Given that the magnetization for these small values is due to the four flat minibands and that there is a low number of associated states (4 per moiré cell), the total magnetization (6) for small values of U is small and thus seen more clearly in the inset of Fig. 2(a) than in the main panel. Indeed, for U/t 0 1.5, the total magnetization per moiré cell remains below 2 = 4 · 1/2, consistent with it coming mainly from the four flat minibands.
An alternative perspective is given by the maximum magnetization (7) that is shown in Fig. 2(b). Here, one can firstly observe the onset of magnetization around U c1 more clearly than in the main panel of Fig. 2(a). For comparison, the main panel of Fig. 2(b) also includes the result for a single layer with the same intralayer hoppings as in the twisted bilayer system. One observes firstly that additional long-range hoppings within each layer reduce the critical value slightly to U MFT c /t ≈ 2.09 as compared to the nearest-neighbor result U MFT c /t ≈ 2.23 [44]. In the region U/t 0 2, the magnetization of the bilayer system is slightly enhanced with respect to the single-layer case, as might be expected thanks to the additional intralayer couplings. However, the transition to full magnetization necessarily involves AB and BA stacking regions (see Fig. 1(a)) that are geometrically frustrated. Consequently, one expects a complex magnetic state in this transition region. A full analysis of the transition to a fully magnetized system is beyond the scope of the present work, but we note that convergence is delicate also in this second transition region, as exemplified by the outlier at The most important finding in the present context is that magnetism arises in the effective twisted bilayer model at the magic angle for Coulomb interactions U that are an order of magnitude smaller than for decoupled single graphene layers. It should be noted that the q = 0 magnetic solution considered here only corresponds to a local, but not the global minimum of the energy such that the true critical value of U 1.08 • c1,MFT is probably even smaller than the result (8) (U 1.08 • c1,MFT ≈ 0.15 t 0 according to the RPA analysis of appendix A). Figure 2 also includes two examples for larger twist angles θ eff = 1.30 • and 1.47 • (green and blue data, respectively). Many of the preceding remarks also apply to these two cases such that we focus on their peculiarities. Remarkably, the case θ eff = 1.30 • yields an even smaller U 1.30 • c1,MFT /t 0 ≈ 0.21 than for θ eff = 1.08 • . Actually, while the velocity at the K point only vanishes at the first magic angle θ = 1.08 • , the minibands have a very small bandwidth over the entire range until θ = 1.30 • . However, when one goes to θ eff = 1.47 • , the critical value of the onsite Coulomb repulsion increases to U 1.48 • c1,MFT /t 0 ≈ 1.0. This is still significantly smaller than the critical value of a single layer U MFT c /t ≈ 2.09, but clearly larger than in the two other cases, as expected for minibands close to the Fermi level that now have both a finite Fermi velocity and a significant bandwidth.
For a more detailed discussion of the magnetic state found above U c1 but before the system becomes completely magnetic, we show in Fig. 3 results for the θ eff = 1.08 • system and a representative value of the on-site Coulomb interaction U/t 0 = 1. The top panel shows the spatial structure of the magnetization pattern that we find to be localized in the AA stacking region. Thus, in this region the magnetic state of the twisted bilayer system resembles that of AA stacked bilayer graphene, but at a significantly lower value of U than would be required for the simple AA system to become magnetic. A different perspective of this magnetic pattern is provided by the lower right panel of Fig. 3 that presents a diagonal line cut of the magnetization. The lower left two panels of Fig. 3 show the spin-resolved local density of states (LDOS) in the AA and AB stacking regions. Interestingly, in the AA region one finds two peaks in the LDOS at low energies that are absent in the AB stacking region. The presence of these peaks correlates with the magnetic state, thus rendering scanning tunneling SciPost SciPost Phys. 11, 083 (2021) microscopy (STM) experiments a promising candidate for the detection of such a magnetic state.

Dynamical mean-field theory (DMFT)
Even though MFT has been shown to be remarkably successful to qualitatively describe static [42] and dynamic properties [43] in the semi-metallic phase of single-layer graphene, it is known to become quantitatively less accurate for larger values of U . For example, the transition to the antiferromagnetic insulator in the nearest-neighbor hopping case is found at U MFT c /t ≈ 2.23 in MFT [44] while more sophisticated and accurate methods place it at a larger U c /t ≈ 3.8 [45][46][47][48].
DMFT [34] takes local charge fluctuations into account and thus improves the quantitative treatment of the on-site Hubbard interaction. Indeed, already single-site DMFT shifts the estimate of the critical point to the range U DMFT c /t = 3.5, . . . , 3.7 [35], i.e., remarkably close to the most accurate estimates [45][46][47][48]. Following previous work, we employ here a real-space version of DMFT [36]. DMFT maps the lattice Hamiltonian Eq. (3) onto a set of quantum impurity problems via the local Green's function for site i inside the moiré supercell [34] HereĤ 0 is the single-particle part of Eq. (3). The main approximation is that the local self-energy matrix for spin projection σ, Σ r σ (z), that plays the role of a dynamical mean field, depends only on frequency z, but not on momentum k. Eq. (9) can be used to define a collection of N c single-impurity Anderson models, that we solve here with the numerical renormalization group (NRG) [49][50][51] and iterate until self-consistency is reached [52,53]. We refer to Refs. [35,36] for details on the procedure and just mention two peculiarities for the present case. Firstly, Eq. (9) requires evidently a combination of integration over the moiré Brillouin zone while at the same time solving coupled problems for the N c atoms inside the moiré supercell. Secondly, even if the band structure of Fig. 1(b), (c) is almost particle-hole symmetric, there is no strict particle-hole symmetry in the present case in contrast to previous work [35,36]. Consequently, the chemical potential needs to be adjusted appropriately during each iteration in order to ensure half filling. Since the chemical potential enters into Eq. (9) in a non-linear fashion, this renders the numerical problem even more challenging, thus limiting DMFT not only to the rescaled system, but also the number of U -values considered. Figure 4 presents some DMFT results for the rescaled system with θ eff = 1.08 • , i.e., N c = 868. Comparison of the DMFT results for the magnetization versus U/t 0 in Fig. 4 with the MFT results of Fig. 2 shows qualitatively similar behavior. At a technical level, the DMFT results are a bit more noisy. This is due to the logarithmic frequency discretization inherent to NRG [51], to DMFT being generally numerically more expensive, and in particular the difficulty to adjust the chemical potential appropriately. Nevertheless, the main quantitative difference remains that the critical U c of a single layer is pushed to larger values as compared to simple MFT, and so is the phenomenon of a magnetization arising in the AA stacking region of the twisted bilayer system. Nevertheless, also DMFT clearly detects a magnetization in the twisted system for values of the local Coulomb interaction down to U/t 0 = 1, amounting to a reduction of the critical value as compared to the single-layer system by at least a factor 3.5 at θ eff = 1.08 • . The inset of Fig. 4 shows an example of the spatial magnetization pattern. This is again very similar to the MFT result shown in the main panel of Fig. 3, just the value of U/t is renormalized to larger values, namely from 1 for the MFT example to 2.5 of the DMFT example. Note that U/t 0 = 2.5 would give rise to a bulk magnetic state within MFT while the DMFT result in the inset of Fig. 4 is still clearly localized in the AA stacking region. Overall, DMFT confirms the qualitative conclusions derived from MFT; it just provides a quantitatively more accurate account of the local Coulomb interaction U .

Non-scaled system
We will now present some results for the non-scaled system. The scaling trick has allowed us to apply the quantitatively more accurate DMFT, but the size of the moiré cells of the non-scaled systems will exceed those accessible to DMFT such that we focus on static MFT in the present section. We use the same parameters as in section 4.1 (convergence criterion 10 −6 , 9 × 9 k-grid). Figure 5 shows MFT results for the total magnetization per moiré cell as a function of U/t 0 at rotation angles θ = 1.08 • , θ = 1.30 • , and θ = 1.47 • . The corresponding moiré cells contain N = 11164, 7804, and 6076 carbon atoms, respectively. At first sight, the behavior is very similar to that found in the inset of Fig. 2(a) for the rescaled system (the smaller number of data points is due to the significantly enhanced computational effort). In According to Ref. [40], in the given normalization, these values should correspond to those found in the rescaled system. This works out more or less for the case θ = 1.47 • where in both cases, the critical U/t ratio is close to 1. However, the values for U 1.08 • c1,MFT and U 1.30 •

c1,MFT
in the non-scaled system are bigger than those we might have expected from the rescaled case. Indeed, the order of the discrepancy corresponds to another factor Λ such that U scales with Λ 2 and not just with Λ. A possible interpretation of this observation is the following: Λ actually also appears in the scaling of the linear length [40]. Now the magnetic instability at the angles θ = 1.08 • and 1.30 • is related to a state localized in the AA region, see, e.g., top panel of Fig. 3. Thus, the area of the relevant spatial region scales with Λ 2 , accordingly the number of contributing local on-site repulsions also scales with Λ 2 such that U should also scale with Λ 2 rather than Λ in the cases where the physics is controlled by localized states. In spite of this additional factor, it remains true that U 1.30 • c1,MFT /t 0 ≈ 0.55 < U 1.08 • c1,MFT /t 0 ≈ 0.85, and that there is still a significant reduction by factors of 4 respectively 3 with respect to the critical value U c for a single graphene layer. In light of the preceding observations, we suggest that not only the non-interacting bandwidth, but also the size of the moiré cell matter. While both the θ = 1.08 • and 1.30 • bilayers have a small bandwidth, the moiré cell of the latter is smaller, and this appears to result in a smaller critical value of U c1 . The size of the moiré cell is smallest for θ = 1.47 • among the three cases studied, but the value of U c1 is biggest, most likely because in this case the minibands closest to the Fermi energy are no longer flat. Nevertheless, even in this case one observes emergence of magnetism for values of U that are about a factor 2 smaller than would be needed to render a single layer antiferromagnetic.
To conclude this discussion, let us have a closer look at the spatial structure of the resulting magnetic states. Figure 6 shows the spatial magnetization profile for non-scaled moiré unit cells with angles θ = 1.47 • and at the first magic angle θ = 1.08 • . For illustration purposes, we consider a value of U just above the first critical point U c1 , i.e. U/t 0 = 1.50 and 1.00, respectively. Like for the recaled system shown in the top panel of Fig. 3, we find an antiferromagnetic pattern that is localized in the AA region. However, thanks to the improved spatial resolution, we can now observe a clearer separation of the magnetic regions between neighboring cells. For slightly larger U , the magnetic region grows, but the structure remains qualitatively similar as in Fig. 6.

Conclusions and perspectives
We have investigated the onset of magnetism in charge-neutral "magic-angle" twisted bilayer graphene with numerical real-space static and dynamical mean-field approaches. In the rescaled system we found that localized magnetic states appear in the twisted bilayer system for values of the on-site Coulomb repulsion U that are an order of magnitude smaller than those needed to render a single layer magnetic. We then showed that the non-scaled system exhibits qualitatively similar behavior. The reduction is less impressive (up to a factor 4 in the cases investigated), but still remarkable. We note that this finding is consistent with a recent diagrammatic real-space mean-field study [54] that focused on two selected values of U .
The rescaling proposed in Ref. [40] actually reproduces the flat minibands close to the Fermi level very well, compare Fig. 1. Our results therefore demonstrate that the band structure is not the only factor that matters. Indeed, the corresponding states are localized in AA stacking regions. This suggests a scaling of the critical U c with area rather than linear size, as is indeed roughly consistent with our findings for θ = 1.08 • and 1.30 • . A more quantitative analysis would involve computation of the Coulomb matrix elements with respect to the Wannier functions [11,15,[55][56][57][58] of the rescaled and non-scaled systems, respectively. However, such an analysis goes beyond the scope of the present work.
A side effect of the observation that the spatial extent of the localized states also matters is that smaller unit cells favor magnetism over bigger ones. Indeed, we find onset of magnetism for θ = 1.30 • for smaller values of U c than for the first magic angle θ = 1.08 • . The system with θ = 1.47 • has an even smaller unit cell than that with θ = 1.30 • , but at this larger angle there is no longer any really flat band close to the Fermi level such that here the value of U c is found to be larger. A related point is that magic angles are usually defined via a vanishing Fermi velocity [27][28][29][30][31] while in fact it may be more relevant that the entire minibands are narrow. Indeed, the latter criterion is satisfied over the entire range θ = 1.08 . . . 1.30 • such that the smaller unit cell can then give rise to a lower U c at the upper boundary of this range of angles.
It should be noted that in our mean-field investigations we have focussed on antiferromagnetic solutions that are periodic over moiré cells. However, the RPA analysis of ap- pendix A suggests that there are other competing instabilities, and indeed the mean-field self-consistency loop sometimes converges to other solutions. In particular, the true lowest-energy state might be modulated in real space and exhibit an internal ferromagnetic structure, like in the case of an electric bias between the two layers [40]. Should this indeed be the case, this can only further reduce the value of the U c for the onset of magnetism such that our estimates are in fact upper bounds. The main conclusion that twisting leads to a significant reduction of the critical U c for the appearance of magnetism is thus unaffected by the assumptions on the nature of the ground state.
Another point to note is that we find a stronger reduction of the critical interaction U c at charge neutrality than a previous RPA investigation [37]. This difference can be traced to a different tight-binding model at the starting point. Indeed, the authors of Ref. [37] have implemented the corrugation of Ref. [55] that takes a modulation of the distance between the two layers in different stacking regions into account. However, other factors may also be relevant in an experiment such as strain when the bilayer is deposited on a substrate. In the same spirit, Coulomb interactions should actually be long-range [59], at least for freestanding bilayers, since atomically thin layers cannot screen the Coulomb repulsion between electrons. Still, screening will depend on the actual substrate and may thus depend on the exact experimental conditions. Even other factors such as spin-orbit interactions that are sufficiently weak to be usually negligible in graphene may matter in the present situation given the significant reduction of the kinetic energy scale in the twisted bilayer system. Thus, which of several competing instabilities finally wins in an experimental realization may depend on a number of factors; here we have simply demonstrated that a magnetic instability (possibly an antiferromagnetic one) is one of the competitors in charge-neutral twisted bilayer graphene.
The macroscopic magnetization of a ferromagnetic state can be detected, e.g., via the Hall effect [32]. Antiferromagnetic or almost ferromagnetic, but modulated spiral states are more difficult to detect experimentally since they do not give rise to a macroscopic moment. In bulk systems, one would resort to (neutron) scattering to detect such states, but in the present nanoscopic setting this may not be feasible. The best option may thus be scanning tunneling spectroscopy (STS) experiments [60] in order to detect the corresponding characteristic features in the local density of states (see lower panels of Fig. 3). In fact, the corresponding signatures might already have been observed in recent STS experiments [61][62][63]. However, the latter samples are subject to heterostrain [64,65] which also gives rise to a splitting in the electronic density of states. An unambiguous detection of a magnetic state would thus require a detailed investigation of the variation of the tunneling spectrum with the different stacking regions.
Returning to theoretical questions, an alternative approach would be via low-energy continuum models in the spirit of Ref. [30]. One reason why we have rather used the rescaled model [40] here is that, as illustrates Fig. 1, it reproduces the band structure well over a wide range of energies and not just the flat minibands close to the Fermi level. However, in the range of intermediate values of U/t where mainly the flat minibands contribute to the magnetism, effective low-energy models would have the advantage of being more amenable to numerical approaches [14,[66][67][68][69] such that we suggest the investigation of magnetism by such methods as a topic for further studies.
A further interesting issue that goes beyond questions accessible to low-energy effective models would be the full phase diagram of the twisted bilayer systems up to larger values of U/t. Indeed, the results underlying Fig. 2 suggest that there is no single simple transition to a bulk magnetized system, but that this transition actually proceeds via several intermediate states in the region U/t 0 ≈ 2. Given that magnetic interactions in the AB and BA stacking regions are geometrically frustrated (compare Fig. 1(a) A Noninteracting susceptibility of the rescaled system In this appendix, we provide an RPA analysis of the noninteracting susceptibility that is similar in spirit to Ref. [37]. However, here we focus on the rescaled system with θ eff = 1.08 • .
We adopt the multiorbital RPA approach to study the instability of the paramagnetic state [70,71].
The multiorbital spin susceptibilities tensor can be expressed in terms of the Matsubara spin-spin correlation function: with the Matsubara frequency ω, the imaginary time τ and spin operatorsŜ at orbitals s, t. The noninteracting (zero-order) susceptibility is just a simple bubble diagram involving two Green's functions. Using the spectral representation of the Green's functions, this can be expressed as where µ, ν are band indices, a s µ (k) and E µ (k) are the µ-th eigenvalue and eigenvector of the noninteracting Hamiltonian, respectively, and n F is the Fermi-Dirac distribution function.
The Coulomb interaction can then be included at the mean-field level and one arrives at a so-called "RPA" (or "Stoner", see, e.g. Refs. [26,72]) formula for the interacting susceptibility χ(q, iω) = χ 0 (q, iω) where in the paramagnetic state we can use χ 0 according to Eq. (11). According to Eq. (12), the static susceptibility χ(q, iω = 0) diverges whenever U equals one of the eigenvalues of the tensor χ 0 (q, iω = 0) −1 . One can use this identity to determine the mean-field critical U c , and indeed, the critical U c of an infinite graphene sheet was originally determined in this manner [44]. The value of q and the corresponding eigenvector yield information about the expected magnetic state for U > U c . The tensor of Eq. (11) is symmetric, but computing all N 2 c entries for a fixed q is timeconsuming since each of them involves a sum over reciprocal space and two sums over all energy levels. In order to keep the CPU time manageable, we have limited the sum µ,ν to states that are close to the Fermi energy. The latter approximation is physically justified since the ground-state ordering should be dominated by the quasi-flat bands close to the Fermi energy. A similar approximation to the Matsubara sums has also been used in Ref. [37] except that we use here a more radical sharp cutoff. Nevertheless, we have checked that taking the 50 to 100 states closest to the Fermi energy into account is sufficient to yield no visible truncation effects; we used 200 states to be on the safe side. Since we use a finite grid for the integration over the moiré Brillouin zone, the sum Eq. (11) consists strictly speaking of a finite number of poles. In order to smooth these out, we introduce a broadening parameter and evaluate Reχ 0 (q, iω = iη) such that we obtain a Lorentzian broadening of width η at iω = 0. Apart from the truncation in energy space, the momentum grid, and the broadening parameter η, the result for χ 0 (q, iω = 0) also depends on temperature T . T = 10 −8 t seems to be sufficiently low to ensure ground-state physics. However, there is a delicate balance between broadening parameter η and the grid in reciprocal space. If η is too large, it will smear out any peaks and thus reduce the values of χ(q, iω = 0). On the other hand, for a too small value of η, the momentum discretization will become visible. We found that the combination η = 5 · 10 −5 t and a uniform 9 × 9 grid of points (k x , k y ) yield a good compromise such that we will present results for these parameters here. Figure 7(a) shows the distribution of the leading eigenvalue of the static susceptibility tensor χ 0 (q, iω = 0) in the moiré Brillouin zone. In contrast to single layers and AA-stacked bilayer graphene that prefer a single type of ordering at q = 0, in the present case the maximal eigenvalue of χ 0 (q, iω = 0) is rather flat in reciprocal space. The global maximum is neither at q = Γ nor at the two symmetry-related points K and K , but rather at another point q max at the boundary of the first Brillouin zone. The values are max χ 0 (q, iω = 0) = 4.35378/t 0 , 4.99189/t 0 , and 6.61892/t 0 , for q = Γ, K, and q max , respectively. According to the discussion around Eq. (12), this predicts a critical value U c = 0.229686 t 0 for a q = 0 state and globally U c = 0.151082 t 0 , but for a state with a spatial modulation with a wave vector q max over moiré cells.
Panels (b)-(d) of Fig. 7 show the corresponding eigenvectors of the susceptibility tensor. At the Γ point (panel (b)), one observes a staggered sign change between nearest neighbors with the maxima located in the AA region. This corresponds to the periodic antiferromagnetic state that we have investigated in the main text. The analogous result for the eigenvector at the K (K ) point is shown in Fig. 7(c). Here we find a ferromagnetic solution in each moiré cell with the maximum again in the AA region, but the corresponding value of q implies that the corresponding state should be accompanied by a tripling of the unit cell in real space. Finally, Fig. 7(d) shows the eigenvector at q max . The local structure inside a moiré cell is still ferromagnetic, but exhibits a stronger internal modulation. Furthermore, the corresponding MFT solution should be modulated with a wave vector q max in real space. Examination of further values of q reveals an antiferromagnetic internal structure close to the Γ point while the ferromagnetic internal arrangement is predominant in other regions of the Brillouin zone.  Figure 7: (a) Distribution of the largest eigenvalue χ 0 (q, 0) of the susceptibility tensor for the rescaled system with θ eff = 1.08 • . The white hexagon denotes the first Brillouin zone. Panels (b), (c), and (d) show the spatial profile of the largest eigenvector of the static susceptibility tensor for q = Γ, K, and q max , respectively. We used η = 5 · 10 −5 t, a uniform 9 × 9 grid to evaluate the sum k , and a total of 200 states around the Fermi level for each µ and ν.