Extended Hubbard model in undoped and doped monolayer and bilayer graphene: Selection rules and organizing principle among competing orders

Performing a leading-order renormalization group analysis, here we compute the effects of generic local or short-range electronic interactions in monolayer and Bernal bilayer graphene. Respectively in these two systems chiral quasiparticles display linear and biquadratic band touching, leading to linearly vanishing and constant DOS. Consequently, the former system remains stable for weak enough local interactions, and supports a variety of ordered phases only beyond a critical strength of interactions. By contrast, ordered phases can nucleate for sufficiently weak interactions in bilayer graphene. By tuning the strength of all symmetry allowed local interactions, we construct various cuts of the phase diagram at zero and finite temperature and chemical doping. Typically, at zero doping insulating phases (such as charge-density-wave, antiferromagnet, quantum anomalous and spin Hall insulators) prevail at the lowest temperature, while gapless nematic or smectic liquids stabilize at higher temperatures. On the other hand, at finite doping the lowest temperature ordered phase is occupied by a superconductor. Besides anchoring such an organizing principle among the candidate ordered phases, we also establish a selection rule between them and the interaction channel responsible for the breakdown of linear or biquadrtic chiral nodal Fermi liquid. In addition, we also demonstrate the role of the normal state band structure in selecting the pattern of symmetry breaking from a soup of preselected incipient competing orders. As a direct consequence of the selection rule, while an antiferromagnetic phase develops in undoped monolayer and bilayer graphene, the linear (biquadratic) band dispersion favors condensation of a spin-singlet nematic (translational symmetry breaking Kekul\'e) superconductor in doped monolayer (bilayer) graphene, when the on site Hubbard repulsion dominates in these systems.


I. INTRODUCTION
Multiband electronic materials constitute a rich landscape harboring a plethora of competing phases, among which spin-and charge density waves, nematicity, unconventional superconductivity are the most prominent and commonly occurring ones, leaving aside the mysterious quantum spin liquids. Indeed, the global phase diagram of strongly correlated electronic materials, such as cuprates, pnictides, iridates, ruthenates, and heavyfermion compounds, displays a confluence of at least some of these phases. In experiments, these phases can be realized as the temperature, pressure, and chemical doping is tuned in the system. However, owing to the complexity of the band structure, it is often challenging to construct a unifying minimal model for these correlated materials that captures the salient features of their global phase diagram. In this regard emergence of a new class of materials, semimetals, featuring band touching at isolated points in the Brillouin zone (BZ), opened a new frontier. Even though semimetals are fascinating phases of matter from the perspective of their topological properties [1][2][3], here we focus to address the role of * szabo@pks.mpg.de † bitan.roy@lehigh.edu strong electronic interactions in such gapless electronic materials. Specifically, we consider two paradigmatic representatives of two-dimensional semimetallic systems, namely monolayer graphene (MLG) and Bernal-stacked bilayer graphene (BLG), where noninteracting itinerant electrons display linear and biquadratic band touchings at the corners of the hexagonal BZ, respectively [4,5].
In what follows, restricting ourselves to local or momentum-independent repulsive four-fermions interactions on monolayer and Bernal bilayer honeycomb membrane of carbon atoms, we accomplish the following. (1) We establish a correspondence between a given local fourfermion interaction that ultimately causes the breakdown of either the linear or the biquadratic nodal chiral Fermi liquid and the broken symmetry phase(s) it supports in the ordered state. Throughout this paper, such a correspondence is referred to as the selection rule. (2) In addition, we demonstrate an organizing principle among competing phases as the temperature or chemical doping is tuned in the system, based on a generalized energyentropy argument. And finally, (3) we highlight the key role played by the normal state band structure in selecting an ordered phase from a soup of incipient competing orders, otherwise chosen by the selection rule. The whole construction is summarized in Fig. 1. Specifically, in this regard, MLG and BLG provide a unique opportunity as these two systems are indistinguishable from the point of view of the point or space group symme-FIG. 1. A schematic flowchart depicting the selection rules and organizing principle. Consider a noninteracting system, described by the single-particle Hamiltonian H0 containing a set of Hermitian matrices {Hi}, a four-fermion interaction with coupling g M and a fermionic order parameter coupled with the conjugate field ∆O, all determined by the symmetry of the system. Here Mi and Oi are Hermitian matrices of the same dimension as Hi. Then, an ordered phase is boosted by the interaction if (a) Oi ≡ Mi or (b) these two set of matrices maximally anticommute. The selected ordered phases are then arranged along the temperature axis according to the organizing principle, quantified by AH and CH , respectively counting the pairs of matrices from the sets {Oi} and {Hi} that mutually anticommute and commute. Order parameter matrices Oi that maximally anticommute (commute) with Hi yield gapped (gapless) quasiparticles inside the ordered phase and inhabit comparatively lower (higher) temperatures, following the energy-entropy argument. Finite chemical doping promotes superconductivity at the lowest temperature even from repulsive electronic interactions, as the underlying Fermi surface can only be gapped by them, also in agreement with the organizing principle. The nature of the nucleated paired state follows the selection rule. try, although they display distinct normal state band structures, see Fig. 2. It is worthwhile mentioning that recently the author(s) proposed a set of selection rules and organizing principle among competing orders for (a) only the on-site Hubbard repulsion in a two-dimensional slow Dirac liquid in twisted bilayer graphene in Ref. [6] and (b) generic repulsive quartic interactions for threedimensional strong spin-orbit coupled biquadratic effective spin-3/2 Luttinger fermions in Ref. [7] that can be relevant for 227 pyrochlore iridates and half-Heuslers, for example. The findings reported in this work are com-pletely in agreement with the ones in these two previous works [6,7]. Therefore, the present discussion besides providing new insights into the global phase diagram of the extended honeycomb Hubbard model in (un)doped MLG and BLG that remained unnoticed despite the existing vast literature on this topic , more importantly extends the jurisdiction of the existing selection rule and the organizing principle to a completely different class of semimetallic compounds. Altogether these analyses suggest a general applicability of the proposed selection rule and organizing principle among competing orders in generic correlated multiband systems, which can be predicted a priory from the symmetry of the system, the noninteracting band structure Hamiltonian, and the (anti)commutation relations among matrices characterizing the four-fermion interaction channels and the order parameters, see Fig. 1. Subsequently, these predictions can be anchored by performing an unbiased RG analysis, as we show here and also in Refs. [6,7].
To briefly and schematically summarize our main findings, here we focus on generic short-range or local, but repulsive four-fermion interactions and construct various cuts of the phase diagram for MLG and BLG at zero and finite temperature and chemical doping by performing an unbiased leading-order renormalization group (RG) analysis in the spirit of appropriate expansions. Typically, these cuts show excitonic orders at zero chemical doping, among which charge density wave, antiferromagnet, quantum anomalous and spin Hall insulators, and various nematic and smectic liquids are the dominant ones. Often insulators (such as, antiferromagnet) are realized at the lowest temperature, while their nodal counterparts (such as, nematic and smectic orders) nucleate at a relatively higher temperature. On the other hand, as the chemical doping is tuned away from the band touching points some pairing order sets in at the lowest temperature that maximally gaps the underlying Fermi surface, following the spirit of the Kohn-Luttinger mechanism [65,66].
To relate our findings to a more experimentally feasible scenario, we construct the phase diagram of (un)doped MLG and BLG, when electrons therein interact with each other only through on-site Hubbard repulsion. Results are shown in Figs. 3(a) and 3(b), respectively. At zero doping, both systems support antiferromagnetic ordering. However, at finite doping while a spin-singlet nematic pairing develops in MLG at the lowest temperature, in doped BLG, the paired state represents a translational symmetry breaking spin-singlet Kekulé superconductor. It is an example of a commensurate Fulde-Ferrell-Larkin-Ovchinikov pairing [67,68], also known as pair-density-wave. Our predictions can be experimentally tested on cold atomic systems, for example, where both MLG and BLG have been engineered, and both doping and the strength of the Hubbard repulsion can be tuned efficiently [69][70][71][72]. In addition, electronic BLG supports various ordered phases near the half filling that include both insulating and gapless (possibly nematic) The b sites (belonging to sublattice B) are generated by δ1, δ2 and δ3. (b) The BZ of MLG, a hexagonal structure rotated by 90 degrees compared to the real space lattice. Six Dirac points are located at its corners, of which two sets are inequivalent, marked with blue and red circles. The representative valleys are chosen to be at K and −K, respectively. (c) Bernal stacked honeycomb bilayer, with the subscript denoting the layer index i = 1, 2 of the sites. Each a1 and b2 sites overlap. The eigenstates of the high energy split-off bands reside dominantly on these overlapping dimer sites. (d) After integrating out the a1 and b2 sites, the resulting structure is another honeycomb lattice, demonstrating the identical point group symmetry of MLG and BLG.
states [73][74][75][76][77][78]. Therefore, increasing the carrier density by applying an external gate voltage one can, at least in principle, induce superconductivity in such systems, where our predictions can be tested directly.

II. EXTENDED SUMMARY
In this paper, we examine the effects of short range electronic interactions in MLG and BLG at zero, as well as finite temperature and chemical potential. Our aim in this undertaking is two-fold. On the one hand, we present a standalone perturbative analysis within an unbiased leading-order Wilsonian momentum-shell RG framework. On the other hand, we apply a set of selection rules, proposed by the authors in Ref. [7], that organize the nucleating phases along various axes of the phase dia- The white region represents disordered chiral Fermi liquid or semimetal at finite or zero doping, respectively, and the gray region depicts an ordered state. Along the phase boundary between them, the system supports an antiferromagnet at and near, and a superconductor (SC) away from the half-filling (µ = 0). Nature of the SC is mentioned in each panel. Realizations of distinct paired states in MLG and BLG solely stem from the distinct normal state band structures in these two systems. Inside the ordered state, one expects a regime of coexistence between antiferromagnet and the adjacent SC, due to their internal O(5) symmetry, with pure phases on either side of it. Our RG analysis, however, cannot capture such coexistence, which nevertheless can be demonstrated from a standard mean-field analysis [6]. Here, µ > 0 (µ < 0) correspond to electron (hole) doped systems.
gram, e.g., the interaction strength, temperature, and chemical doping, and thereby further demonstrate the versatility and applicability of these rules in general correlated multiband electronic systems. The different band structures but identical underlying point group symmetries of MLG and BLG makes the study of these two systems in tandem specifically suited for this purpose, about which more in Sec. II A. In this section we give an outline of our findings, including the quintessential technical details. Interested readers are invited to consult Secs. III-VI for more detailed discussion. Our analysis is built on simple tight binding models of MLG and BLG, consisting of a single and double layer of carbon-based honeycomb lattice, respectively. In MLG [see Fig. 2(a)], we account for the nearest neighbor (NN) hopping only. The BZ of MLG is itself a honeycomb structure, rotated by 90 degrees compared to the realspace lattice, see Fig. 2(b). The NN hopping results in six Dirac cones located at the six corners of the BZ, only two of which are inequivalent, while the rest are related to them by reciprocal lattice vectors. In Bernal stacked BLG [ Fig. 2(c)], besides the intralayer NN hopping in each layer, we couple the two honeycomb sheets via direct interlayer hopping terms. The associated BZ is analogous to that of MLG, but the band structure now displays biquadratic touchings at the six corners of the BZ, again only two of which are inequivalent. Furthermore, the additional layer degree of freedom gives rise to two gapped bands that are localized to the overlapping dimer sites. These so called split-off bands can be integrated out to obtain an effective description of the low lying biquadratic band touching. The resulting effective lattice is another honeycomb lattice, see Fig. 2(d). Consequently, MLG and BLG, possessing identical D 3d point group symmetry, are described by equal number of lowenergy degrees of freedom, while their distinct normal state band structures allow us to investigate its role in the nature of the ordered phases.
Our main object of interest is the low-energy continuum description of interacting spin-1/2 electrons in MLG and BLG, captured by their respective action S = S 0 + S int , assuming the schematic form with r being position, τ imaginary time, d the spatial dimensionality, and Ψ and Ψ † are independent conjugate Grassmann variables. The chemical potential term µ is set to zero unless otherwise mentioned. The operatorĥ(k) in the momentum space describes the band structure, which in the continuum formalism consists of two copies of a linear (quadratic) band touching points in MLG (BLG) (see Sec. IV A), after accounting for the two inequivalent valleys. Local or momentum-independent, electronic interactions are described by the quartic terms in S int , where Γ µν are Hermitian matrices (with µ, ν = 0, 1, 2, 3), g µνλρ are coupling constants and the summation over repeated indices is assumed. Power counting provides the scaling dimension of the coupling constants to be [g µνρλ ] = z − d, where the dynamical critical exponent z = 1 (2) for MLG (BLG). Therefore, sufficiently weak contact interactions are irrelevant in MLG in the RG sense. However, the corresponding linearly dispersing fermions undergo a continuous quantum phase transition, when the interactions become sufficiently strong. In comparison, local interactions are marginal in BLG and the quadratically dispersing chiral quasiparticles display "BCS-like" instabilities in the presence of infinitesimal interactions.
We construct the interacting Lagrangian containing all symmetry-allowed local or intra-unit cell four-fermion terms. We show that momentum-independent electronic interactions are described by only 9 linearly independent quartic terms. See Appendix B. On the other hand, spontaneous symmetry breaking is described via 27 local order parameters. Of them, 18 are excitonic, which are grouped into 9 spin singlet and 9 spin triplet orders, see Table I. The pairing sector accommodates in total 9 local pairings, out of which 5 are spin singlets and the remaining 4 are spin triplets, see Table II.
Following the spirit of the expansion around the lower critical dimension d c = z, with = d − z, we carry out a one-loop Wilsonian RG analysis in the individual interaction channels in MLG and BLG. We identify the dominant instabilities, thereby constructing various cuts of the global phase diagram, that are displayed in Figs. 7, 8, 9 and 10. We claim that the selection of ordered phases for a given dominant interaction channel is not arbitrary, rather in general obeys certain algebraic principles, which we discuss in Sec. II A. See also Fig. 1.
Finally, in Sec. II B, we exemplify our findings by focusing on the microscopic extended honeycomb Hubbard model on MLG and BLG, and carry out analogous analyses in three distinct interaction channels, the on-site (U ), NN (V 1 ), and next-nearest-neighbor (V 2 ) components of the Coulomb repulsion. Instructive cuts of the phase diagram of this model are shown in Figs. 4-6.

A. Selection rules
In Ref. [7], the authors introduced a set of directives that organize possible broken symmetry phases in the global phase diagram of interacting (effective) spin-3/2 fermions in a three-dimensional (3D) Luttinger (semi)metal [79], displaying bi-quadratic band touching in the normal state. 1 These consist of (I) a set of selection rules that allow only two types of fermion bilinear order parameters to be boosted by a certain four-fermion interaction channel, and (II) a generalized energy-entropy argument that organizes these ordered phases along the temperature axis. We here demonstrate the general applicability of these rules by verifying them in a vastly different system. While the Luttinger semimetal is realized in 3D strong spin-orbit coupled systems, MLG and BLG are constituted by two-dimensional honeycomb membranes of carbon atoms. In these systems, due to the low atomic number of carbon, the spin-orbit coupling is negligibly small. In fact we completely neglect spin-orbit coupling throughout this paper. Here we quote these two rules for the sake of completeness.
Let the band structure be described by the Hamiltonianĥ, containing Hermitian matrices H i . Consider a local interaction channel and an arbitrary but symmetry allowed local order parameter schematically written as respectively. Here g M is the coupling constant, and ∆ O is the conjugate field to the fermion bilinears. Moreover, we define the following quantities. Let  6) and (7), in the presence of only the on-site repulsion (U > 0). Each cut is shown in the (U, t) plane, where t is dimensionless temperature and U is the strength of on-site repulsive interaction, for various fixed values of dimensionless chemical potential µ [see Eq. (33)]. The shaded region corresponds to the ordered phase, while colored lines indicate phase boundaries of second order phase transitions from the disordered phase (white region), which at finite doping represents chiral Fermi liquids. Along such phase boundaries the blue (red) lines represent excitonic (superconducting) orders, whose nature we explicitly highlight in each panel. The effects of only U > 0 are displayed in (a) and (b) for MLG and BLG, respectively. Here U parameterizes the bare dimensionless onsite interaction strength. In both cases, the adjacent superconductor (SC) and excitonic phases form composite O(5) order parameters, shown in Eq. (62).
respectively be the number of anticommuting and commuting matrix pairs between an order parameter and a four fermion term. Also, let A H and C H respectively be the number of anticommuting and commuting matrix pairs between an order parameter and the single-particle Hamiltonian. We can summarize these relations as 2 Then, various cuts of the global phase diagram in the presence of local interactions is organized according to the following two rules.
(I) Among the available ordered phases, the interaction channel coupled by g M maximally boosts nucleation of the ones satisfying (II) Among the phases selected by (I), an ordered phase is energetically (entropically) favorable if We name the two "or" conditions in (I) selection rules. They are operative at zero as well as at finite temperature 2 We write all matrices as Kronecker products of Pauli matrices as Γµν... = σµ ⊗ σν . . . , therefore any two matrices either commute or anticommute. and chemical potential, and only depend onĥ indirectly, in the sense that the interaction terms by construction preserve the microscopic symmetries of the noninteracting Hamiltonian. Of crucial importance is the fact that MLG and BLG are endowed with the same microscopic symmetries. Therefore (I) cannot distinguish between the two systems. However, (II) can as we show for the extended honeycomb Hubbard model in Sec. II B. Besides (I) having exactly two ways [(a) and (b)] of boosting an ordered phase by a given four-fermion interaction, it leads to the following consequence regarding the adjacent or competing phases. Namely, two broken symmetry phases that reside next to each other in a generic cut of the global phase diagram as we tune the interaction strength, temperature or chemical doping, involving the following two sets of matrices say According to (I) it has to be that O  On the other hand, we name the statement in (II) organizing principle. It arranges the phases, selected by (I), along the temperature axis. Note C H = 0 implies a fully and isotropically gapped state, which causes maximal gain of the condensation energy and is, therefore, favorable at the lowest temperature. On the other hand, C H > 0 results in gapless states at low energies, which, hence is a configuration of higher entropy. Therefore, such states are typically more prominent at higher temperatures. As such (II) is a generalized energy-entropy argument, that goes beyond the binary distinction of whether an order is gapped or gapless, and capable of organizing more than one ordered phases with C H > 0, but varying A H , as the temperature is varied in the system. A schematic representation of the selection rules and the organizing principle is shown in Fig. 1.
Chemical potential plays an important role in promoting superconductivity from pure repulsive interactions by enhancing carrier density. In the above outlined algebraic framework, the chemical doping termμ anticommutes with all pairing orders. Therefore, for a K-component superconducting order parameter A H → A H + K aŝ h →ĥ −μ at finite doping. By contrast, all the excitonic orders commute withμ and thus for an L-component excitonic order parameter A H → A H − L asĥ →ĥ −μ. Therefore, at finite doping pairing states get energetically favored at the lowest temperature, while excitonic orders being more entropic are favored at higher temperature. Still, the competing excitonic and pairing orders for a given interaction channel follow the selection rule (I).
Note that (I) and (II) still leave room for degeneracy among the selected order parameters. For example, when two bilinears gap the quasiparticle spectra and fully anticommute with the matrices in a given interaction channel, then temperature cannot distinguish between them. If, however, one of them is a pairing phase, it will be favored at low temperatures when µ > 0, as it optimally gaps the underlying Fermi surface [selection rule (II)]. We demonstrate the applicability of these rules and display the composite order parameters for the various interaction channels for MLG and BLG in Sec. V. The adjacent orders are tabulated in Tables III and IV. B. Extended honeycomb Hubbard model We exemplify the selection rules and organizing principle outlined in the previous section via one of the simplest but instructive microscopic models for correlated electrons on the honeycomb lattice, the extended Hubbard model, captured by the Hamiltonian Here H 0 describes a collection of noninteracting itinerant electrons on MLG and BLG, where it respectively produces linear and quadratic band touchings, and Here, n σ (R) is the number operator on the site at R with spin projection σ =↑, ↓. The sites located at A (B) span the A (B) sublattice of the bipartite honeycomb lattice, while · · · and · · · denote all pairs of NN and next-nearest-neighbor (NNN) sites respectively. Note in the case of BLG these sites reside on the emergent honeycomb lattice, resulting after integrating out the dimer sites. The terms H U , H V1 , and H V2 then respectively describe the on-site, NN, and NNN interactions with the interaction strengths U , V 1 , and V 2 . We address the effects of repulsive U , V 1 , and V 2 within a RG framework separately in MLG and BLG. At zero temperature and chemical potential, and in the presence of only repulsive on-site interaction (U > 0, but V 1 = V 2 = 0), we find the nucleation of the antiferromagnetic (AFM) ordering in MLG [11,23,30], as well as in BLG [58,80]. See Figs. 4(a) and 4(b), respectively. The induced carrier density upon elevating the chemical potential from the band touching points promotes condensation of electrons into Cooper pairs, hence giving rise to superconductivity. Sufficiently low temperatures but µ > 0 accommodates E g nematic pairing in MLG. However, with µ > 0 the dominant superconductor in the presence of repulsive Hubbard-U in BLG is the singlet Kekulé pairing phase.
Setting U = V 2 = 0, but increasing the strength of repulsive NN interaction (V 1 > 0) at t = µ = 0 results in a charge density wave (CDW) ordering in both MLG and BLG, see Figs. 5(a) and 5(b). Often in the context of BLG, the CDW phase is also called the layer polarized state. Keeping the temperature sufficiently low and extending the Fermi surface by setting µ > 0 gives rise to s-wave and f -wave pairing phases (appearing in a degenerate fashion) in both systems. Therefore, deep inside the paired state these two systems are expected to sustain s + if or f + is superconductor that optimally gaps the underlying Fermi surface [81].
Finally, setting U = V 1 = 0, but increasing the strength of NNN repulsion (V 2 > 0) at t = µ = 0 we find the nucleation of a quantum spin Hall insulator (QSHI) state in both systems, see Figs. 6(a) and 6(b). At the same time, finite chemical doping and sufficiently low temperatures result in an s-wave pairing phase both in MLG and in BLG. Presently, there is an ongoing debate regarding the fate of topological insulators in MLG resulting from NNN repulsion at the half-filling [22,25,32,48,82]. All the exact numerical diagonalizations have been performed for spinless fermions [25,32,48], and the realization or absence of quantum anomalous Hall insulator still remains unsettled. Here, from an unbiased RG analysis we show that at least for spin-1/2 fermions NNN repulsion supports QSHI in MLG. Similar RG analysis for spinless fermions clearly suggest the presence of quantum anomalous Hall insulator in half-filled MLG for strong enough NNN repulsion, which will be discussed in a forthcoming publication.
The resulting broken symmetry phases in the phase diagram of the extended Hubbard model, as well as the similarities or distinctions among the orders realized in MLG and BLG, are in accordance with our proposed selection rules. Namely, the s-wave and f -wave pairing phases, adjacent to the CDW phase in both systems with a finite NN repulsion, fully anticommute with the CDW order parameter, and form with it O(3) and O(4) composite vectors, respectively, while also opening an isotropic spectral gap in both MLG and BLG. Similarly, the s-wave pairing and QSHI orders, appearing adjacent to each other in the presence of NNN repulsion, constitute an O(5) composite order parameter. Like the s-wave pairing, QSHI ordering also gaps the spectrum in both MLG and BLG. On the other hand, the AFM order parameter, which is the dominant excitonic phase in both systems with finite on-site repulsion, fully anticommutes with the E g nematic superconductor (SC) and the singlet Kekulé SC, and forms O(5) composite vectors with them. The singlet Kekulé pairing fully commutes with the Dirac Hamiltonian in MLG. Therefore in the presence of finite chemical doping a partially anticommuting (with the Dirac Hamiltonian) E g nematic pairing is favored in this system. In contrast, the singlet Kekulé SC represents a fully gapped phase in BLG [8], and is thus favored over the partially anticommuting E g nematic SC. The phase diagrams of the extended Hubbard model share some common features. For example, with increasing chemical doping (1) any ordering sets in at stronger coupling, (2) the requisite strength for the repulsive interaction for the onset of any excitonic order gets pushed toward stronger coupling, and (3) range of interaction over which pairing phase is realized increases.
More detailed and quantitative analysis of the extended honeycomb Hubbard model is presented in Sec. VI. On a technical note, even though the terms H U , H V1 , and H V2 translate to linear combinations of linearly independent interaction terms in the continuum formalism, the effect of these quartic terms can be captured qualitatively by increasing the strength of the local interactions in the AFM, CDW, and QSHI channels, respectively, see Sec. V. This is because these coupling constants are the ones that diverge dominantly under coarse grain while increasing the bare strengths of U , V 1 , and V 2 of the extended Hubbard model.
Finally, we note that our conclusions for the Hubbardlike model is consistent with the findings from nonperturbative numerical analysis on honeycomb lattice [83,84]. Specifically, a functional RG analysis with only NN Coulomb repulsion among spinless fermions (thus accommodating only spin-triplet pairings) reported charge density wave at and near half-filling, and f -wave pairing away from it [83], see Fig. 5. A quantum Monte Carlo simulation found nucleation of s-wave pairing by doping honeycomb monolayer, residing in the vicinity of a quantum spin Hall insulator [84], see Fig. 6. Thus the proposed selection rule and organizing principle, resting on (anti)commutation relations among matrices in the interaction channel, order parameter and noninteracting Hamiltonian, should be decisive in determining the nature of ordering in correlated multiband systems.

C. Organization
The rest of this paper is organized as follows. In Sec. III we review the lattice descriptions and the corresponding tightly bound electron models of MLG and BLG. In Sec. IV A we establish a low-energy continuum theory of the above models by means of Fourier expansion around the band touching points. In doing so, we carry over the microscopic symmetries of the honeycomb lattice so that the two descriptions are symmetry equivalent. We then leave the realm of single-particle physics and introduce short range electron-electron interactions in Sec. IV C, which we address in a perturbative fashion via Wilsonian momentum-shell RG and the expansion schemes. A review of the possible broken symmetry phases is then provided in Sec. IV D. Section V is dedicated to addressing the various interaction channels one-by-one and demonstrating the realization of ordered states and their compatibility with the selection rules. Finally we revisit the problem of the extended honeycomb Hubbard model in greater details in Sec. VI. A summary of our findings is presented in Sec. VII. Additional technical details are relegated to the appendices.

III. LATTICE MODELS
In this section we introduce single-particle tight binding descriptions of electrons on the honeycomb monolayer and bilayer. Our model includes only NN hopping for MLG. But, in order to capture the coupling between two layers, we incorporate interlayer hopping terms in BLG, besides the intralayer one. In what follows we restrict ourselves to carbon-based honeycomb lattices, and as such will completely neglect spin-orbit coupling, which is extremely small for carbon atoms [4].

A. Monolayer graphene
The simplest tight-binding Hamiltonian for an isolated sheet of graphene describes a bipartite lattice consisting of two interpenetrating triangular sublattices A and B with the NN hopping amplitude (t) as [85,86] where a(A) are fermion annihilation operators on the A sublattice, generated by linear combinations of the primitive lattice vectors v 1 = √ 3a(1, 0) and v 2 = a/2( √ 3, 3), with a being the distance between neighboring atoms. On the other hand, b(B) are fermion annihilation operators on the sites belonging to the B sublattice that are obtained as B = A + δ i , with δ 1 = a/2(− √ 3, −1), δ 2 = a/2( √ 3, −1) and δ 3 = a(0, 1), see Fig. 2(a). The reciprocal vectors are given by q 1 = 2π/a(1/ √ 3, −1/3) and q 2 = 4π/a(0, 1/3). The BZ has the same hexagonal structure as the real space lattice, only rotated by 90 degrees. For convenience we shift the origin of the BZ to the middle of the hexagon (the Γ point), see Fig. 2(b). The lack of the inversion symmetry of the honeycomb lattice about any given site (see Fig. 2(a)) results in Fermi points in the band structure of H MLG instead of an extended Fermi surface, as obtained on a square lattice, for example. The spectra of H MLG read where and f (k) vanishes linearly at the six corners of the BZ [85]. Thus, MLG k supports six Dirac points, around which the low energy dispersion is linear and isotropic [87]. However, only two of these band touching points are inequivalent, which we take to be at ±K, with K = ± 4π 3 √ 3a (1, 0). The rest are connected to these two by primitive reciprocal vectors [86].

B. Bilayer graphene
One generates the Bernal stacked BLG by adding another graphene sheet on top of Fig. 2(a) and shifting it by e.g. −δ 3 , see Fig. 2(c). Due to the two layers the minimal model has twice as many degrees of freedom when compared to MLG. We denote the layer index (i) with a subscript on the annihilation (a i and b i ) and creation (a † i and b † i ) operators, where i = 1, 2. The tight-binding Hamiltonian including intralayer and interlayer hopping terms is of the form [8,80] where Here H describes the intralayer NN hopping, where B ± i = A ± δ i , while the remaining terms encompass the interlayer ones. The terms H ⊥ AA and H ⊥ BA represent remote interlayer hopping processes, which we will drop from now on. The direct hopping between the dimer sites (a 1 and b 2 ) is captured by H ⊥ AB . After setting t AA = t BA = 0, the spectra of H BLG describe four bands with the dispersions Since f (k) vanishes linearly at ±K, BLG 1,k describes two quadratic band touchings, whereas the two BLG 2,k bands are gapped everywhere, also called the split-off bands. Note that H ⊥ BA produces trigonal warping that splits each quadratic band touching into four Dirac points [9]. In Sec. VII, we qualitatively discuss the role of trigonal warping on the phase diagrams of interacting electrons in BLG.
In this paper we study these two systems in terms of their low-energy effective theories. As the continuum models are well established within the scientific literature [4], in presenting the respective single particle descriptions we will only quote the cornerstones of the derivations.

IV. LOW-ENERGY EFFECTIVE THEORY
We formulate the effective continuum theories by retaining the low energy degrees of freedom of both lattice models, Eqs. (8) and (12). As the split-off bands in BLG are at higher energies, and the eigenfunctions of BLG 2,k are dominantly localized on the overlapping a 1 and b 2 sites, we can project out these bands with relative ease [4,8,80]. After integrating out the a 1 and b 2 sites, the resulting lattice assumes another honeycomb structure, see Fig. 2(d). As such, the effective BZ has the same structure as that of MLG, with the important distinction that the band touching at the six corners are now quadratic, in comparison to the linear band touching in MLG. Moreover, the microscopic symmetries of these two models (reflections, translation, and time reversal) are identical, about which more in a moment. Consequently, the part of the action describing local electron-electron interactions takes the same form in these two systems, see Sec. IV C. This allows us to investigate the role of the normal state band structure in emergent phases within the same symmetry group.
Due to the small atomic number of carbon atoms we neglect the spin-orbit coupling in both MLG and BLG, and introduce the spin degree of freedom as a mere doubling of the Hamiltonian. Besides spin, we will also introduce a particle-hole index using the Nambu doubling, which facilitates the inclusion of both excitonic and superconducting phases within a unified framework. Note that the split-off bands now being projected out, the two effective models have the equal number of low energy degrees of freedom.

A. Non-interacting models
Let us add a spin index to the fermionic degrees of freedom and expand the fields in terms of their Fourier components as where τ and ω are respectively imaginary time and frequency, r = a and b for electrons on the A and B sublattices, respectively, with the spin orientation s =↑, ↓.
We retain the modes near the two inequivalent valleys at ±K, and write r s (ω, ±K + k) = r ±K s (ω, k). The eightcomponent spinor in Fourier space is then structured according to where r v s are annihilation operators for fermions on the A and B sublattices with valley index v = ±K, and spin projection s =↑, ↓, and denotes transposition. Finally, we introduce the particle-hole degree of freedom as the "outermost" index and write a sixteen-component Nambu doubled spinors as where in the lower block we absorbed the unitary part of the time reversal operator U = Γ 210 [see Sec. IV B], so that Ψ ω,k and Ψ Nam transform identically under all symmetries. For the sake of brevity we suppress the subscript "Nam" from now on, and unless otherwise mentioned, thus Ψ ≡ Ψ Nam and Ψ † ≡ Ψ † Nam . To capture the symmetry properties of MLG and BLG, and subsequently develop the interacting models describing these two systems, we write the matrices acting on the sixteen-dimensional spinors as Γ µνλρ = η µ ⊗ σ ν ⊗ τ λ ⊗ α ρ where {α}, {τ }, {σ} and {η} are four sets of Pauli matrices that operate on the sublattice or layer, valley, spin and particle-hole degrees of freedom, respectively. Here, µ, ν, λ, ρ = 0, 1, 2, 3, and the index 0 always corresponds to the unit matrix. In Eq. (17), Γ 210 = σ 2 ⊗τ 1 ⊗α 0 . Note that having integrated out the overlapping a 1 and b 2 sites in BLG, the sublattice and layer description is now redundant, as all the remaining a and b sites reside on layers 2 and 1, respectively. Therefore, from now on we refer to this degree of freedom simply as sublattice. Also note that by virtue of writing all sixteen-dimensional matrices in terms of Kronecker products of Pauli matrices, any two such matrices either commute or anticommute.
The non-interacting effective low-energy theory of MLG is composed of linearly dispersing chiral relativistic Dirac fermions, that persist all the way up to a nonuniversal ultraviolet momentum cutoff Λ ∼ 1/a. In comparison, the corresponding effective theory of BLG is constituted by quadratically dispersing chiral quasiparticles (also up to a ultraviolet momentum cutoff Λ). We write the Hamiltonians for MLG and BLG respectively aŝ where v = ta is the Fermi velocity and m = t ⊥ /(2t 2 a 2 ) is the effective mass. The momentum-dependence is described by the p-wave and d-wave cubic harmonics, that are related to the spherical harmonics Y m l with l = 1 and l = 2, respectively, such that We also introduced a chemical potential (µ) term, which facilitates tuning the Fermi energy to and away from the band touching points. Therefore, µ plays the analogous role as the gate voltage. In this paper we only consider electron like doping (µ > 0). All the conclusions are equally germane for hole-doped (µ < 0) systems. The chemical doping develops an extended (one-dimensional) Fermi surface from the (zero-dimensional) Fermi points, and increases the carrier density, which is conducive for the condensation of electrons into Cooper pairs.
The above models then constitutes the following noninteracting Euclidean (imaginary time) action with j =D and L respectively for MLG and BLG. Here d is the spatial dimensionality, Ψ † (τ, r) and Ψ(τ, r) are the inverse Fourier transform of the spinors Ψ † and Ψ, and τ is the imaginary time.

B. Symmetries and action
The microscopic symmetries of the honeycomb monolayer and bilayer are identical and correspond to the D 3d point group. Here we cast these symmetries in terms of reflections and rotations, and show how they manifest in the low-energy theory. We augment the point group symmetries with translation invariance, discrete time reversal symmetry, as well as continuous spin SU(2) and pseudospin SU p (2) rotational invariance.
The two reflection symmetries of the honeycomb lattice are sublattice (S) and valley (T ) reflections [10]. They respectively take the form and Γ3013, Γ3020 smectic charge density wave g s Γ0s13, Γ0s20  Reflection S exchanges the sublattices A↔B, but leaves the valleys invariant. Whereas under T , K ↔ −K, and hence it exchanges the Dirac points, but does not affect the sublattice degree of freedom. Both the linear and quadratic band touching points are rotationally invariant at low energies. Let us denote the orbitals as v i = p i or d i for i = 1, 2. Then, rotations by an angle θ in orbital space take the form R(θ) = exp(iθΓ 0033 /2), and a rotation of the Hamiltonian corresponds to an ordinary vector rotation of the two-component vector For example, under rotation by θ = π/2 The corresponding rotation of momentum axes by an angle θ = −π/(2l) takes the form where the angular momentum l = 1 and 2 for p and d orbitals, respectively. The requisite rotation by the angle θ = −π/(2l) for the invariance of the Hamiltonian indicates that respectively the linear and biquadratic band touching represent momentum space vortices of vorticity l, and the wave functions of the valence and conduction bands are eigenstates of orbital angular momentum l. Translation by a vector R acts on the electron fields as Therefore the matrix transforming the 16-component Nambu spinors under translation has to pick up opposite signs at the valleys at ±K. A translation of the spinor Ψ is then written as with Γ 0030 being the generator of translations. More intuitively, the above transformation in real space reads The reversal of time is implemented by the antiunitary operator T = U K. Here U = Γ 0210 is a unitary operator and K is complex conjugation. For spin-1/2 electrons, we indeed find T 2 = −1.
The Hamiltonian operators in Eqs. (18) and (19) are manifestly invariant under SU(2) spin rotations, generated by Γ 0s00 , with s = {1, 2, 3}. Finally, the pseudospin SU p (2) transformations are generated by Here Γ 1003 and Γ 2003 rotate between an excitonic (such as, charge density wave) and two (real and imaginary) components of a paired state (such as, s-wave pairing). While the number operator Γ 3000 rotates the U(1) superconducting phase [88]. See Tables I and II. Having established the non-interacting continuum model and the operative symmetries, in the next section we build on this foundation and examine the effects of electron-electron interactions.

C. Electron-electron interactions
We incorporate short range or local electron-electron interactions by adding all symmetry-allowed momentumindependent four fermion terms to the action S j 0 . The schematic form of the interacting part of the action in the presence of such a local quartic term is where Ψ ≡ Ψ τ,r and Ψ † ≡ Ψ † τ,r , M and N are 16 dimensional (in the Nambu doubled basis) Hermitian matrices, and g M N is the corresponding coupling constant.
Note from Eq. (21) the scaling dimension of the spinors is [Ψ † τ,r ] = [Ψ τ,r ] = d/2, such that S j 0 is scale invariant for j =D and L. Consequently, a coupling constant of local quartic terms scales as [g M N ] = z − d, where z is the dynamical scaling exponent, determining the relative scaling between energy and momentum according to k ∼ |k| z . Therefore z = 1 (2) in case of Dirac (Luttinger) fermions. Consequently, S int remains scale invariant. As d = 2, sufficiently weak short range electron-electron interactions among two-dimensional Dirac fermions are irrelevant, and the Dirac cones remain stable up to a critical strength of interactions, where, however, the system undergoes a quantum phase transition into an ordered phase. On the other hand, in a two-dimensional system of Luttinger fermions local quartic interactions are marginal, and the quadratic band touching point is unstable in the presence of infinitesimal interactions. In terms of the density of states (DOS), the Dirac system exhibits linearly vanishing DOS, namely ρ(E) ∼ |E|, while in the case of Luttinger fermions the low-energy DOS is constant i.e., ρ(E) ∼ constant. The scaling of the DOS determines the (ir)relevance of local interactions.
Throughout we neglect the long range tail of the Coulomb interaction. In an undoped Dirac liquid, the long range Coulomb interaction is a marginally irrelevant coupling (due to the vanishing DOS), irrespective of its bare strength, that otherwise causes only a logarithmic enhancement of the Fermi velocity [11,[89][90][91][92][93][94][95][96][97]. On the other hand, in undoped Luttinger liquid long-range Coulomb interaction gets screened due to the finite DOS, resulting in local density-density interaction, thus only renormalizing the bare strength of the coupling g s 1 [12] [see Eq. (32)]. At finite temperature this approximation is further justified due to the thermal screening [98], while at finite doping conventional Thomas-Fermi screening sets in. Nevertheless, long-range Coulomb interaction can shift the phase boundaries between the ordered and disordered phases, without altering the nature of the former one, as has been shown for both two-and threedimensional Dirac liquids [99,100].
All the local quartic terms appearing in the interacting Lagrangian must transform as scalars under all operative symmetries introduced in Sec. IV B, which severely restricts the number of couplings constants to 18. The interacting Lagrangian (L int ) containing all symmetryallowed local quartic terms reads [10] where   Tables III and IV. See also Sec. V A.
Visibly, L int separates into spin singlet (L sing int ) and spin triplet (L trip int ) parts, containing the matrices σ 0 and σ s (with s = 1, 2, 3), acting on the spin index, respectively. The sublattice and valley degrees of freedom (accompanied by α and τ matrices, respectively) repeat for the two irreducible representations of the SU(2) spin algebra, while the particle-hole index (η matrix) is exactly the opposite between singlet and triplet channels. Since there is only one operative SU(2) algebra in the spin sector, this suggests a redundancy in the above expression of the interacting Lagrangian. This expectation can be verified using the Fierz identity for eight-dimensional Hermitian matrices (note the particle-hole degree of freedom in this sense is an artificial doubling), which reveals that there are altogether only nine linearly independent four fermion terms. Without the loss of generality, we choose these to be the spin singlet interactions (g s i ). For a detailed presentation of the Fierz reduction see Appendix B. Throughout, here we focus on repulsive electron-electron interactions, which in our notation corresponds to g j µ > 0 for µ = 1, . . . , 9 and j = s, t. To shed light on the structure of the global phase diagram of interacting two-dimensional Dirac and Luttinger fermions, we perform Wilsonian momentum shell RG analysis at zero, as well as finite temperature (T ) and chemical potential (µ). We already established the scaling of the coupling constants to be g j µ = z − d, which pins the lower critical dimension of the corresponding theories at d = z and facilitates an expansion around the marginal dimensionality with = d − z [101][102][103]. Notice that the physical values of are 1 and 0 for twodimensional Dirac and Luttinger fermions, respectively. However, our conclusions are insensitive to the exact value of , except the nonuniversal locations of the phase boundaries between the ordered and disordered phases.
To proceed further, we introduce an ultraviolet momentum cutoff Λ, that replaces the (hexagonal) lattice Brillouin zone with a spherical one around two inequivalent valleys in the continuum theory. Moreover, we introduce the dimensionless temperature and chemical potential. These two quantities respectively for Dirac (D) and Luttinger (L) systems are defined as From now on we omit the tilde and denote the dimensionless chemical potential by µ. Furthermore, for brevity we suppress the D and L indices. The scaling dimensions of these two quantities are [t] = [µ] = z, and hence they are relevant parameters. The RG flow equation or the β function of a quantity describes its behavior under coarse grain as we integrate out the fast Fourier modes within a Wilsonian momentum shell Λe − < |k| < Λ, where is the logarithm of the RG scale. The flow equations for dimensionless temperature, chemical potential and coupling constants can be summarized respectively as The coupling constants appearing in the RG flow equations are also dimensionless, obtained by taking For the exact form of the functions H i jk (t, µ) see Appendix C. Eq. (35) implies that all orderings take place at critical couplings g s i ∼ , which is consistent with our observation of local interactions being irrelevant and marginal in MLG and BLG, respectively. Therefore, in a Dirac system all orderings set in beyond a finite strength of the interactions, while a Luttinger system is conducive for ordering even for sufficiently weak interactions.
Our methodology of detecting phase transitions is as follows. As temperature and chemical potential "grow" under the RG transformation, they provide an infrared cutoff when either one of them becomes of order of unity. The renormalized or scale dependent temperature and chemical doping take the closed analytic forms t( ) = t(0) exp[z ] and µ( ) = µ(0) exp[z ], respectively. Thus, with finite bare values (at = 0) of t and µ, denoted by t(0) and µ(0), respectively, we must terminate the RG flows of the coupling constants at a scale For a given interaction channel (g i ), we integrate the differential equations in Eq. (35) between 0 ≤ ≤ * for increasing g i (0). Its smallest bare value for which at least one coupling constant diverges (not necessarily g i ) indicates a phase transition to an ordered phase and contributes a data point (t, µ, g * i ) to the phase boundary, where g * i is the critical interaction strength in the given interaction channel. However, to reveal the type of instability and the nature of the ordered phase we have to add order parameter fields to the Lagrangian and investigate their fate under the coarse grain, which we do next. It is worth pointing out that here we determine the infrared cutoff ( * ) for the RG flow of coupling constants by setting the renormalized temperature t( ) = 1 and chemical potential µ( ) = 1. However, the nature of the ordered states is insensitive to such choice, as long as t( ) ∼ 1 and µ( ) ∼ 1. Only the nonuniversal details, such as the transition temperature (t c ), interaction coupling at which ordering switches from pairing to excitonic, depend on it.
Γα013, Γα020  (2) rotations, respectively. Here + (−) means even (odd), whereas 0 (1) means the bilinears transform as scalar (vector) under the corresponding rotation. In the TR column a symbol outside (within) parentheses corresponds to α = 1 (2). The order parameters marked by the same roman numeral in the SUp (2) column are related to each other by pseudospin SUp(2) rotations, which always relates a pairing order to an excitonic one, see Table I. The last two columns display whether the order parameter fully gaps fermions in MLG and BLG, where () means gapped (gapless) spectrum.

D. Broken symmetry phases
To identify the nature of the broken symmetry phase in an unbiased fashion we introduce all symmetry-allowed local fermion bilinears to the action, which assume the schematic form ∆ M (Ψ † M Ψ), where ∆ M is the corresponding conjugate field and M is a sixteen-dimensional Hermitian matrix. The ordered state is characterized by the expectation value of the fermion bilinear Ψ † M Ψ = 0. Note while the RG analysis of the coupling constants g i can be performed using eight-dimensional matrices (without invoking the Nambu doubling), one has to utilize the Nambu basis in the renormalization of order parameter fields to account for the superconducting channels. The effective single-particle Hamiltonian containing all symmetry allowed local orders reads where h exc = h sing exc + h trip exc , with and whereas Here α = 1, 2, ∆ p i,1 = ∆ p i cos φ, ∆ p i,2 = ∆ p i sin φ, and φ is the U(1) superconducting phase. According to H local , 27 order parameters are organized into three groups, each containing nine entries. Namely h sing exc and h trip exc contain nine spin singlet and nine spin triplet excitonic bilinears, respectively, while h pair accommodates five singlet and four triplet pairing orders. The properties of fermion bilinears e.g. their irreducible point group representations, physical meanings, and transformations under various discrete and continuous symmetries are displayed in Tables I and II, where we also collect the notation that is used throughout this paper.
Besides the 9 independent coupling constants g s i we also renormalize the 27 order parameter fields from H local and arrive at their β-functions taking the schematic form where j =D and L respectively for Dirac and Luttinger fermions in MLG and BLG. Here we absorb the contribution from the scaling dimension (z) of the conjugate fields ∆ into its β-function. For the exact form ofβ j Notice that along some phase boundaries multiple order parameter fields diverge in a degenerate fashion. However, as we cut the multi-dimensional space of couplings in a very specific way (e.g. along the axes g µ i ), this does not necessarily indicate an enlargement of symmetry among these orders at the governing quantum critical point, which will be discussed in a forthcoming paper. Rather, this suggests that the respective order parameters simultaneously diverge with a specific phase locking of the internal degrees of freedom. For example, when the divergences of the conjugate fields for both s-wave and f -wave pairings are degenerate, the ordered state is expected to support s + if or f + is pairing [81], which produces a maximal gap in the quasiparticle spectra.
To facilitate further discussion on order selection in two-dimensional Dirac and Luttinger fermions, let us briefly outline some properties of the available broken symmetry phases. In this section we quote the symbol of the respective orderings, with which they are referred to in the phase diagrams (Figs. 7-10). The symbols, together with nomenclature and various symmetry properties of fermion bilinears are collected in Tables I and  II.

Particle-hole or excitonic orders
Due to the neglected spin-orbit coupling the excitonic order parameters (just like the interaction terms) fall into two distinct categories. Namely, for each spin singlet order, where all matrices come with the σ 0 matrix, there exists a spin triplet analog appearing with the σ s matrices, where s = {1, 2, 3}, operating on the spin index. We indicate this distinction by s (singlet) and t (triplet) superscripts on the corresponding conjugate fields. Since the Hamiltonians in Eqs. (18) and (19) are oblivious to spin or invariant under the rotation of the spin quantization axis, this distinction does not affect the emergent fermionic quasiparticle spectra inside the ordered state. Note, however, that any triplet ordering breaks SU(2) spin rotational symmetry and is hence accompanied by two massless Goldstone modes.
The fermionic density and spin density (respectively denoted by A s 1g and A t 1g ) do not break any discrete lattice symmetries, while the chiral and spin-chiral chemical potential ( denoted by A s 1u and A t 1u , respectively) break valley reflection symmetry (T ). All four orders commute with both the Dirac and Luttinger Hamiltonians, resulting in a trivial renormalization of their conjugate fields at zero temperature, see Appendix A. Consequently, none of them is realized in the ordered phase at t = 0.
Both the Dirac and Luttinger systems altogether accommodate six mass orders in the particle-hole channel [8,105]. These phases introduce an isotropic gap in the quasiparticle spectrum and are thus favorable at zero and low temperature. All of them fully anticommute with the noninteracting Hamiltonian. The quantum Hall states break both sublattice (S) and valley (T ) reflection symmetries. The quantum anomalous Hall insulator (QAHI), denoted by A s 2g , additionally breaks time rever-  sal symmetry, while it is restored in its spin triplet counterpart, the quantum spin Hall insulator (QSHI), denoted by A t 2g [50,54,106]. The charge density wave (CDW) and antiferromagnet (AFM) order parameters break sublattice reflection, but preserve valley reflection symmetry [10,11]. They are respectively denoted by A s 2u and A t 2u . The QAHI, QSHI, CDW and AFM orderings represent massive phases in both Dirac and Luttinger liquids.
Both models accommodate two additional masses in the family of Kekulé orders.
The Kekulé orderings break translational symmetry of the honeycomb lattice and result in the enlargement of the unit cell, while preserving the rotational invariance [55]. In MLG, the Kekulé valence bond solid (KVBS) and spin-Kekulé solid (sKS) are fully gapped phases, while the Kekulé current (KC) and spin-Kekulé current (sKC) order parameters fully commute withĥ D . Exactly the opposite is true in BLG, namely KC and sKC represent spin singlet and spin triplet masses for Luttinger fermions, respectively, whileĥ L fully commutes with the KVBS and sKS order parameters [8]. The KVBS, sKS, KC and sKC phases are respectively denoted by A s 1k , A t 1k , A s 2k , and A t 2k . All rotational symmetry breaking phases are gapless, and due to the higher density of states at low energies they are more entropically favorable than their gapped counterparts. Such orderings in MLG merely shift the location of the Dirac points, while in BLG it splits each quadratic band touching into two linear Dirac cones [8]. Consequently, all the nematic and spin nematic orders cause power-law suppression of the DOS for Luttinger fermions, yielding ρ(E) ∼ |E| at low energies. Four nematic phases, transforming under the E g and E u representations of the D 3d point group, preserve translational invariance. Two of them are spin singlets (denoted by E s g and E s u ), while two spin-nematic order parameters (denoted by E t g and E t u ) additionally breaks the spin-rotational invariance. In addition to rotational invariance, translational symmetry is also broken in the two smectic phases, the smectic charge density wave and smectic spin density wave orderings, respectively denoted by E s k and E t k . The smectic phases also produce gapless, but anisotropic Dirac points in the ordered phases.

Particle-particle or superconducting orders
When we add chemical doping to these nodal systems, instead of Fermi points, they develop extended Fermi surfaces around the two independent band touching points. The chemical potential term is of the form µΓ 3000 , which commutes with all excitonic orders, therefore they cannot gap such metallic systems. On the other hand, all pairing orders include η α = {η 1 , η 2 } in the particle-hole sector, and hence they can conceivably anticommute with h j (k) when µ is finite. Therefore, low temperatures (when a gapped spectrum is favorable due to the gain of condensation energy) and finite µ are conducive for superconducting orders, as the increased carrier density is conducive for the condensation of electrons into Cooper pairs. The Nambu basis altogether accommodates nine local or momentum-independent pairing orders, which we mark with a superscript p in Table II. The singlet swave (denoted by A p 1g ) and triplet f -wave (denoted by A p 1u ) [41] superconductors (SCs) represent fully gapped phases in MLG and BLG. The singlet and triplet Kekulé SCs break translational symmetry. The singlet Kekulé SC (denoted by A p 1k ) represents a mass order in the Luttinger system (BLG), but fully commutes with the Dirac Hamiltonian [8]. In the exact opposite way, the triplet Kekulé SC (denoted by A p 2k ) gaps a Dirac liquid, but fully commutes withĥ L [16]. Two nematic superconductors, transforming under the E g and E u representations of the D 3d group (denoted by E p g and E p u , respectively), break rotational symmetry and produce gapless quasiparticle spectra. Furthermore, the smectic pairing (denoted by E p k ) breaks rotational and translational invariance, and also leads to gapless bands. While the nematic and smectic pairings produce point nodes in the spectra of emergent Bogoliubov-de Gennes quasiparticles, the remaining two pairings transforming under the A 2g and A 2u representations (denoted by A p 2g and A p 2u , respectively) produces Fermi surfaces inside the ordered phases.

V. PHASE DIAGRAMS
In this section, we scrutinize multiple cuts of the global phase diagram of interacting fermions in MLG and BLG by systematically increasing the strength of interactions in one channel at a time, both at zero and finite chemical potential. The fact that the form of contact interactions in MLG and BLG is identical and they only differ in their noninteracting Hamiltonians gives us an opportunity to further dissect the selection rules outlined in Sec. II A and assess the imprint of normal state band structure on the global phase diagram. This is because selection rule (I) only depends on S int while (II) only on S 0 . Consequently, (I) cannot distinguish between MLG and BLG, and therefore any deviation must be rooted in (II).
Our methodology, involving Wilsonian momentum shell RG analysis and expansion around the lower critical dimension, is as outlined in Secs. IV C and IV D. Each phase diagram is displayed in the (g i µ , t) plane, where g i µ is a dimensionless (bare) coupling constant with µ = 1, . . . , 9 and i = singlet (s) and triplet (t), and t is dimensionless temperature. We scan for the phase boundaries by keeping temperature constant and increasing the bare interaction strength, and detect the phase boundary and the nature of symmetry breaking between disordered and ordered phases, marked in the phase diagrams respectively as white and gray regions. Respectively, at zero and finite doping the disordered phase represents chiral nodal Fermi liquids (with point nodes) and regular Fermi liquids (with extended Fermi surface). Phase boundaries from the disordered state into an excitonic (superconducting) order are indicated with blue (red) lines. The RG procedure is only equipped to detect the divergence of some conjugate field that directly couples with a fermion bilinear, indicating onset of the corresponding ordered state and hence a phase transition. But, it does not tell us about potential regions of coexistence of adjacent phases deep inside the ordered phase. Also note that the critical interaction strength and the transition temperature are non-universal quantities.
Four cuts of the phase diagram for each interaction channel (for MLG and BLG, at zero and finite chemical potential) are shown in Figs. 7-10. All phase transitions and the nature of the supervector order parameters of adjacent phases are summarized in Tables III and IV. Notice that we always cut the global phase diagram along CC Monolayer graphene (MLG) Bilayer graphene (BLG) Low t High t Symmetry Low t High t Symmetry Low t High t Symmetry Low t High t Symmetry 2⊗O (5) 1⊗O (3) 1⊗O (5) 1⊗O (4) A one certain axis in the space of the coupling constants, which is generally not how a physical system behaves. However, unless fine tuned, among multiple running couplings there will always be one that diverges the fastest. The selection rules are then determined by the strongest diverging coupling constant. However, effective quartic interactions in a certain channel can in principle be tuned in Determinantal quantum Monte Carlo (DQMC) simulations [107], obtained by integrating out bosonic order parameter fluctuations. For brevity, when quoting a set of matrices we will use α = {1, 2} in the particlehole or Nambu sector of the Hilbert space, indicating off-diagonal or pairing orders, and s = {1, 2, 3} in the spin sector, indicating a spin triplet order parameter.

A. Quartic interactions: mass channels
For quartic interactions in the four excitonic mass channels that are common across the two systems, selec-tion rule (Ia) is the decisive at µ = 0, namely O i = M i , where O i (M i ) are the matrices involved in the order parameter of the broken symmetry phase (four fermion term of the given interaction channel). At µ > 0 and low temperatures, we observe nucleation of "adjacent" superconducting phase(s). Next we discuss these cases.
(1) For quartic interaction in the QAHI channel M = Γ 0033 , which fully commutes with all superconducting masses. These orders are therefore not available to condense into. Instead, in both systems we observe nucleation of the E u and E g nematic SCs (denoted by nematic SC 1 and nematic SC 2 , respectively, in Tables I and II) and the smectic SC phases. These pairing order parameters indeed fully anticommute with M [selection rule (Ib)]. The high temperature phase is the QAHI, in accordance with selection rule (Ia). The multicomponent order parameters combining QAHI and various nematic and smectic SCs in this case constitute the following su- The corresponding cuts of the phase diagram are displayed in Fig. 7a.
(2) For quartic interaction in the CDW channel M = Γ 3003 , which fully anticommutes with two superconducting masses [selection rule (Ib)], the singlet s-wave and triplet f -wave pairings, which form following the composite order parameters with CDW: Indeed we observe degenerate nucleation of these two pairing phases at low temperatures in MLG and BLG, when µ > 0. The adjacent excitonic ordering is the CDW [selection rule (Ia)], which sets in at higher t where superconductivity is destroyed, see Fig. 7b. The absence of order differentiation between Dirac and Luttinger systems with quartic interaction in the CDW channel is due to the fact that all dominant order parameters (one excitonic and two superconducting) fully anticommute witĥ h D as well as withĥ L . At zero doping, we only find CDW ordering in both MLG and BLG [selection rule (Ia)]. In Sec. VI we demonstrate the same result for NN repulsion.
(3) For quartic interaction in the QSHI channel M = Γ 3s33 , which fully anticommutes with two superconducting order parameters. One of them, the s-wave pairing, is a mass order in both systems, while the singlet Kekulé SC gaps only Luttinger fermions. Accordingly, in the Dirac system at low temperatures and finite chemical potential only s-wave SC can be observed. On the other hand, we find degenerate nucleation of s-wave and singlet Kekulé pairings in BLG. The nucleating orders and their differentiation through the underlying band structure follows from the selection rules (Ib) and (II), respectively. The high temperature phase in either case is the adjacent QSHI order, which follows from selection rule (Ia). The relevant cuts of the phase diagram are shown in Fig. 7c. The supervector order parameters formed by adjacent phases read At zero doping we always find QSHI in MLG and BLG, according to selection rule (Ia).
(4) For quartic interaction in the AFM channel M = Γ 0s03 , which fully anticommutes with the E g nematic and singlet Kekulé SCs. Latter one represents a mass order for Luttinger fermions, but fully commutes with the Dirac Hamiltonian. Therefore Kekulé SC is not a viable candidate for Dirac fermions at low temperatures. Rather, the maximally anticommuting pairing for the linearly dispersing Dirac electrons is the E g nematic SC.
The composite order parameters in these two cases read Therefore, at finite µ the low temperature paring state in MLG (BLG) is singlet E g nematic (Kekulé) SC, which follows from selection rule (Ib). The high temperature phase at finite µ in both systems is the AFM [selection rule (Ia)], see Fig. 7d. At zero doping we always find AFM in both MLG and BLG. In Sec. VI we again support this observation in the phase diagram of the honeycomb Hubbard model with only on-site repulsion.

B. Quartic interactions: nematic channels
Nematic order parameters break rotational symmetry and they do not introduce a gap in the Luttinger or Dirac spectrum. Therefore other gapped phases might be more favorable at low temperatures when we tune the quartic interactions in the nematic channels, especially if they satisfy the selection rule (Ib). Also note that in general spin fluctuations may induce charge ordering, while the opposite process is comparatively suppressed. This is because in the spin sector of the Hilbert space σ i σ j (with i, j = 1, 2, 3) can conceivably equal to σ 0 when i = j, but σ 0 σ 0 = σ i . Here two Pauli matrices are appearing from the interaction vertices, and their product results from the corresponding one-loop Feynman diagrams (see Fig. 13). As a result, the spin singlet nematic interactions at zero chemical doping or high temperatures follow selection rule (Ia). In contrast, the appearing excitonic phases for interactions in their spin triplet counterparts follow selection rule (Ib). Next we systematically discuss these outcomes.
(1) For quartic interaction in the E g nematic (nematic 1 ) channel M = (Γ 3001 , Γ 3032 ), which fully anticommutes with the s-wave and f -wave pairings, and forms the following supervector order parameters: These two pairing phases are fully gapped in both MLG and BLG, and indeed we observe degenerate appearance of them with increased chemical doping in both systems [selection rule (Ib)]. At µ = 0 or at high temperatures we find nucleation of the E g nematic order [selection rule (Ia)] in both systems, see Fig. 8a.
(2) For quartic interaction in the E u nematic (nematic 2 ) channel M = (Γ 0031 , Γ 0002 ). These interaction matrices anticommute with the singlet Kekulé SC, which gaps the Luttinger fermions, but not the Dirac fermions. Therefore, in BLG at low temperatures we observe the singlet Kekulé SC together with the degenerate QSHI and AFM orders [selection rule (Ib)]. Adding carrier density (µ > 0) favors the singlet Kekulé pairing at sufficiently low temperatures. The high temperature phase is the E u nematic order, favored by the gain in entropy, see Fig. 8b. At the same time, in MLG we find only the E u nematic phase to be the dominant broken symmetry phase for zero and finite µ [selection rule (Ia)]. The adjacent phases then form the following composite order parameters: (3) For quartic interaction in the E g spin-nematic (spin-nematic 1 ) channel M = (Γ 0s01 , Γ 0s32 ). The maximally (but not fully) anticommuting pairing mass is the f -wave pairing. Among the excitonic masses, the CDW and KC orders fully anticommute with M . Note that KC is a mass order in BLG, but commutes with the Dirac Hamiltonian. Accordingly, at µ = 0 or at sufficiently high temperature KC and CDW appear together in BLG, whereas in MLG we only find CDW order in this parameter regime. The f -wave pairing can be observed at µ > 0 and low temperatures in both MLG and BLG, while higher t destroys the superconductivity and the system instead enters into the respective excitonic phases, see Fig. 8c. Note that the E g spin-nematic ordering is completely absent in our range of parameters, thus all appearing broken symmetry phases follow selection rule (Ib). The corresponding composite order parameters constructed from the adjacent phases read (4) For quartic interaction in the E u spin-nematic (or spin-nematic 2 ) channel M = (Γ 3s31 , Γ 3s02 ). These interaction matrices fully anticommute with the s-wave pairing and CDW orders, which indeed appear simultaneously on the phase diagram of MLG and BLG at zero doping, see Fig. 8d. Furthermore, M also anticommutes with the KVBS order parameter, which represents a massive phase for the Dirac fermions, but commutes with the Luttinger Hamiltonian. Accordingly, KVBS accompanies CDW and s-wave pairing in MLG, but not in BLG. As all these phases represent mass orders and they fully anticommute with M , temperature does not lift the degeneracy among them. However, by adding chemical potential we select the pairing phase: s-wave SC. As only selection rule (Ib) is operative in these phase diagrams, we do not obtain any composite order parameter.

C. Quartic interactions: smectic channels
Smectic orders break rotational and translational symmetries, but they do not support any mass gap. Selection rule (Ia) is, therefore, not conducive for a fully gapped phase, and low temperature phases are typically inhabited by mass orders, obeying selection rule (Ib).
1. For quartic interaction in the smectic CDW channel M = (Γ 3010 , Γ 3020 , Γ 3013 , Γ 3023 ). The low temperature region of the phase diagrams in Fig. 9a is occupied by the massive phases that fulfill selection rule (Ib): the swave SC and QSHI in both MLG and BLG. While at zero doping they nucleate in a degenerate fashion, finite doping favors the s-wave pairing at low temperatures. At high temperatures, we observe the smectic CDW phase at both zero and finite µ in MLG and BLG, following selection rule (Ia). The adjacent phases then form the following supervector order parameters: 2. For quartic interaction in the smectic spin density wave (SDW) channel M = (Γ 0s10 , Γ 0s20 , Γ 0s13 , Γ 0s23 ), which fully anticommutes with the QAHI mass order. As such we observe nucleation of this phase in both systems at zero chemical potential for low, as well as high t [selection rule (Ib)]. On the other hand, there is no fully anticommuting available pairing phase, and by elevating the chemical potential we induce the only partially anticommuting smectic SC. Once superconductivity is destroyed, in MLG the system goes back into the QAHI phase. However, in BLG it is the chiral density order parameter that develops a finite expectation value at high t. This is possibly due to the high chemical potential (µ = 0.95), required to induce superconductivity, which causes the system to act like a metal, as opposed to a semimetal. The relevant cuts of the global phase diagram are displayed in Fig. 9b. The corresponding composite order parameters are

D. Quartic interactions: Kekulé channels
Finally, we focus on the quartic interactions in the Kekulé channels. Note that all Kekulé orders (including both bond and current) at least break the translational symmetry, but preserve the rotational invariance. 1. For quartic interaction in the KVBS channel M = (Γ 3011 , Γ 3021 ). The high temperature phase in both systems is KVBS itself, fulfilling selection rule (Ia). Note Low t High t Symmetry Low t High t Symmetry Low t High t Symmetry Low t High t Symmetry  that for the Dirac fermions KVBS is a mass, therefore at µ = 0 in MLG it persists all the way down to the lowest temperature. However, the KVBS order parameter fully commutes withĥ L and in BLG this phase is replaced by the competing massive phases, namely the CDW, AFM and s-wave SC orderings at µ = 0. Since M fully anticommutes with the s-wave pairing order parameter, finite µ selects this phase in MLG as well as BLG, see Fig. 10a. The corresponding supervectors are 2. The quartic interaction in the KC channel displays analogous behavior to that in the KVBS channel, but the mass nature of the phase is flipped between MLG and BLG. Here M = (Γ 0012 , Γ 0022 ) and at high temperatures the dominant instability is KC in both systems [selection rule (Ia)]. In BLG the KC is a massive phase and at µ = 0 it persists down to the lowest temperature. In MLG at zero chemical doping the low temperature regime is occupied by phases that follow selection rule (Ib). At the lowest temperature and µ = 0 we observe degenerate nucleation of CDW, AFM and f -wave pairing states, all representing mass orders with C H = 0. Increasing the temperature the broken symmetry phase becomes the E g nematic, for which C H = A H . See section II A for the definitions of C H and A H . As the KC order parameter fully commutes with the Dirac Hamiltonian (A H = 0), this ordering only occurs at high temperatures. The onset of these orderings along the temperature axis is consistent with selection rule (II). In both systems, finite chemical potential and low temperatures favor the adjacent f -wave pairing, following selection rule (Ib). The corresponding cuts of the phase diagram are displayed in Fig. 10b, and the composite order parameters are 3. For quartic interaction in the spin Kekulé solid (KS) channel M = (Γ 0s11 , Γ 0s21 ). At zero chemical potential, following selection rule (Ib), we observe nucleation of the excitonic CDW order in both MLG and BLG. For µ > 0 and low temperatures the Dirac, as well as the Luttinger fermions condense into the smectic superconductor, whereas intermediate t destroys the pairing phase and gives rise to the E u nematic ordering, for which C H = A H . Furthermore, at even higher temperature BLG displays the spin KS phase, according to selection rule (Ia), that maximally commutes with the Luttinger Hamiltonian, A H = 0. Once again, onset of these phases by increasing temperature is in accordance with selection rule (II). For the phase diagrams in this interaction channel see Fig. 10c. The composite order parameters in these cases are 4. For quartic interaction in the spin KC channel M = (Γ 3s12 , Γ 3s22 ), which fully anticommutes with the s-wave pairing and CDW order parameters. At µ = 0 we observe degenerate nucleation of these two phases for the entire range of temperature. Note that for both order parameters {O i , M j } = 0 and C H = 0 (i.e., they represent fully gapped phases), and thus they are not distinguished by temperature. Only by setting µ > 0 we select the s-wave SC in both MLG and BLG, see Fig. 10d. In the absence of an adjacent high temperature phase we do not obtain composite order parameters in this case. Note that this is due to the fact that only selection rule (Ib) is operative in this case for both MLG and BLG.

VI. HONEYCOMB HUBBARD MODEL
In the preceding sections we scrutinized the honeycomb monolayer and bilayer in the presence of generic short range repulsive electron-electron interactions. We demonstrated through numerous examples the validity of the selection rules and organizing principle (see Sec. II A), by tuning the strength of interactions in one specific channel at a time at zero and finite temperature and chemical doping. However, the bare values of the interaction couplings in a physically realistic microscopic model would in principle be non-zero simultaneously in multiple interaction channels. To substantiate our findings, we here consider one of the simplest yet realistic microscopic models for correlated electrons on the honeycomb lattice, the extended Hubbard model.
In Sec. II B, we present the Hamiltonian containing on-site, NN and NNN repulsion, see Eqs. (6) and (7), while the resulting phase diagrams for MLG and BLG are shown in Figs. 4, 5 and 6. Here we give a detailed walkthrough on how to translate between the Hubbard Hamiltonian, written in terms of a lattice model, and our effective low-energy description. The three distinct interaction terms (on-site, NN and NNN repulsion) are addressed separately in three subsequent sections. In Bernal BLG as the split-off bands are dominantly localized on the a 1 and b 2 sites and here we consider only density-density interactions, one can integrate out fermionic density on these two sites to arrive at an effective or renormalized extended honeycomb Hubbard model defined on the low-energy sites, at the cost of renormalized strengths of the corresponding coupling constants [80]. Such trivial renormalization is, however, omitted here.

A. On-site repulsion (U )
To make a connection with Eq. (7), we first make use of the Hamman decomposition to rewrite H U as [11,108] where n(R) = r † σ (R)r σ (R) is the number operator, m = r † σ (R)σ σσ r σ (R) is the magnetization, with r † σ (R) and r σ (R) being fermionic creation and annihilation operators with spin projection σ =↑, ↓ on the site at R, and summation over repeated indices is assumed. The first (second) term in the Hamman decomposition corresponds to the total   (staggered) density, while the third (fourth) term to the total (staggered) magnetization. Recalling the Fourier expansion of the fermionic fields in Eq. (15) and the spinor structure in Eq. (16), they are written as where m s (R) is the s component of the magnetization (s = 1, 2, 3) on the site at R.
Upon squaring the above quantities we neglect the oscillatory contributions, as any position (and thus wave-number) dependent term will be less relevant in the RG sense (as can be shown via power counting). Then the on-site repulsion in the continuum theory takes the form Hence the initial condition (or bare value) of the coupling constants as a function of U is We then use the Fierz constraints to rewrite the spin triplet four fermion terms as a linear combinations of spin singlet terms, eliminating g t i , and subsequently rescale according to U a 2 /16 → U . Doing so we arrive at and all other linearly independent coupling constants g s i ( = 0) = 0 for i = 2, 3, 4, 6, 7, 8.
To study the phase diagram of the honeycomb Hubbard model in the presence of on-site repulsion we apply the same methodology as in Sec. V. But instead of scanning along one certain axis in the space of coupling constants, we use Eq. (61) and tune U until we detect a phase transition. The (U, t) cuts of the phase diagram for various values of µ are displayed in Figs. 4(a) and 4(b) for MLG and BLG, respectively. At µ = t = 0 both the Dirac and Luttinger fermions display antiferromagnetic (A t 2u ) ordering, which is to be expected on a half-filled bipartite lattice in the presence of repulsive Hubbard interactions. For finite chemical doping the induced pairing phase for linearly dispersing Dirac fermions is the E g nematic superconductor [ Fig. 4(a)], while for the quadratic Luttinger fermions it is the A 1k singlet Kekulé pairing [ Fig. 4(b)]. Note, the adjacent excitonic and pairing phases in both cases fully anticommute with each other and form O(5) supervector order parameters given by One can immediately recognize the pattern where two adjacent and fully (or partially) anticommuting phases fulfill the selection rules (Ia) and (Ib). Selection rule (Ia) in Eq. (4) requires that the kernel of the interaction term and the bilinear order parameter are the same. Since the quartic interactions are written in the particle-hole basis, the four-fermion term in the AFM channel is the dominant interaction with U > 0. Indeed the phase diagrams with finite repulsive Hubbard interactions accommodate the same ordered phases as those corresponding to the purely antiferromagnetic interaction channel [ Fig. 7d]. Consequently, the resulting composite order parameters V U MLG and V U BLG are the same as V AFM 1 and V AFM 2 in Eq. (68), respectively. As mentioned in Sec. II B, the different pairing orders arising in MLG and BLG are due to the fact that the singlet Kekulé pairing is a gapped phase only in BLG, while in MLG the maximally anticommuting pairing order parameter is the E g nematic superconductor. Such a distinction in the paired state at low temperature solely stems from the differences in the normal state band structures in these two systems, in agreement with selection rule (II).

B. Nearest neighbor repulsion (V1)
Next we turn our focus to the H V1 term in Eq. (7), which models (with V 1 > 0) a finite ranged Coulomb repulsion between electrons on the NN sites. We derive the initial conditions similarly to those for the on-site repulsion. The Hamman decomposition in this case reads In our basis the second term takes the form while the number operator on the A sites can be written as We again neglect the oscillatory terms when taking the square of the density and the staggered density, and arrive at the expression After rescaling as 3V 1 a 2 /4 → V 1 we obtain for the bare values of the coupling constants to be Following analogous logic to the on-site repulsion case, we find that the quartic term in the CDW channel is the dominant interaction that fulfills the selection rule (Ia), while the adjacent pairing orders are the s-wave and fwave superconductors, satisfying the selection rule (Ib). This can be anchored by comparing V V1 1 and V V1 2 with Eq. (43), as well as Figs. 5(a) and 5(b) with Fig. 7b. Note the appearing pairing phases fully gap both systems, and therefore we do not see an order differentiation between MLG and BLG, unlike the situation for the on-site repulsion. Below the transition temperature, we expect the pure superconducting state to display either s + if or f + is symmetry that maximally gaps the underlying Fermi surface [81].

C. Next-nearest neighbor repulsion (V2)
Finally, we examine the effect of the H V2 term in Eq. (7), that describes NNN repulsion. Each site on the honeycomb lattice has six NNNs, and the repulsion in this case acts between sites belonging to the same sublattice. The six vectors separating NNNs are ±v i , i = 1, 2, 3, where v 1 and v 2 are the primitive lattice vectors defined in Sec. III A, and v 3 = v 2 −v 1 [Fig. 2]. Once again we decompose the quartic terms and write The first factor of the product contains the density and staggered density, the form of which we already derived in Eqs. (55) and (56), respectively, while the two terms in the second factor can be written as When performing the multiplication we again neglect any oscillatory terms, and write the continuum Hamiltonian as Rescaling 3V 2 a 2 /4 → V 2 , the initial conditions of NNN repulsion in the continuum formalism read g s 1 ( = 0) = g s 5 ( = 0) = −2V 2 , g s 9 ( = 0) = and g s i ( = 0) = 0 for i = 2, 3, 4, 6, 7, 8. At zero temperature and chemical potential we find quantum spin Hall insulator (A t 2g ) phase in both MLG and BLG, see Figs. 6(a) and 6(b) respectively. At finite chemical doping and sufficiently low temperatures, the adjacent pairing phase in both systems is the s-wave superconductor. These phases constitute the composite order parameter Note that both nucleating orders gap the Dirac, as well as the Luttinger fermions. The dominant interaction channel in the presence of NNN repulsion is the QSHI, and indeed the phase diagrams are analogous to Fig. 7c, while the composite order parameter V V2 is the same as V QSHI 1 in Eq. (44). Note, however, in contrast to the pure QSHI interaction channel, repulsive NNN interactions do not support singlet Kekulé pairing in BLG. This is likely due to the smectic charge density wave component in the corresponding initial conditions (g s 9 ).

VII. SUMMARY AND DISCUSSION
Here we analyze a system of spinful fermions on monolayer and bilayer graphene, interacting via short range (or momentum-independent) Coulomb repulsion. We start by writing down the respective lattice models featuring tightly bound electrons and systematically deriving the low-energy continuum theories, that describe two copies of linear and biquadratic band touching, respectively. Retaining the Fourier components in the vicinity of the band touching points, we arrive at analogous descriptions for MLG and BLG, featuring the same number of degrees of freedom. Incorporating the contact interactions then involves constructing all quartic terms allowed by symmetry. As the point group describing the symmetry transformations of the lattice models is the same D 3d group for MLG and BLG, the interacting Lagrangian in Eq. (32) is identical in the two systems. We address the effect of short range electronic interactions in a perturbative fashion, via Wilsonian momentum-shell RG analysis and the expansion scheme up to the one-loop order.
The central distinction in the continuum theories lies in the band structure, and the fact that the band touching points are linear (quadratic) in MLG (BLG). The purely linear (quadratic) dispersions result in the dynamic scaling exponent z = 1 (2) for Dirac (Luttinger) fermions at the non-interacting fixed point, and hence = d − z = 1 (0) in the expansion. Since the coupling constants of contact interactions scale as such interactions are irrelevant (marginal) in MLG (BLG). However, this only sets the physical value of in the RG scheme and therefore the location of the phase boundaries. More importantly, the k-linear functions (k x and k y ) describing Dirac fermions in MLG are odd under inversion and transform under the E u representation of the D 3d point group. In comparison, the d-wave harmonics (k 2 x − k 2 y and 2k x k y ) are even under inversion and thus transform under the E g representation. To obtain an A 1g quantity, they are multiplied by matrices from the same irreducible representation in the Hamiltonian. This algebraic difference is the origin of order differentiation between MLG and BLG, as in the case for the purely on-site Hubbard repulsion, which brings us to the main objective of this paper.
In Sec. II A we outline a set of selection rules proposed in Ref. [7], and argue that while the selection rules (Ia) and (Ib) in Eq. (4) allow for two distinct ways for an ordered phase to be promoted via a certain contact interaction term, the selection rule (II) in Eq. (5) is a generalized energy-entropy argument and organizes these phases along the temperature axis. The example of MLG and BLG is specifically suited to demonstrate these principles, since only (II) relies on the underlying band structure, while (I) is oblivious to it. Therefore, the ultimate ordered phases for a given interaction can in principle be different in MLG and BLG. We demonstrate the validity of the selection rules and organizing principle by increasing the strength of interactions in the individual channels and identifying the nature of symmetry breaking, see Tables III and IV. By adding chemical doping and forming an extended Fermi surface, we also obtain superconductivity from repulsive electronic interactions. We show that the nucleating pairing phases also obey the selection rules, and that the adjacent excitonic and pairing phases consequently constitute composite order parameters, that form an enlarged O(N) algebra.
Besides the individual interaction channels, we also consider a microscopic description of interacting electrons on the honeycomb lattice, the extended Hubbard model, containing the on-site, the NN and the NNN components of the Coulomb repulsion. See Figs. 4, 5 and 6 for the corresponding phase diagrams. Even though in this case we increase the strength of interactions in multiple channels at once, the dominantly diverging channel plays the role of M in the selection rules (Ia) and (Ib). Then, our findings for the extended Hubbard model are in complete agreement with the selection rules. Namely, for repulsive Hubbard U at half filling (µ = 0) we find antiferromagnetic orderings in both systems, while chemical doping (µ > 0) gives rise to E g nematic pairing in MLG and A 1k singlet Kekulé superconductor in BLG. Either of these two pairings forms O(5) composite order parameter with the antiferromagnet [see Eq. (62)]. In the presence of the NN repulsion (V 1 > 0), we find charge density wave ordering when µ = 0 and simultaneous nucleation of both s-wave and f -wave pairing for µ > 0 in both systems. The s-wave and f -wave pairings respectively form This observation is in agreement with Ref. [83], where the authors carried out a non-perturbative functional RG analysis for the NN repulsive interaction. Note that only the f -wave pairing was found in this study, which, however, considered spinless fermions in MLG that does not permit s-wave pairing, due to the requisite antisymmetry property of the electronic wave function. On the other hand, a quantum Monte Carlo study in Ref. [84] found quantum spin Hall insulator and s-wave superconductor phases, respectively at zero and finite doping. The dominant four-fermion interaction in this case is expected to be in the QSHI channel after integrating out the bosonic degrees of freedom. See Fig. 7c.
In light of these qualitative agreements with nonperturbative numerical results, and considering their purely algebraic nature (based on (anti)commutation among matrices), we believe that the validity of our proposed selection rules goes beyond one-loop RG calculations. In Appendix A we show how they manifest in a one-loop RG calculation at zero and finite temperature and zero chemical doping. Possible future investigations could target systems with existing experimental and/or numerical results, e.g. twisted bilayer graphene [109][110][111][112], Weyl semimetals [113][114][115], nodal-loop semimetals [116,117], to name a few.
We close the discussion by qualitatively addressing the role of the trigonal warping [stemming from t BA in Eq. (12)] on the phase diagrams of interacting BLG. Note that trigonal warping splits the biquadratic band touching near each valley into four Dirac points, yielding linearly vanishing DOS, namely ρ(E) ∼ |E| at the lowest energy [9]. Consequently, all orderings in BLG strictly take place at finite coupling, which is, however, much smaller than those in MLG, as t BA t. Recently, it has been shown that even when one begins with purely biquadratic band touchings in BLG, quantum fluctuations via self-energy corrections can split them into Dirac points (analogous to the trigonal warping), thereby deferring onset of any ordered phase to finite coupling [118,119]. Even then our conclusions remain qualitatively valid for BLG. For example, at zero doping when the interaction strength is stronger (weaker) than the scale of the trigonal warping, the ordered phases are determined by the selection rules for Luttinger (Dirac) fermions. At finite doping (setting an infrared cutoff for the RG flow), when the chemical doping is above (below) the scale of trigonal warping, the emergent superconducting phases at the lowest temperature are qualitatively similar to the ones we show for Luttinger (Dirac) fermions. As such for repulsive on-site Hubbard model, a doped BLG is expected to support singlet E g nematic (Kekulé) pairing for small (large) doping at the lowest temperature [13].
This anticipation can be germane (at least qualitatively) in twisted BLG near the magic angle, where a cascade of ordered states (including both excitonic and superconducting) have been observed by tuning the filling of the valley and spin degenerate nearly flat bands [109][110][111][112]. In particular, as the system approaches the magic angle, the Fermi velocity of slow Dirac fermions in twisted bilayer graphene gradually vanishes [120][121][122][123][124][125], that in turn enhances the importance of the higher gradient terms, similar to the ones we considered here for Bernal BLG [124,125], which may therefore play crucial roles in determining the ultimate nature of the ordered phases. In the future, our formalism can be extended to shed light on this enigmatic correlated system. In the above outlined scenario there are two sources of one-loop correction to an order parameter field, the bubble and vertex diagrams [see Fig. 11], denoted by B(t) and V (t), respectively, where t is the dimensionless temperature [see Eq. (33)]. The contribution from the bubble diagram reads where Performing the Matsubara sum and integrating over the momentum shell Λe − < |k| < Λ, we obtain  where dim(M ) is the dimension of M , and The functions f E (t) and f S (t), as well as their arithmetic average are shown in Fig. 12. Note that f E (0) = 1 and f S (0) = 0, and thus the energetically favorable fully gapped phases dominate at low temperatures. On the other hand, at higher temperatures f E (t) decreases and f S (t) increases, hence giving way to phases with higher entropy gain, a manifestation of the organizing principle (II). In this context recall that while {O, H} = 0 leads to fully and isotropically gapped spectra in the ordered phase, [O, H] = 0 yields gapless quasiparticles therein. A partially anticommuting (A H = C H ) order parameter acquires a contribution that is the average of the former two cases. The vertex correction yields the following contribution where g(k) is defined in Eq. (A2). After performing the Matsubara sum and the shell integral this contribution reads where the positive (negative) sign follows if O and M mutually anticommute (commute).
The above results at zero temperature can be summarized as follows.
(a) An order parameter O that fully commutes with the Hamiltonian acquires identically zero anomalous dimension, irrespective of the interaction channel. These findings are also summarized in Table V. Note since k ∼ |k| z the quantity Λ d / Λ ∼ Λ d−z = Λ . Indeed at the lower critical dimension where = 0 the cutoff dependence drops out from B(t) and V (t).

Appendix B: Fierz reduction of local interaction terms
In this appendix we present the Fierz reduction of 18 momentum-independent quartic terms, that describe all the local quartic interactions on monolayer and bilayer graphene. We report the existence of 9 linearly independent four fermion terms, which we choose, without the loss of generality, to be the spin singlet interactions. Since we write all quartic terms in the particle-hole basis, we here refrain from using the Nambu basis, and Ψ and Ψ † are eight-component spinors structured according to Eq. (16).
We start by organizing the interaction terms into a vector as X = X s , X t , where and X s and X t contain the spin singlet and triplet quartic terms, respectively. The Fierz identity for quartic terms of eight-component Grassmann variables reads [10,80] [ where M and N are eight-dimensional Hermitian matrices, Γ a span a basis in the space of such matrices, and the summation over repeated indices is assumed. Here we apply the above identity for contact interactions, for which x = y. Using this identity we write each element in X as a linear combination of the remaining ones. We can summarize these linear connections with the matrix equation F X = 0, with the eighteen-dimensional Fierz matrix (B4) The number of independent quartic terms is given by 18 (dimension of F ) -9 (rank of F )=9. We can extract the relevant equations by reordering the columns of F and performing row reduction. Choosing the spin singlet interactions as the independent ones we obtain nine equations which we can write compactly in matrix form as However, when 9 linearly independent couplings are chosen in the spin singlet channels, we do not generate any quartic term in the spin triplet channel via coarse grain. Nevertheless, the above linear relations allow us to set the initial condition for any quartic interaction in the spin triplet channel in the RG equations, even though they are expressed in terms of the quartic interactions in the spin singlet channel, and that way construct various cuts of the phase diagrams shown in Figs. 7-10.

Appendix C: Renormalization group (RG) flow equations
Here we provide further details to the RG analysis at finite temperature and chemical doping. Finite temperature causes the previously continuous Matsubara frequency to take on discrete values, ω → ω n = (2n + 1)πT , with n ∈ Z (integers) for fermions, where T is the dimensionfull temperature. Using the compact notations Ω ± = iω ± µ, the four distinct loop integrals read Λf ω (t, µ)l, where 2 k = v 2 1 (k) + v 2 2 (k), with v i (k) = p i (k) and d i (k) for Dirac and Luttinger fermions, respectively. See the definitions of these cubic harmonics in Sec. IV A of the main text, and Λ = |k|=Λ is the given quasiparticle spectrum evaluated at the ultraviolet momentum cutoff Λ. The functions of dimensionless temperature and chemical potential [see Eq. (33)] are of the form . (C2) Note, in the limit of zero chemical doping (but at finite temperature) they become On the other hand, in the limit of both zero temperature and zero chemical doping the functions take on the values In the rest of this appendix we display the flow equations of various coupling constants and order parameter conjugate fields. Let us recall the schematic form of them from Eqs. (35) and (41) in the main text, To differentiate between the two systems, we use the notation H i jk (t, µ) = M i jk (t, µ) and B i jk (t, µ) for MLG and BLG respectively. These functions are obtained by computing the Feynman diagrams (c)-(f) in Fig. 13. All the terms that are not explicitly shown are zero. In case of the conjugate fields we write out the full β functions,β D ∆ k (β L ∆ k ) for MLG (BLG), which are obtained by computing the Feynman diagrams from Fig. 11. For brevity, we suppress the t and µ dependence of the functions from Eq. (C2).