Functional renormalization group for a large moir\'e unit cell

Layers of two-dimensional materials arranged at a twist angle with respect to each other lead to enlarged unit cells with potentially strongly altered band structures, offering a new arena for novel and engineered many-body ground states. For the exploration of these, renormalization group methods are an appropriate, flexible tool that take into account the mutual influence of competing tendencies. Here we show that, within reasonable, non-trivial approximations, the functional renormalization group known from simpler two-dimensional systems can be employed for the large-unit cell moir\'e superlattices with more than 10.000 bands, remedying the need to employ ad hoc restrictions to effective low-energy theories of a few bands and/or effective continuum theories. This provides a description on the atomic scale, allowing one to absorb available ab-initio information on the model parameters and therefore lending the analysis a more concrete quantitative character. For the case of twisted bilayer graphene models, we explore the leading ordering tendencies depending on the band filling and the range of interactions. The results indicate a delicate balance between distinct magnetically ordered ground states, as well as the occurrence of a charge modulation within the moir\'e unit cell for sufficiently non-local repulsive interaction.


I. INTRODUCTION
In recent years the field of two-dimensional materials has made major experimental and theoretical leaps, which led to many fascinating discoveries and have broadened our spectrum on available phases of matter in these highly controllable structures. Examples of these novel findings include superconducting [1,2] or magnetic [3] phases realized down to the monolayer limit and, related to that, the discovery of quantum anomalous hall behavior in thin films [4][5][6][7][8][9][10][11], with potentially far-reaching technological applications in the realm of spin-tronics and quantum computing.
Recently, the twist angle between two sheets of material stacked atop each other was added as a further interesting research direction. In these twisted structures it was shown that, by the emerging huge real space moiré supercells (tiny Brillouin zone), kinetic energy scales can be reduced drastically, giving rise to prominent interaction effects. In Refs. 12 and 13 it was experimentally demonstrated that using this route of control two sheets of twisted bilayer graphene can be tuned to exhibit insulating as well as superconducting behavior. This exciting finding has spurred an enormous wave of experimental [14][15][16][17][18][19][20][21][22][23] as well as theoretical [24][25][26][27][28][29][30][31][32][33][34] research. Part of the excitement is founded in the fact that these twodimensional systems can be controlled with relative ease using a backgate, strain or the value of the twisting angle as a parameter. This allows one to access correlation regimes which might be otherwise difficult to access in a structure with a chemically straightforward composition (in this case graphene). The astonishing similarity of the phase diagram of twisted bilayer graphene to those of the cuprates has raised the hope that we might be able to understand the mechanism driving high-temperature superconductivity using the much simpler twisted graphene system. From a theoretical side a first approach is to concentrate on the correlation physics if we restrict ourselves to the low-energy bands near the Fermi surface. However, the validity of such an approach is difficult to assess. Often, when considering twisted van der Waals materials (like in the case of twisted bilayer graphene), there are no band gaps separating the lowest from higher energy bands and such a separation is unclear. Furthermore, topological arguments might obstruct the construction of such a simple low-energy theory [35]. In addition, a recent theoretical study shows that the interplay of correlations might be very subtle and fragile [36]. This calls for a different vantage point where, in contrast, we do not want to restrict ourselves to effective low-energy theories neglecting most of the many back-folded bands arising from the moiré supercell. When considering such a theory we immediately face the problem that at small twist angles many thousand bands need to be kept. The sheer number of bands makes theoretical descriptions very cumbersome, especially when one tries to include interaction effects. Recent advances in this direction include Refs. 37 and 38, which treat correlation effects on the randomphase approximation (RPA) or mean-field level.
Here, we want to add to this by establishing a more sophisticated tool, the so-called functional renormalization group (fRG), applied to a Hamiltonian keeping the many bands in the moiré Brillouin zone without reducing to effective low-energy theories. The fRG is a versatile method capable to describe a plethora of interacting electron systems [39][40][41][42][43]. Yet, solving the full fRG equa-arXiv:2002.11030v1 [cond-mat.str-el] 25 Feb 2020 tions after a truncation at the two-particle vertex (Γ (4) ) is in general numerically feasible only for systems with a few orbitals per unit cell, since the vertex function itself scales with the number of orbitals to the fourth power. In addition, the similarly rich dependence on momentaor the unit cell positions when formulated in real spaceremains a numerical challenge. Fortunately, for the latter dependence, reasonable simplifications can be formulated, with justifiable restrictions to short-ranged fermion bilinears [40,[44][45][46][47]. Furthermore, at least for a larger set of questions, additional approximations simplifying the orbital dependence can be made [48,49]. Here we employ these two approximation steps to the fRG equations in order to keep the numerical effort low enough to treat systems with more than ten thousand orbital sites per unit cell. We show that the treatment of twisted graphene bilayers close to the so-called 'magic angle' using these approximations, resolving all individual carbon sites in the large moiré unit cell, yields similar results to what we know from our previous study using the RPA of the crossed particle-hole channel [38].
This work extends such methodology to nonlocal interactions and the coupling to other channels beyond the RPA, and of course justifies the use of the RPA for the dominant instability a posteriori. We summarize the main results of this paper in the tentative phase diagram obtained with this method for a Hubbard interaction in Fig. 1. We find two magnetic orderings: First, antiferromagnetic order on the atomic scale with a sign change of the order parameter around the AA regions and, second, ferromagnetic order. The latter is only found at small interaction strengths and for fillings close to the van Hove singularities of the material's flat bands at low critical scales. Nodal antiferromagnetism is present for all fillings that show an instability at larger couplings. With increasing nonlocal contributions to the interaction, the instability gets dominated by the repulsion among electrons outside the AA regions of the moiré unit cell, indicating an interaction-induced charge redistribution.
The rest of this paper is structured as follows: We will shortly introduce the tight-binding model we use to describe 'magic angle' twisted bilayer graphene in Sec. II. Thereafter we discuss the method and approximations to it in Sec. III. Sec. IV shows the results of our simulations. Finally, we finish with some concluding remarks in Sec. V.

II. MODEL
We set up the moiré unit cell by constructing the superlattice vectors from lattice vectors of the honeycomb lattice of one single graphene sheet, l 1 = ( √ 3/2, 3/2, 0) and l 2 = ( √ 3, 0, 0). The first superlattice vector can be written as L 1 = nl 1 + ml 2 with n and m integers defining the twist angle. The second superlattice vector L 2 is rotated by 60 degrees with respect to L 1 . One of the layers is a honeycomb lattice with Bravais lattice vectors l 1 FIG. 1. Tentative phase diagram from the Γ-point IOBI fRG for twisted bilayer graphene at angle 1.05 • using a Hubbard interaction. The squares mark nodal antiferromagnetic order, the circles ferromagnetic order. For non-diverging flows a gray square is shown. The color indicates the critical scale ΛC that roughly corresponds to the critical temperature of the phase transition. The orange rectangle shows the parameter region which simulations with longer ranged interactions have been carried out for (cf. Fig. 6). At strong interactions, the nodal antiferromagnetic phase is favored, whereas at smaller U , the system is susceptible to ferromagnetic ordering if its filling is close to the van Hove singularities of the flat bands. and l 2 , the other one is shifted vertically and rotated by the twist angle θ = arccos m 2 +n 2 +4mn 2(m 2 +n 2 +mn) around the AA site. We implement corrugation effects by varying the interlayer distance in the supercell as described in Ref. 32. Choosing n = 31, m = 32 sets θ = 1.05 • and leads to 11908 sites in one supercell. The hopping parameters are taken from Ref. 50. Once we have set up the Hamiltonian for one specific Bloch momentum k of 'magic angle' twisted bilayer graphene, we diagonalize the matrix and obtain the spectrum b (k) and the orbital makeup u ob (k). The low energy part of the non-interacting band structure shows four (eight including spin) flat bands around charge neutrality (cf. Fig. 2). From the spectrum and orbital makeup, we can construct the free (Matsubara) Green's function as a matrix in orbital space: As electron-electron interaction, we employ the Ohno ansatz [51,52] with an exponential cutoff that can be motivated by screening of the substrate [53]. Its functional form in real space reads where we introduced the Ohno radius r O , the decay width and the limit µ O → 0, i.e. a Hubbard interaction. Note that there are more specific descriptions for the non-local interaction, also including the environmental screening [54]. For this qualitative, first fRG study we concentrate on the simple form Eq. (2) and mainly study the dependence on µ O . The interaction is of density-density type and can initially be written as a vertex function with components in the D channel only.

III. METHOD
We use the truncated unity approximation to the fRG equations which is described in more detail in Refs. 45 and 47, building essentially on earlier works, mainly Refs. 40 and 44. The main idea of this scheme is to write the interaction vertex as a sum of the following three channels: direct particle-hole (D), crossed particle-hole (C) and particle-particle (P ) channel. Each part can be understood as interaction between fermion bilinears of the corresponding types. The spatial structures of these bilinears can then be expanded in basis functions with specific symmetries in the moiré Brillouin zone. This expansion is truncated in its length. As a simple but nontrivial approximation we here consider only the in-cell contribution, i.e. the two fields in the fermion bilinears in pairing, charge and spin channels are in the same unit cell. Given the large extent of the moiré unit cell, this already captures quite some spatial dependence and hence may be a tolerable approximation. In addiotn, we then employ a Γ-point approximation of the vertex function in the momentum argument of the three channels. This means that we ignore all momentum dependence of the interaction within the small moiré Brillouin. One can convince oneself that is a reasonable approximation as long as the interaction decays significantly in real space through the moiré unit cell. Due to the Γ-point approximation, we are sensitive with respect to order parameters that vary spatially within the moiré unit cell but that do not change when translated into other unit cells. Again, due to the large unit cell, this does not render the analysis trivial, but it excludes, e.g., density waves that enlarge the moiré unit cell. These can be investigated by allowing other wavevectors besides Γ, which we however postpone to later work. In addition, we neglect self-energy contributions as well as the frequency dependence of the vertex function and truncate the fRG equations at the four point vertex Γ (4) . Since we want to treat SU (2) symmetric systems, we only need to consider the flow equations for the symmetrized vertex function V Λ where Λ is the flow parameter. The three channels obey coupled differential equations in the flow parameter Λ and carry orbital indices explicitly: with the particle-particle and particle-hole loops The projection operators P reduce to unity without further approximations since they affect the momentum dependencies we disregarded. In the zero-temperature limit, the Matsubara sums in Eq. (4) become integrals.
As we only need the differentiated loops d/dΛ L PP/PH,Λ , it is in some cases possible to evaluate the integral analytically. For the sharp cutoff we employ in our simulations, we use the Green's function G Λ o1o2 (k) = G (0) o1o2 (k) Θ(|k 0 | − Λ) and are able to trivially carry out the frequency integrals in the derivative of L PP/PH,Λ via the resulting delta function.
The approximation we make in order to be able to treat systems with a large number of orbitals is to only allow orbital bilinear interactions in each of the three interaction channels (IOBI approximation) [48]. This lowers the complexity of each channel of the vertex to only be number of orbitals to the power of two instead of four. This simplification reads The full vertex is composed by simply adding up the three channels: V Λ still is a fourth rank tensor in orbital space in the IOBI approximation, even though each channel is reduced to be a matrix in the orbital indices. The projections of the full vertex V Λ to a channel are now simply given by the restrictions of orbital bilinearity: Insertion of the channel projections into Equation (3) yields the flow equations in the static IOBI Γ-point approximation: +V P D ,ΛLPH,ΛV P C ,Λ +V P C ,ΛLPH,ΛV P D ,Λ .
The quantities are all matrices and connected by matrix products. The projected vertex function's matrices read By the IOBI approximation, the loops are constrained to be matrices as well and take the following form: where the operationÂ •B is an element-wise matrix product of the matricesÂ andB.

IV. FUNCTIONAL RENORMALIZATION GROUP RESULTS
Starting from a value Λ ini = 10 eV for the frequency cutoff, we integrate the flow to Λ 0 = 0. The value of the frequency step dΛ as a function of Λ is adaptively chosen in the window given by two envelope functions as shown in Fig. 3. Quite generally, the fRG flow leads to strong coupling, i.e. a rapid growth of certain components of the flowing interaction. We draw physical conclusions from this by estimating the maximum eigenvalue of each interaction channel for each step [55] for the intraorbital bilinear matricesP Λ ,Ĉ Λ andD Λ . If one of the eigenvalue estimates surpasses a critical value (which we set to 10 3 eV), we stop the flow and do an eigendecomposition of the three channel matrices. The eigenvector corresponding to the maximum eigenvalue indicates the orbital character of the order parameter associated to the divergence. An example of several fRG flows at fixed filling and interaction strength for different interaction ranges is shown in Fig. 4. Increasing the range of the interaction, µ O , leads to a divergence in the D channel compared to the magnetic instability in the C channel found for short ranged interactions.
The three types of ordering found as run-away flows in this work are shown in terms of their leading eigenvectors in Fig. 5: (a) ferromagnetic order throughout the unit cell with some residual antiferromagnetism on the carboncarbon-bond scale in the AB regions, (b) nodal antiferromagnetic order with a sign change of the antiferromagnetic order parameter on the atomic scale around the AA regions and (c) charge-modulated states that form a honeycomb lattice of the AB and BA regions. This last instability is characterized by a leading eigenvector in the density-channel. It should be analyzed in its pure form for µ O = 1 where the sign change in the eigenvector occurring for smaller µ O has disappeared. Then the charge modulation arises from the leading eigenvector that has largest weight in the AB-regions. This represents a strong growth of the effective electronic repulsion in those regions, which should then push the charges away from there into the AA-regions.
The instabilities of magnetic type occur for shorterranged interactions with smaller µ O and agree with the types of ordering we found in our previous RPA study for the magnetic susceptibility due to a pure onsite interaction [38]. Additionally, the values of Λ C are similar in their order of magnitude to the critical temperatures found in RPA. Increasing the range of the interaction leads to charge modulated (CM) states as dominant instability. In the context of inheriting the instabilities FIG. 6. Dependence of the phase diagram of magic angle twisted bilayer graphene on the range of the interaction µO. Squares mark nodal antiferromagnetic order, circles ferromagnetic order, plus symbols CM states without sign change and crosses CM states with sign change. Gray squares indicate that the vertex did not diverge before reaching Λ = 0. For longer ranged interactions, the critical scale rises while the magnetic orderings get replaced by charge-modulation orderings.
from AA and AB stacked graphene bilayers, charge density wave states with opposite charge density modulation on A and B carbon sublattices are expected at longer ranged interactions [56,57], but here the CM states fea- Ferromagnetism is present for fillings close to the van Hove singularities of the flat bands and small interaction strengths. Increasing U or doping away from the van Hove singularities leads to nodal antiferromagnetism. For longer ranged interactions, we find that the system is susceptible to a charge modulation (CM) instability. For the fillings and interaction strengths indicated by the orange square in Fig. 1, we carried out the simulations with longer range interactions. The results are summarized in a series of additional tentative phase diagrams in Fig. 6. We observe that for the cases where magnetism is found, the critical scale Λ C is almost independent of whether the interaction is longer ranged (µ O = 0.1, 0.2). Starting with µ O = 0.5, we observe that CM ordering becomes relevant for large U (8 eV). While the quantitative location of this AF-to-CM transition has to be checked in more elaborate numerical studies, our current finding indicates a high sensitivity of the interacting system with respect to non-local interaction, as also stated in Ref. [36]. For µ O = 0.5 and using Eq. (2), the nearest-neighbor repulsion V 1 is only about one third of the onsite repulsion U . The experimentally studied systems should have larger V 1 /U -ratios [54]. Additionally, the critical scale takes much larger values for this instability (see Fig. 7). At even longer ranged interactions (µ O = 0.7, 1.0), the magnetism subsequently makes way for the CM instability for large U . One exception is the charge neutral point, where no instability is found for U = 4 eV for short ranged interactions but at µ O = 0.7, a CM state is found at low Λ C . For µ O = 1.0, the critical scale is very similar across the whole parameter range tested in our simulations.

V. CONCLUSION
We use the functional renormalization group to describe interactions and their ordering tendencies in large moiré unit cells. The results for twisted bilayer graphene at magic angle with pure onsite interactions agree with our previous RPA study. Using the fRG allows us in addition to include long-range interactions. Here, we find charge modulated states besides magnetic orderings and investigate the competition between the different ordering tendencies.
Until now, only s-wave correlations are allowed by our approximations. Even in this approximation scheme, we are able to show how different magnetic order parameters emerge, some of which are inherited from nontwisted graphene bilayer systems (antiferromagnetism on the atomic scale). Additionally, the longer ranged Ohno interaction that is accessible in fRG leads to a competition between direct-and crossed particle hole channel, i.e., charge and magnetic fluctuations. As a strong competitor of the magnetic order we find a intra-moiré-cell charge-modulated state for which the Coulomb repulsion in the AB-regions between the AA-spots flows to strong coupling. This instability has not been discussed before and deserves further studies. In particular one should understand whether this state can help explain the phenomenology of the twisted-bilayer graphene systems, e.g. by depleting the low-energy density of states between the AA regions even further and therefore localizing electrons in nearly isolated AA islands. Furthermore, it might well be that, once the charge redistribution has been accounted for, e.g. by self-energy terms, the previous magnetic instabilities become relevant again beyond the parameter range found in this work.
The main goal of this work was to demonstrate that one can set up functional renormalization group calculations within useful approximations even for such systems with O(10 4 ) bands. We plan to extend the method to allow for more quantitative comparisons and also to treat other pairing bilinears in the P channel and momentum dependencies in the other channels. This will allow us to study, e.g., unconventional superconductivity and other bond ordering phenomena on the moiré scale in the IOBI fRG for systems like 'magic angle' twisted bilayer graphene.

VI. ACKNOWLEDGEMENTS
The German Science Foundation (DFG) is acknowledged for support through RTG 1995. Most simulations were performed with computing resources granted by RWTH Aachen University under project rwth0496. In addition, the authors also gratefully acknowledge the computing time granted through JARA on the supercomputer JURECA at Forschungszentrum Jülich [58].