Spontaneous symmetry breaking in a honeycomb lattice subject to a periodic potential

Motivated by recent developments in twisted bilayer graphene moir\'e superlattices, we investigate the effects of electron-electron interactions in a honeycomb lattice with an applied periodic potential using a finite-temperature Wilson-Fisher momentum shell renormalization group (RG) approach. We start with a low-energy effective theory for such a system, at first giving a general discussion of the most general case in which no point group symmetry is preserved by the applied potential, and then focusing on the special case in which the potential preserves a $D_3$ point group symmetry. As in similar studies of bilayer graphene, we find that, while the coupling constants describing the interactions diverge at or below a certain"critical temperature"$T=T_c$, it turns out that {\it ratios} of these constants remain finite and in fact provide information about what types of orders the system is becoming unstable to. However, in contrast to these previous studies, we only find isolated fixed rays, indicating that these orders are likely unstable to perturbations to the coupling constants. Our RG analysis leads to the qualitative conclusion that the emergent interaction-induced symmetry-breaking phases in this model system, and perhaps therefore by extension in twisted bilayer graphene, are generically unstable and fragile, and may thus manifest strong sample dependence.


I. INTRODUCTION
Moiré superlattices of various kinds have recently become a topic of great theoretical and experimental interest, spurred in large part by the discovery of superconductivity at a surprisingly high transition temperature 1,2 in twisted bilayer graphene (tBLG) with the twist angle close to some small "magic angle". In addition, other experiments 3,4 have uncovered a Mott-like correlated insulator phase in tBLG, and various experimental and theoretical investigations of this and other, similar, systems have been undertaken. For example, a gate-tunable correlated insulator phase was also found in trilayer graphene moiré superlattices 5 , and a crystal field-induced gap in encapsulated twisted double bilayer graphene has been observed 6 . A consensus that electronelectron interaction, enhanced strongly by the flatband nature of the moiré system, is the driving mechanism producing various symmetry-breaking phases in twisted bilayer graphene has emerged although the possibility that superconductivity itself may arise from electron-phonon interactions cannot be ruled out [7][8][9][10] . One work proposes a number of potential correlated insulating phases in tBLG 11 . Another, also considering tBLG, considers the interplay between van Hove singularities and the symmetry-breaking phases, and proposes different types of s-wave superconductivity 12 . Yet another 13 considers tBLG at two different fillings, 1 and 2 electrons per moiré unit cell, finding a ferromagnetic stripe phase in the former case and an approximately SU (4) symmetric insulating state in the latter. Two other works consider other moiré superlattice systems, with one 14 describing the existence of Chern bands in twisted double bi-and trilayer graphene and in hexagonal boron nitride in the presence of an applied electric field, while another 15 investigates a number of models of generic moiré systems. The sub-ject is exhaustive with hundreds of theoretical papers proposing different mechanisms for symmetry-breaking ground states using many different approximations and models, and we do not attempt a review of the subject, only mentioning a few representative recent publications where renormalization group (RG) techniques were used to study the moiré ground states [16][17][18][19][20][21][22] .
We will theoretically investigate, using extensive RG techniques, a related system to those described abovea honeycomb lattice, such as that formed by graphene, subject to a periodic potential. While we expect differences with tBLG, we expect that similar physics should arise generically in both systems. Our goal is to obtain universal effects of interactions in the most general situation, necessitating leaving out many important, but nonessential, realistic complications in the experimental moiré systems such as strain, substrates, disorder, band structure details. Such a generic theory using the minimal model of our work could be the starting point for more realistic future theories. We therefore restrict ourselves to the seemingly simple model of a honeycomb lattice subjected to an external periodic potential to mimic the essential features of twisted bilayer graphene. As we would see, even this simple system is extremely difficult to handle from the perspective of dealing with electronelectron interactions because of the greatly reduced symmetries in the interacting lattice Hamiltonian.
We begin with a low-energy effective theory, which consists of two Dirac cones. In general, the periodic potential could completely eliminate all point group symmetry from the system, leaving only (reduced) translational, time reversal, and spin SU (2) symmetries. In such a case, we may see terms that shift the Dirac cones away from their usual positions at the corners of the Brillouin zone in addition to a mass term that opens a gap (there are experimental indications of such Dirac point gaps in tBLG samples). Our main focus, however, will be on potentials that preserve a D 3 point group symmetry, in which case only the mass term may appear without any shifts in the Dirac cones. We then construct the four-fermion interaction terms that the symmetries of the system allow; we find that 22 such terms are allowed purely by symmetry, and that this number may be reduced to 10 by the use of Fierz identities. The interaction problem facing us is therefore formidable involving in general a 22-parameter RG flow even within this minimal model of this moire system. Even the reduced problem of a 10-parameter flow has never before attempted in the graphene literature.
We employ a finite-temperature Wilson-Fisher momentum shell RG procedure in this work [23][24][25] . In such a procedure, we rewrite the partition function for our system as a path integral in terms of anticommuting Grassmann fields and impose a momentum cutoff Λ on the electronic modes. Next, we divide these modes into "fast" modes, which are those within a thin shell near the momentum cutoff, and "slow" modes, which are the remaining modes. We then integrate out the "fast" modes perturbatively to one-loop order and rescale the "slow" modes and the various coupling constants to recover an action of the same form as that which we started with. This procedure yields corrections to the various terms in the action, which yield a set of differential equations that we call RG equations describing how the various coupling constants evolve as we integrate out modes. We show that, if we set the temperature to a "critical temperature" T = T c , the constants describing the fourfermion interactions diverge exponentially, but that ratios thereof tend to finite values. Therefore, we find that the coupling constants tend toward "fixed rays" rather than fixed points. These ratios in fact contain information about what symmetry-breaking phases the system is unstable to. To determine which phases these are, we note that, as we integrate out electronic modes, we also generate contributions to the free energy, which is simply given by F = −k B T ln Z, where Z is the partition function. We use this fact to calculate the free energy in the presence of "source terms" corresponding to various symmetry-breaking orders; we may then calculate susceptibilities towards these orders by taking second derivatives with respect to the source term coefficients. We find that, just above the critical temperature, any divergent susceptibility does so as a power of T −T c , with the power related to the coupling constant ratios. If a given susceptibility diverges, then we say that the system is unstable to the associated symmetry-breaking order. Our work is a highly nontrivial (because of the considerably reduced symmetry of the system leading to multi-parameter RG flow) generalization of the earlier RG work on bilayer graphene 26,27 . The imposition of the additional periodic potential considerably complicates the technical aspects of the RG flow compared with these earlier works, and leads to some qualitative differences in the results as discussed below.
We find that, in contrast to similar studies of bilayer graphene 26,27 , the only fixed rays that appear here are isolated. A number of these fixed rays correspond to multiple instabilities. As a result of these two facts, we expect that the results of integrating the RG equations will be sensitive, even qualitatively, to changes in initial conditions. This is not surprising, given the diverse orders found in the literature in tBLG. A necessary conclusion is that the symmetry-breaking phases in the system would be fragile and highly sensitive to all the details of the specific sample, and there could be considerable sample-to-sample variations in the observed phase diagram (or even in the same sample under thermal cycling). We emphasize that this finding of "unstable symmetry breaking" in our minimal model , arising from just pristine interaction effects, could only be much more complex in real tBLG samples where many nonessential realistic effects (e.g. strain, phonon, disorder, substrate) would come into play well beyond our effective low-energy RG theory. The rest of this work is organized as follows. In Sec. II, we introduce the system and our low-energy effective theory. We then describe and implement our RG procedure in Sec. III. In Sec. IV, we describe how we obtain the fixed rays, and then describe how we determine what instabilities they correspond to in Sec. V. We give our conclusions in Sec. VI, and provide further technical details of the calculation in the appendices.

II. SYSTEM AND MODEL
We consider here electrons on a tight-binding honeycomb lattice subject to a periodic potential. Such a lattice, in the absence of the periodic potential, possesses a D 6 point group symmetry, along with time reversal, spin SU (2), and discrete translational symmetries. In the "worst-case" scenario, the potential removes all point group symmetries of the honeycomb lattice, leaving only time reversal, spin SU (2), and translation symmetries (the last of these reduced by the applied potential). We assume that the applied potential is commensurate with the honeycomb lattice. If the potential is arranged in such a way as to place a maximum or minimum at a lattice site, however, then the full system will possess a D 3 point group symmetry; we will in fact focus on this case after a brief discussion of the case in which we have no point group symmetry at all.
We will adopt a low-energy effective theory of this system, which consists of two Dirac cones (valleys ± K), a sublattice degree of freedom (A/B), and spin (↑/↓). We also include all "mass" terms allowed by the symmetries of the system. We list the valley and sublattice components of all bilinears that may be formed, along with what representations they transform under both with and without the D 6 point group symmetry, in Table I. Some of the bilinears transform nontrivially under translations; these instead transform under representations of the D 3 "small group" of the wave vector, K. For the sake of a self-contained presentation, we also list the representations 28 of D 6 and D 3 and their associated characters in Tables II and III, respectively. Here, each matrix is of the form τ i σ j s k , where the τ matrix operates on the valley pseudospin, the σ matrix on the sublattice pseudospin, and the s matrix on the actual spin. We should note here that the form of the low-energy theory is completely independent of the exact details of the applied potential, though determining the values of the constants appearing therein requires a more detailed calculation.
The model Hamiltonian for the case with no point group symmetry is similar in form to the low-energy effective theory for a honeycomb lattice, but includes additional terms that transform trivially, i.e., under the A+ representation, under the remaining, reduced, symmetry group of the system with an applied periodic potential. This Hamiltonian is (1) The first two "mass" terms, the 1σ x 1 and τ z σ y 1 terms, simply displace the Dirac cones away from ± K, while the third, the 1σ z 1 term, opens a gap in the cones. We have so far covered the noninteracting Hamiltonian; we now turn our attention to the interaction terms. All such terms must, as with the noninteracting terms, be invariant with respect to the symmetries of the system. We may form such invariant terms by use of the generalized Unsöld theorem; all of these will have the form, where the matrices S i,j both belong to the same "row" of a given representation 28 . In the D 6 symmetric case, there is only one matrix per "row" for each representation; for example, is the only interaction term that we can associate with the E 2 + representation. On the other hand, in the case of no point group symmetry, we see that there are three matrices in the same "row" of, say, the A+ representation. There are therefore six distinct ways to form interaction terms within this representation. Overall, in the case with no point group symmetry, we find that there are 54 allowed interaction terms, while, in the case of D 3 point group symmetry, the number is reduced to 22. These numbers may be reduced further through the use of Fierz identities, which are, for 8 × 8 matrices, of the form, where the Λ n are all possible matrices of the form, τ i σ j s k , and we omitted the position dependence of the operators ψ for the sake of brevity (we assume that all are at the same position). As may be seen, these identities relate the various interaction terms to one another. Using these identities, we may reduce the number of independent couplings to 22 in the case of no point group symmetry and 10 in the case of D 3 point group symmetry. Even though we will not consider the D 6 case here, we note that such a symmetry allows 18 interaction terms, which may be reduced to 9 using Fierz identities.

III. RENORMALIZATION GROUP (RG) PROCEDURE
We now turn our attention to describing the Wilson-Fisher momentum shell RG technique that we will employ 24 . This technique is as follows. We begin by writing down the partition function for our system as a path integral: where the ψ are now Grassmann numbers corresponding to the coherent states of the corresponding operators and the action S is given by where H is the interacting Hamiltonian written in "normal order" (i.e., Hermitian conjugates to the left of nonconjugates) and τ is an "imaginary time". The next step is to integrate out electronic modes in momentum shells. We divide the fields into "slow" modes, denoted by ψ < , and "fast" modes, denoted by ψ > , where the fast modes are those modes with momenta within the shell, Λe −δ ≤ k ≤ Λ, and δ is a small change in a scale parameter used to characterize how many modes have been integrated out. Finally, we rescale all momenta of the remaining, slow, fields to k = ke δ to restore the region of integration over momentum to k ≤ Λ, and then rescale the fields and constants, thus recovering the original overall form of the action. This procedure may be done exactly for the noninteracting system, but must be done perturbatively once we introduce interactions. We may express this renormalization of constants in the form of differential equations (RG equations), as we will see shortly.
We now apply this procedure to the system under consideration and, from this point on, we will specialize to the case of D 3 point group symmetry for the A1+  (2), time reversal, and translations), those of the reduced symmetry group with only D3 point group symmetries, and those with no point group symmetries. The + or − in each representation name indicates how the bilinears transform under time reversal (even or odd, respectively). Corresponding to each of these "charge" representations are "spin" representations, in which the spin part of the matrix is s k , k = x, y, z; these transform oppositely with respect to time reversal to the "charge" representations (e.g., the A1 spin representation is odd under time reversal). , along with how they transform under time reversal (A "+" means the representation is even, while "−" means that it is odd).
Representation E 2C3 C 2 TABLE III: Representations of the group D3, along with how they transform under time reversal (A "+" means the representation is even, while "−" means that it is odd).
sake of relative simplicity. In this case, v F x = v F y and k 0x = k 0y = 0, as the associated terms no longer transform trivially under the symmetries of the system. We first determine how the various constants determining the theory rescale at "tree" level, i.e., at lowest order in the perturbative expansion. Performing a Fourier transform, we find that the noninteracting part of the action S 0 is where kω is a shorthand for the Matsubara sum (here, the sum is over fermionic Matsubara frequencies, ω n = (2n+1) π β , where n is an integer) and momentum integral, The interaction terms, on the other hand, are given by where 1234 is a shorthand for the integrals and sums appearing in this expression along with delta functions expressing momentum and energy conservation, and ψ(n) = ψ( k n , ω n ).
If we now perform the procedure summarized above, we find that the various constants in our theory rescale at tree level as follows: The last may be rewritten in terms of temperature as T = T e δ . These may also be recast as differential equations; letting x = x( + δ ) and x = x( ), where x is any one of the above constants, we may easily show that We now turn to one-loop corrections. These are the highest-order corrections that will appear within our RG analysis, as multi-loop corrections will be of order (δ ) k , k > 1, and thus will vanish in the resulting RG equations. We obtain contributions to m, µ, and g SU at this order; we show the diagrams corresponding to the m and µ corrections in Fig. 1 and those for the corrections to g i in Fig. 2. Before we begin determining these corrections, we need the bare Green's function for the system G 0 ( k, ω); it is given by A. One-loop corrections to mass and chemical potential We will begin with the contributions to m and µ, and with the "tadpole" diagram, Fig. 1a. Both of these contributions are first order in the interaction terms, and come from those terms containing two "slow" and two "fast" modes. This corresponds to terms of the form, Using Wick's theorem and the fact that G 0 ( k, ω) = ψ( k, ω)ψ † ( k, ω) , we find that This simplifies to where the notation, > k ω , simply means to integrate only over the "fast" momenta. These sums and integrals may easily be evaluated; we obtain and where and If we now substitute this into the expression for ∆S(tadpole), we find that We see that this represents a correction to one of either the chemical potential term or the mass term, depending on whether S i,1 = 1 8 or S i,1 = 1σ z 1. In fact, the traces in this expression will only be nonzero if S i,2 = 1 8 or S i,2 = 1σ z 1, meaning that only those terms corresponding to the A 1 + representation give a nonzero contribution via the "tadpole" diagram. We now consider the "sunrise" diagram, Fig. 1b, which corresponds to terms of the form, Diagrams representing the one-loop corrections to the four-fermion couplings gi. The solid red lines represent "fast" modes, the solid black lines "slow" modes, and the dashed lines interactions.
Evaluating the averages as before, we get Using the formulas given earlier for the integrals over the "fast" momenta, this becomes Unlike the "tadpole" contribution, the "sunrise" contribution will yield nonzero terms proportional to the coupling constants coming from all representations, not just A 1 +. The full RG equations for µ and m are of the form, When we evaluate the coefficients, however, we find that all of the B z µ,i = 0 and that all of the B 0 m,i = 0. This simplifies our later calculations, as we will see below. This also implies, as may be seen from the form of the equations, that, if we set either µ or m to zero, then we will never generate them, i.e., they will not renormalize to nonzero values.

B. One-loop corrections to four-fermion coupling constants
We next determine the corrections to the four-fermion coupling constants g SU , which are depicted in Fig. 2.
Evaluating these five contributions, we find that the RG equations for the four-fermion couplings all have the form, In arriving at this form, we evaluate integrals and sums of the form, We present the results of doing so, in the form of the functions Φ a , in Appendix A. We list the expressions obtained from each of the diagrams in Fig. 2 in Appendix B.

IV. FIXED RAYS
We now determine what the various outcomes of integrating the RG equations are. We start by showing that, if the temperature is tuned to what we will call the critical temperature, T = T c , integrating the RG equations for the four-fermion couplings g i will result in at least one of said couplings diverging exponentially. However, it turns out that ratios of these couplings tend toward fixed values; we will show later that these "fixed rays" in the space of the g i tell us what symmetry-breaking orders the system is unstable to. If the temperature is above this critical temperature, then the g i all saturate to finite values as → ∞. If, on the other hand, the temperature is below the critical temperature, then one or more of the g i will diverge to infinity at some finite value of .
In all of these cases, we need to first solve for the fixed ratios themselves. To do this, we first derive the equations for ratios of the g i with one of the couplings, which we will call g r (which of course is assumed to diverge), Doing this, we obtain We find that the behavior of the equations for large depends on how rapidly T , µ, and m increase, and breaks down into three cases: 1. T runs faster than µ and m.
Which of these three cases we consider determines the form of the equations that we have to solve to determine the fixed rays and what instabilities they represent. We now discuss each case in turn.
Case 1: In this case, the temperature T increases more quickly than µ or m for large . We will determine the large behavior of µ, m, and the g i under this assumption. In the limit of large , the Φ a functions that decrease the most slowly are In addition, the functions F 0 and F z appearing in the equations for µ and m behave as follows.
We first determine the asymptotic behavior of the fourfermion couplings g i for T = T c . The equations will all take the form, whereĀ ijk = A 2,+ ijk + A 2,− ijk . By the definition of T c , we know that, for T = T c , g i ( → ∞) → ∞. We now substitute in g i ( ) = g i,0 e δg , obtaining or (δ g + 1)g i,0 e δg = Λe (2δg−1) 8πT c jkĀ ijk g j,0 g k,0 .
This equation is satisfied if δ g = 2δ g − 1, or δ g = 1. This is consistent with our assertion that g i diverges to infinity as → ∞. We now determine the constants, g i,0 . We mentioned earlier that ratios of any two divergent g i tend to finite values; we will rewrite the above equation in terms of these ratios. If we now rewrite the above equation in terms of these ratios ρ i , we obtain, after simplification, Solving for g r,0 , we obtain where All other g i,0 may be obtained simply by multiplying g r,0 by the appropriate ratio. Next, we consider the equations for µ and m. In the limit of large , these become If we substitute the asymptotic expressions for g i obtained above into these equations, we get where We see that the equations for µ and m reduce to a pair of first-order linear differential equations. Our earlier results imply that B z µ = B 0 m = 0, so that these equations are decoupled. We just need one of either µ or m to increase more slowly than T (or even decrease); if the other increases more quickly, then that means that the corresponding parameter must be set to zero to obtain that outcome.
Case 2: Next, we consider the case in which the g i increases more slowly than µ and m, and in which |µ| > |m|, i.e., the chemical potential is outside the gap. In this case, the most slowly-increasing of the Φ a functions will be We will assume for now that the g i , µ, and m all increase exponentially for large ; we will show here that this assumption is consistent. For concreteness, we assume that g i ( ) = g i,0 e δg , µ( ) = µ 0 e η , and m( ) = m 0 e η . In this case, the equations for g i become Substituting our ansatz for g i into this equation, we get This equation is satisfied if we let δ g = 2δ g − η, or δ g = η.
We may now solve for g i,0 with a similar procedure to the previous case. Rewriting the above equation in terms of the ratios ρ i and simplifying, we get If we now denote the sum over j and k by A r (µ 0 , m 0 ), we get As before, we may obtain the other g i,0 by multiplying the above by the appropriate ratio.
With this result, we can now consider the equations for µ and m. In the limit of large , only F 0 is nonzero: We then obtain If we now make our earlier ansatz and substitute the expression for the g i obtained earlier, we obtain We note, however, that B 0 m,i = 0, so that the second equation simplifies to We thus find that, unless we set m 0 = 0, we must take η = 1, violating our assumption that T increases more slowly than µ or m. If we take m 0 = 0, then the equation for µ becomes However, our expression for A r (µ 0 , m 0 ) reduces to We expect the most likely outcomes of actual integration of the RG equations to fall within the previous case, and thus we will not treat this case further here. In concluding this, we are guided by a similar study of bilayer graphene with a band gap opened by, for example, an applied electric field undertaken in Ref. 26. We also note that the requirement that m 0 = 0 would imply a complete absence of a periodic potential, i.e., we are dealing with just a honeycomb lattice.
Case 3: Finally, we consider the case in which, once again, µ and m increase more rapidly than T for large , but this time |µ| < |m|, i.e., the chemical potential is inside the gap. In this case, the most slowly-decreasing of the Φ a are If we now make the same ansatz as before, taking g i ( ) = g i,0 e δg , µ( ) = µ 0 e η , and m( ) = m 0 e η , then the equations for the g i become Once again, this equation is satisfied if δ g = η. We may solve for this in terms of the fixed ratios as before, obtaining or, denoting the sum on j and k by A r (µ 0 , m 0 ), Now we consider the equations for µ and m. In this case, only F z is nonzero for large : The equations then become If we now make our earlier ansatz and substitute the expression for the g i obtained earlier, we obtain We note, however, that B z µ,i = 0, and thus the above becomes Similarly to the previous case, we conclude that we must set µ 0 = 0 in order to obtain consistency with our assumptions in this case. if we do so, however, we find that A r (µ 0 , m 0 ) simplifies to For similar reasons as in the previous case, we will save further treatment of this case for future work.

V. ANALYSIS OF FIXED RAYS
We now describe how we determine what symmetrybreaking phases each fixed ray represents an instability towards. To do this, we calculate the susceptibility of the system as a function of temperature near the "critical temperature". We can, in turn, do this by noting that, as we integrate out electronic modes in our RG analysis, we produce a multiplicative constant contribution to the partition function that we have been ignoring so far. These multiplicative constants in fact represent contributions to the free energy of the system. Our basic strategy for determining the susceptibilities is as follows. We start by adding "source terms" to the action, which have the form, where the matrices M i run over all possible 8 × 8 matrices of the form, τ i ⊗ σ j ⊗ s k . Note that some of the ∆ may be equal to one another if their associated matrices belong to the same representation of D 3 . These source . 3: Diagrams corresponding to one-loop corrections to the particle-hole (ph) source terms. Solid black lines correspond to "slow" modes, solid red lines to "fast" modes, dashed lines to four-fermion interactions, and wavy lines to a source term vertex.
terms correspond to different "particle-hole" (ph), or excitonic, and "particle-particle" (pp), or superconducting, order parameters. We provide a list of the representations that each of these source terms correspond to, along with what order they represent, in Tables IV (ph) and V (pp). We note that the excitonic states are very similar to those that are possible in bilayer graphene 27 due to the mathematical similarity to that case, though the physical interpretation will be slightly different. In the pp case, only those terms with antisymmetric matrices M i appear. We determine the contributions to the free energy as we integrate out modes up to second order in these source terms. Finally, we can calculate the susceptibilities to various order parameters by calculating second derivatives of the free energy: where f is the free energy per unit area.

A. One-loop RG equations for the source terms
Before we determine the free energy, we must also determine how these source terms renormalize to one loop. At tree level, the renormalized source terms are given by ∆( ) = ∆ 0 e , or We now determine the one-loop corrections, depicted in Figs. 3 and 4. The corrections for the ph source terms are given in Fig. 3. The first diagram, Fig. 3a, yields 4: Diagram corresponding to one-loop corrections to the particle-particle (pp) source terms. Solid black lines correspond to "slow" modes, solid red lines to "fast" modes, dashed lines to four-fermion interactions, and wavy lines to a source term vertex.
The second, Fig. 3b, yields where Now we consider the pp source terms. In this case, we have only one diagram contributing to one-loop renormalization, shown in Fig. 4. This diagram yields Overall, these contributions lead to RG equations of the form, We now consider the behavior of the source terms for large and at T = T c . For reasons stated earlier, we will focus only on the case where T increases  List of representations of D3 and the corresponding particle-hole (excitonic) order parameters. The ± after each representation represents how the corresponding charge order transforms under time reversal. Only the valley and sublattice components of the matrices are shown. We list both the charge and spin variants of each order in the same row; the spin order has the opposite time reversal symmetry to the corresponding charge order (e.g., the ferrimagnetic state is odd under time reversal).

Representation
Matrices Order List of representations of D3 and the corresponding particle-particle (superconducting) order parameters. The letter after each representation name denotes whether the order is a singlet (s) or triplet (t) order. We omit the spin matrix in the list of matrices; it is sy for singlet orders and 1, sx, or sz for triplet orders.
more rapidly than µ and m. As stated before, the Φ a functions that decrease the most slowly in this case are Φ 2,+ ≈ Φ 2,− ≈ Λ 8πT . The equations for the ∆ i then become, after substituting the large expression for the g i , We will thus have, at most, a system of two linear equations describing a given source term. If a given term does not couple to another, then it will simply be given by ∆ i ( ) = ∆ i ( 0 )e ηi , where Otherwise, we simply solve the system of equations using standard techniques.

B. Free energy
Now that we have derived the RG equations, we turn our attention to the free energy. More specifically, we will FIG. 5: Diagrams representing one-loop contributions to the free energy from the (a) particle-hole source terms and (b) particle-particle source terms. The wavy lines represent the source term vertices and the red lines to fast electronic modes.
calculate the contribution to the free energy per unit area from the source terms alone. The diagrams that represent contributions from the source terms is shown in Fig.  5. The first diagram, Fig. 5a, represents the contribution from the ph source terms, and yields the following result for the free energy per unit area from the source terms: The second diagram, Fig. 5b, represents the contribution from the pp terms, and yields With these results, we may now derive the susceptibilities and thus determine which of them diverge for a given fixed ray. More specifically, we will determine their behavior for T close to, and just above, T c . We start by revisiting the equations for the g i . In this case, we may still use the large approximations for the Φ a , but now we assume that g i tends to a very large, but finite, value as → ∞. If we solve the equation for the g r that we divide by to obtain the ratios ρ i in this case, we get where 0 1. We now make use of the fact that, at to rewrite the above as We may now use the fact that, to first order in T − T c , where c r is a constant. We now perform a similar analysis of the equations for the source terms. Doing so, we find that, for a ∆ i that is not coupled to any other ∆ j , where C r is a constant proportional to the constant c r appearing in the equation for g r ( , T ).
If we now substitute these results into the free energy and determine the contribution from the integral on for > 0 (which is where we expect the divergence to come from), we find that it goes as (T − T c ) 2−ηi , and thus the susceptibility, will diverge with the same exponent provided that η i > 2. We note that ∆ i ( = 0 ) ∝ ∆ i ( = 0) simply due to the linearity of the equations giving ∆ i . Therefore, if the condition that η i > 2, or is satisfied, then the corresponding susceptibility exponent diverges, and we claim that the system is unstable to the associated order. A similar analysis for the case of coupled source terms shows that, provided that the system of RG equations yields at least one solution e ηi that satisfies the condition that η i > 2, then we have an instability towards the corresponding order. We note that coupling of source terms only occurs if they correspond to the same representation of the symmetry group of the system.
With these results, we are now ready to determine the leading instabilities of the system that the various fixed rays correspond to. Unlike in the case of bilayer graphene, considered in Refs. 26 and 27, we do not find any one-parameter or more families of fixed rays; all of the fixed rays are isolated. However, we find a very large number (thousands) of solutions. A number of these rays correspond to multiple instabilities simultaneously present. The fact that we only obtain isolated fixed rays, rather than any multiparameter families, indicates that the system is unstable to perturbations in the initial conditions, causing the system to converge to a different fixed ray. These two facts are not surprising, given the diversity of symmetry-breaking states found so far in the related, but not identical, twisted bilayer graphene system.

VI. CONCLUSION
We have investigated the possibility of instabilities to interaction-induced symmetry-breaking orders in a honeycomb lattice away from half-filling, i.e., the chemical potential µ = 0, subject to a periodic potential. For simplicity, we assumed that the periodic potential preserved a D 3 point group symmetry, along with translational symmetry (though reduced from that of the unmodified lattice), time reversal, and spin SU (2). This allows for a mass gap m to be present. We employ a finitetemperature Wilson-Fisher momentum shell RG procedure in this work. We derived the RG equations for the four-fermion coupling constants g i , the chemical potential µ, the mass m, and the temperature T . Also for simplicity, we focused on the case in which T increases more quickly at large RG scaling parameter than µ or m. We then showed that, at some "critical temperature" T c , the coupling constants diverge exponentially, but that ratios thereof remained finite. We finally showed how these ratios could be used to determine what symmetry-breaking orders the system would be unstable to. We found that there were thousands of isolated fixed rays, a contrast to similar studies of bilayer graphene 26,27 , where a twoparameter family of fixed rays was found in addition to only a few isolated rays. In some cases, these rays corresponded to instabilities toward several different orders, which is not surprising given the diverse orders detected in experiments on the related, though not entirely identical, twisted bilayer graphene.
Our detailed multi-parameter RG analysis within a simple minimal lattice model points to the real possibility that moiré superlattices (e.g. twisted bilayer graphene near the magic angle), may manifest "unstable symmetry-breaking" where the symmetry-breaking phases are intrinsically fragile, and the physics depends sensitively on all the details and initial conditions, even excluding the realistic complications of disorder, strain, phonon, substrate, etc. Our work is consistent with there being considerable sample dependence in the observed phenomenology of various correlated exotic phases in tBLG. Such a generic fragile "unstable symmetry breaking" scenario leads us to conclude that it is likely that experimental development would lead to the observation of many exotic ground states in the system. Similar considerations in these previous works on bilayer graphene concerning the exponents apply here as well; we expect that our procedure captures the basic qualitative behavior of the susceptibilities (i.e., whether or not they diverge), even if the exact exponents are not quite correct. We also note that many of these fixed rays correspond to multiple instabilities. Our method does not give further information other than the possibility of these orders emerging. Other methods are required to determine which of these orders actually emerges.