{'fraction_non_alphanumeric': 0.032147228724203836, 'fraction_numerical': 0.012541069856336032, 'mean_word_length': 4.431303942150688, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Interacting electrons on trilayer honeycomb lattices

Few-layer graphene systems come in various stacking orders. Considering tight-binding models for electrons on stacked honeycomb layers, this gives rise to a variety of low-energy band structures near the charge neutrality point. Depending on the stacking order these band structures enhance or reduce the role of electron-electron interactions. Here, we investigate the instabilities of interacting electrons on honeycomb multilayers with a focus on trilayers with ABA and ABC stackings theoretically by means of the functional renormalization group. We find different types of competing instabilities and identify the leading ordering tendencies in the different regions of the phase diagram for a range of local and non-local short-ranged interactions. The dominant instabilities turn out to be toward an antiferromagnetic spin-density wave (SDW), a charge density wave and toward quantum spin Hall (QSH) order. Ab-initio values for the interaction parameters put the systems at the border between SDW and QSH regimes. Furthermore, we discuss the energy scales for the interaction-induced gaps of this model study and put them into context with the scales for single-layer and Bernal-stacked bilayer honeycomb lattices. This yields a comprehensive picture of the possible interaction-induced ground states of few-layer graphene.


I. INTRODUCTION
In the last years, the electronic properties of fewlayer graphene systems have been studied extensively 1 . Experiments on graphene systems with two 2-9 and three layers [10][11][12][13][14][15][16][17] aimed to clarify the role of many-body interactions 18,19 . These are expected to be more effective due to the enhanced low-energy density of states and thus correlated electronic ground states may be realized. For bilayer graphene the experimental studies agree on the correlated nature of the ground state, but disagree on the symmetries of the underlying order, a topic that is matter of current debate. At temperatures below 5K, transport experiments found gap openings of about 2-3 meV [7][8][9] . For trilayer graphene the stacking order is crucial for the electronic properties [20][21][22][23][24][25][26][27][28][29] . In Ref. 19, it was found that neutral trilayer graphene with ABC stacking shows many-body correlations with a pronounced gap of ∼ 6 meV while graphene trilayers with ABA stacking do not show a gap. In view of this experimental situation a better theoretical understanding of the many-body instabilities and the nature of the correlated phases which are candidates for possible trilayer graphene ground states is required.
Here, we employ a functional renormalization group (fRG) approach (for a recent review, see Ref. 30) to address the problem of competing instabilities on honeycomb trilayers in an unbiased way. We explore a region of the phase diagram with a range of interaction parameters with density-density repulsions up to the second nearest neighbor that are motivated by the ab-initio values as proposed in Ref. 31. In a previous study, our RG approach was applied to the AB stacked honeycomb bilayer 32 to investigate the phase diagram as a function of local and non-local interaction parameters and complemented by a more quantitative study for pure onsite interactions combining mean-field theory, functional RG and quantum Monte Carlo (QMC) techniques 33 . These studies clarified the relevance of various phases on the bilayer that have been subject to previous theoretical studies by different authors [34][35][36][37][38][39][40][41][42][43][44][45] , e.g. the antiferromagnetic spin density wave (AF-SDW), two kinds of charge density waves (CDW, CDW 3 ) and a quantum spin Hall state (QSH), depending on the model parameters when the full four-band model with the available ab-initio estimates for the model parameters is used. The combination with QMC simulations 33 showed the compatibility of fRG and QMC methods for the investigated range of parameters, and allowed to gauge the quantitative precision of the fRG data at least partially.
The instability toward an interaction-induced quantum spin Hall (QSH) state is of particular interest in connection with realizing topological insulators in graphene. While the original proposal by Kane and Mele 46 , extending earlier work by Haldane 47 with a mass term due to spin-orbit interactions turned out to too small gaps, it was later argued 48 that second-nearest neighbor repulsions could lead to a interaction-induced mean-field having the same effect as the Kane-Mele mass term. In the single graphene layer, however, a nonzero interaction strength is needed to open any kind of gap. So far there is no experimental evidence for this. In the bilayer system, the same instability occurs for arbitrarily small interactions of appropriate distance dependence (see e.g. Ref. 32), and one basically gets two copies of the Kane-Mele QSH state coupled by the interlayer terms. The sign of the order parameter gives rise to distinct choices with interestingly different properties 49 . While this is an interesting many-body state, the corresponding spin-polarized edge modes would not be topologically protected in a strict sense as the bilayer-doubling now permits time-reversal invariant single-particle terms that would localize the edge states 50 . However, if one finds the same state in the trilayer system, one would again have an odd number of Kramers pairs at the edges or edge states per spin, and the topological protection would keep alive at least one edge mode. In this work we show that the fRG suggests that the trilayer QSH state is not unlikely for realistic parameters, and that it can even occur at sizable energy scales ∼ 10 meV.
This paper is organized as follows. In Sec. II, we introduce simple tight-binding Hamiltonians for ABC and ABA trilayer honeycomb lattices and discuss the resulting band structures depending on the stacking order. Further, we introduce the interacting part of the Hamiltonian with onsite, nearest and next-to-nearest neighbor interaction terms. In Sec. III we describe the fRG method and discuss the patching scheme as well as the employed approximations, together with the relations to other methods. In Sec. IV we analyze the emergent effective interactions from the fRG flow. We classify the instabilities and ordering tendencies in ABA and ABC trilayer honeycomb lattices for interacting electrons. With this information, we can construct the tentative phase diagram from fRG in the space of interaction parameters. We finish with a discussion of the energy scales for the gaps in the electronic spectrum due to the ordering phenomena in Sec. IV D and some conclusions in Sec. V.

A. Basic Hamiltonians
Here, we construct the non-interacting part of the tight-binding Hamiltonians for ABC and ABA stacked trilayer honeycomb lattices. The position vectors of the two-dimensional bipartite lattice structure are called R, with each R having three nearest neighbors at positions that can be characterized by the three vectors δ n with n ∈ {1, 2, 3}. Explicitly, the δ n are given by 2 e x + a 2 e y and δ 3 = −a e y . Here, a is the distance between two neighboring lattice sites and the vectors point from the B-sublattice to the A-sublattice. The tight-binding Hamiltonian for ABC/ABA stacked trilayers is composed out of single layer Hamiltonians for the in-plane hoppings, perpendicular hoppings between different layers and remote hoppings between the layers including a certain planar distance. We denote the different lattice sites by A 1 , B 1 , A 2 , B 2 , A 3 , B 3 , where A and B specify the sublattice and the index numbers the layer as shown in Fig. 1. Accordingly, we define the annihilation operators c o,s, r and the creation operators c † o,s, r of an electron at position r and the layer and sublattice are encoded in the subscript o ∈ {a 1 , b 1 , a 2 , b 2 , a 3 , b 3 }.  The index s =↑, ↓ specifies the electron spin. With these preliminaries we can write down the single-layer tightbinding Hamiltonians in the first two layers, and a stacking dependent Hamiltonian for the third layer, where we have attributed an additional index c to H 3,c , with c = 0 for the ABA stacking and c = 1 for the ABC or chiral stacking. In graphene and few-layer graphene systems the hopping γ 0 = t ≈ 3 eV, see e.g. Refs. 40,51.
Next we introduce the interlayer hoppings, γ 1 , between sites that lie on top of each other and connect adjacent layers, where again we accounted for the different stackings by introducing H ⊥ 23,c . Ab initio values for γ 1 = t ⊥ are available for graphite 26,52 , t ⊥ ≈ 0.4 eV , and ABC trilayer graphene 40 , t ⊥ ≈ 0.5 eV. In Fig. 1, we also show the more remote hoppings γ 2 , γ 3 , γ 4 , γ 5 whose effect will be discussed in the following section. In the present model study we will ignore these terms for most explicit calculations. For the discussion of the energy bands, we introduce the Fourier transform where N is the number of unit cells and k is an element of the first Brillouin zone (BZ). We can write the tightbinding Hamiltonians in this approximation in Fourier space as ) and where d k = n exp(i k · δ n ). The resulting energy bands are depicted in Fig. 2. The dispersion for ABC trilayers is very flat near the Fermi level with cubic wavevector dependence close to the K and K ′ points. Hence one should expect an enhanced role of interactions as compared to ABA trilayers and AB bilayers where the wavevector dependence is quadratic.

B. Remote hoppings in trilayers
In this section, we want to discuss the effect of remote hopping terms on the band dispersion and possible implications for instabilities. For the ABC trilayer system the density of states close to the Fermi energy has a van Hove singularity ∼ ǫ −1/3 due to the cubic band crossing point. This enhances the role of interaction effects, leading to high critical scales for ordering tendencies. In our model study, we therefore choose to ignore remote hoppings that change the topology of the bands below an energy of order 10 meV 40 , which is of the order of the measured gap in the experiment. Of course, it is an interesting question whether these terms affect the nature of the instability in trilayer graphene. However, we find in our computations that instabilities occur already on an energy scale that is higher than or at least comparable to 0.004t ∼ 10 meV. This serves as an a posteriori justification for dropping the remote hopping terms and proves the consistency of our approach.
In ABA trilayer the situation is different. Here, the band structure close to the K-, K ′ -points is separated into a J = 1 (linear) and a J = 2 (quadratic) subband. The remote hoppings induce separate gaps for these subbands, however, with an individual energy shift of the subbands. This destroys the particle-hole nesting and the associated instabilities are suppressed (at least for small interaction terms). Explicitly, the tightbinding Hamiltonian in presence of the important remote hoppings 26,27,52 , γ 2 = −0.02 eV, γ 3 = 0.3 eV, γ 4 = 0.04 eV, γ 5 = 0.04 eV, δ = 0.05 eV, can be written aŝ (10) Within the fRG approach, we will show that also for this system flows toward many-body instabilities occur, however, only beyond a critical interaction strength which, similarly to the single-layer case, most probably is larger than the interaction strength in the real material.

C. Interaction terms
In order to investigate the instabilities that are possible on the trilayer honeycomb lattice for interacting electrons we will take into account a number of different interaction terms, most importantly an onsite or Hubbard interaction U , a nearest neighbor intralayer interaction V 1 and a next-to-nearest neighbor intralayer interaction V 2 . For these interaction parameters ab initio parameters from constrained random phase approximation (cRPA) computations are available 31 and we take those values as a motivation for the investigated range in the phase diagram. The interaction Hamiltonian reads is performed numerically and has to be carried out on H I as well. This adds 'orbital makeup' to the interaction terms in band representation, leading to an additional angular dependence of the interactions near the K, K ′ -points. In order to resolve this dependence we use an angular patching of the interaction terms as explained in Sec. III.

III. FRG METHOD
We employ a functional renormalization group (fRG) approach for the one-particle-irreducible vertices with a momentum cutoff (for a recent review on the fRG method, see 30 ). In this scheme, an infrared regulator with energy scale Λ is introduced into the bare propagator function, Here, ω is the Matsubara frequency, k the wavevector, b the band index and ǫ( k, b) is the single-particle dispersion. As the spin is conserved the free propagator is spin independent. The cutoff function is chosen to enforce an energy cutoff, which regularizes the free Green's function by suppressing the modes with band energy below the scale Λ, For better numerical feasibility the step function is slightly softened in the actual implementation. With this modified scale-dependent propagator, we can define the scale-dependent effective action Γ Λ as the Legendre transform of the generating functional G Λ for correlation functions, cf. Ref. 53. The RG flow of Γ Λ is generated upon variation of Λ. By integrating the flow down from an initial scale Λ 0 , which is in our case chosen as the maximum energy of all bands, to the infrared Λ → 0, one smoothly interpolates between the bare action of the system and the effective action at low energy. In order to limit the numerical effort we employ the following approximations. First, the hierarchy of flowing vertex functions is truncated after the four-point (twoparticle interaction) vertex. Secondly we neglect the frequency dependence, by setting all external frequencies to zero, as we are interested in ground-state properties. The general coupling function, which depends on three momentum and Matsubara frequency indices (the fourth index is fixed by conservation) is thus replaced by the coupling function Third, we neglect self-energy corrections. This approximate fRG scheme then amounts to an infinite-order summation of one-loop particle-particle and particle-hole terms of second order in the effective interactions. It allows for an unbiased investigation of the competition between various correlations, by analyzing the components of V Λ (k 1 ,k 2 ,k 3 , o 4 ) that create instabilities by growing large at a critical scale Λ c 30 . With the approximations mentioned above, this procedure is well-controlled for small interactions. At intermediate interaction strengths we still expect to obtain reasonable results. In any case, the fRG takes into account effects beyond mean-field and random phase approximations.
The fRG calculation is performed in the band basis in which the free part of the Hamiltonian is diagonal. The free Hamiltonian can be diagonalized by a unitary transformation of the form where the index o includes all six sublattices in the three layers and the index b denotes the corresponding bands. In the interaction part of the Hamiltonian the coupling function V Λ (k 1 ,k 2 ,k 3 , o 4 ) in orbital space is chosen so that it reproduces the interaction Hamiltonian of Eq. (11). For application of the fRG in the band basis we have to transform the interaction via Thus, in the band basis the interaction vertex V (k 1 , k 2 , k 3 , b 4 ) acquires a pronounced momentum dependence due to the extra prefactor u o1b1, k1 u o2b2, k2 u * o3b3, k3 u * o4b4, k4 , which is sometimes also referred to as 'orbital makeup'. Previous studies have shown that this has a crucial impact on the fRG flow by allowing for unconventional instabilities 54 .
With the approximations discussed above the fRG flow for the coupling function reads with the particle-particle channel  ). Below, the loop contributions to the flow of the interaction vertex including the particleparticle-diagram (a), the crossed particle-hole-diagram (b) and the direct particle-hole diagrams (c). Right panel: Patching of the Brillouin zone as explained in the text.
the direct particle-hole channel and the crossed particle-hole channel where k = ( k, b) collects the wave vector and the band index. As mentioned above, external frequencies are set to zero The frequency of the second line is fixed by frequency conservation to be −ω in the particle-particle diagram and ω in the particlehole diagrams. V BZ ist the volume of the BZ. The loop kernel is given by where in our approximation self-energy corrections are neglected, i.e., the full propagator is identical to the free propagator.
The wavevector dependence of the interaction vertex is simplified by discretization. The BZ is divided into N patches with constant wavevector dependence within one patch, so that the coupling function has to be calculated for only one representative momentum in each patch. The representative momenta for the patches are chosen to lie close to the Fermi level. The patching scheme is shown in Fig. 3, with N = 24. Each of the four momenta in V Λ (k 1 , k 2 , k 3 , b 4 ) is additionally equipped with a band index. Momentum conservation fixes one of the four wavevectors. Altogether this results in a 6 4 × N 3 component coupling function V Λ .
We start the fRG flow at the initial scale Λ 0 which is in our case chosen as the maximum energy of all bands. We then integrate out all modes of these bands by decreasing Λ. In typical flows some components of the effective interaction vertex become large and diverge at a critical scale Λ c > 0. In this work we use the scale at which the interaction vertex exceeds a value of the order of 10 times the bandwidth as an estimate for the critical scale. The precise choice of this value has only a minor effect on the extracted critical scale, as the couplings grow very fast in the vicinity of the divergence.
The divergence is strictly speaking a (physically meaningful) artifact caused by the neglect of the self-energy in the flow. With self-energy correction a gap would open up or some other modification of the low-energy spectrum would take place, and the flow would be regularized. This is all well known from the Cooper instability in superconductors. Our analysis here tells us in which channel ordering occurs most prominently. The pronounced momentum structure of the vertex near the critical scale can be used to extract an effective Hamiltonian for the low-energy degrees of freedom. This is used to determine the leading order parameter of a given instability. Furthermore, the scale Λ c can be interpreted as an estimate for ordering temperatures, if ordering is allowed by the Mermin-Wagner theorem, or at least as the temperature below which the dominant correlations should be clearly observable. Furthermore, one can understand Λ c as energy scale for the modification of the spectrum, typically by a gap.
In this work we study the flow at temperature T = 0. We find flows to strong coupling with non-zero critical scales Λ c for all choices of non-vanishing interaction terms provided there is a non-vanishing density of states at the Fermi level of the coupled layers.

IV. INSTABILITIES AND PHASE DIAGRAM
A. ABA and ABC trilayer Hubbard model Let us start the description of the fRG results with the case of onsite interactions only, i.e. U > 0, V 1 = V 2 = 0. We limit the study to the charge-neutrality point, i.e. with Fermi points at K and K ′ in the Brillouin zone.  Fig. 3. On the horizontal axis the wavevector k1 can be read off and on the vertical axis we enumerate k2. k3 is fixed on the first patch, k4 then follows from momentum conservation. Here, we see that the sharp vertical structure (k1 = k3) comes with the double magnitude as the horizontal structure (k2 = k3). Middle Panel: Effective vertex function for the orbital combination, where o1 = o3, o2 = o4 = o1 and if o1 ∈ {a1, a2, a3} then o2 ∈ {b1, b2, b3}. Right panel: Effective vertex function for the orbital combination, where o1 = o4, o2 = o3 = o1 and if o1 ∈ {a1, a2, a3} then o2 ∈ {b1, b2, b3}. If in the second and third plot both o1, o2 ∈ {a1, a2, a3} or ∈ {b1, b2, b3} the sign of the vertices changes.

Simplified band structures
First, we want to investigate the simpler band structures when all remote hopping terms are neglected and only t = 0 and t ⊥ = 0. Then, running the fRG flow with pure onsite interaction U = 0 in the ABC and ABA stacked trilayer, we observe an instability toward an antiferromagnetic spin density wave (AF-SDW) with a typical signature of the interaction vertex near the instability as shown in Fig. 4. The leading part of effective interaction corresponding to the clearly discernible sharp structures in wave vector space reads in this case with V oo ′ > 0 and S o q=0 = 1 2 k,s,s ′ σ ss ′ c † o,s, k c o,s ′ , k . The fact that the above Hamiltonian only contains the q = 0 component means that the effective interaction has become infinitely-ranged 55 . The factors ǫ o depend on the orbital, ǫ o = +1 for o ∈ {a 1 , a 2 , a 3 } and ǫ o = −1 for o ∈ {b 1 , b 2 , b 3 }. This sign structure implements the staggering of the interaction within the unit cell appropriate for antiferromagnetic interactions. Note that this parametrization holds in both cases, for the ABC as well as for the ABA stacking.
A mean-field decoupling of H AF results in an AF spin alignment in each layer where a net spin (e.g. 'up') moment is located on the A 1 -, A 2 -and A 3 sublattices, and an opposite net spin ('down') moment on the B 1 -, B 2and B 3 -sublattices.
The critical scale Λ c as a function of the onsite interaction U for ABC and ABA honeycomb trilayers with model hopping parameters t = t ⊥ is shown in Fig. 5 together with the critical scales of single-and bilayer hon- eycomb lattices with the same hopping parameters. This choice of band parameters takes us beyond the regime of realistic parameters for TLG, but pronounces the characteristic features of the ABC and the ABA stacking close to the K, K ′ points and therefore allows to study the differences of the various honeycomb stacks more explicitly. Furthermore, it allows to compare to recent QMC results for the bilayer system 33 . We add a systematic study of the dependence on t ⊥ below. Most importantly, we observe that in the ABC trilayer, the critical scale decreases more slowly as compared to AB bilayer when the onsite interaction is decreased. While in the case of ABA trilayer and AB bilayer with quadratic band crossing points the functional dependence of Λ c (U ) can be fitted by an exponential decay ∼ exp(−α/U ) (cf. inset of Fig. 5), this does not hold for the ABC trilayer case. Instead, at small U , based on the density of states ∼ ǫ −1/3 , one would naively expect a behavior Λ c ∼ U 3 . This is however not reproduced by our data, presumably due to the influence of the high energy sector, i.e. additional bands. We expect that the leading Λ c ∼ U 3 dependence might be recovered at smaller U and thus smaller Λ, which is numerically hard to access due to the rapidly decreasing critical scale. From this analysis we conclude that in ABC trilayers an interaction-induced gap, which we expect to be of the order of the critical scale Λ c , may be considerably larger than in the other structures. This also implies a more stable correlated ground state. All this seems compatible with recent experiments 19 .

Inclusion of remote hoppings in ABA trilayers
When the remote hoppings γ 2 , γ 3 , γ 4 , γ 5 and the onsite energy δ in ABA trilayers are taken into account, the band dispersion is deformed considerably and the particle-hole nesting of the band structure is destroyed 26,27 , see Fig. 2. In this case, one should expect that in the weakly interacting limit the tendency toward instabilities will be strongly reduced. However, it is interesting to study what happens in the case of intermediate and larger onsite interactions.
An accurate treatment of the dispersion with remote hoppings at lowest scales would require the implementation of a new patching scheme to resolve the vicinity of the K and K ′ points, i.e. the fact that the Fermi surface is now not restricted to a single point. This is beyond the scope of this work. Therefore we now stop the flow at the energy scale Λ * at which the bands become non-monotonous. This procedure is routinely done in parquet-and g-ology studies of imperfectly nested bands 56 . Of course it leaves the scales below Λ * unintegrated, but as the dispersion at these lowest scales is not nested and partially gapped, we do not expect significant contributions of these modes to the flow. Hence we expect that the so-obtained estimate for critical scales and U c is already quite good. For the hopping parameters given in Refs. 26,27 the energy scale at which the dispersion is not monotonous any more is given by Λ * ∼ 10 meV ∼ 0.004t.
However, it is important to notice that the remote hoppings already change the dispersion far above the scale Λ * . This explains why for the ABA stacking (here from now on called ABA * stacking when the remote hopping are included), the effect of the remote hoppings is quite drastic. Studying again the case of onsite interactions U only, we observe clearly diverging susceptibilities only for onsite interactions U 2.6t at critical scales well above 0.004t. For smaller interactions, the couplings grow very slowly and no divergences at finite scales above Λ * can be identified. In Fig. 6 we show the fRG results for the critical scales for the ABA * trilayer with remote hopping terms. For comparison, we also show the curves for the single-and bilayer system. As the similarity to the single layer is strong, this analysis suggests a critical onsite interaction U c ≈ 2.6t above which a many-body instability can occur in the ABA * stacking with remote hopping included. We would like to add that the fRG in the present approximation has the tendency to overestimate critical scales. Therefore we would expect the true U c to be slightly larger. For instance, in the single layer system for onsite interactions only, QMC gives U c,QMC ≈ 3.4t for the opening of a single-particle gap 57 , while in the fRG we find U c ≈ 2.6t as well (see Fig. 5).
The experimental study of Ref. 19 does not find a gap for ABA( * ) trilayer graphene. This is compatible with our findings, given the interactions in real material are weaker than this critical value. This is to be expected for consistency, as the critical interaction strength of the single layer and the ABA * -trilayer are roughly the same, and the single-layer remains semi-metallic, too. Even if they are slightly above the threshold, the expected transport gaps would be small and very hard to measure. Therefore, while a precise quantitative picture can- not be obtained with our approximate fRG method and here for onsite interactions only, on a qualitative level we reach consistent conclusions. It would be interesting to resolve better the critical region close to U c,TL for the ABA * stacked trilayer with remote hoppings and analyze the onset of instabilities and their nature. However, this would require a different implementation of our patching scheme which we leave out for future work. In the remainder of this work, we will therefore concentrate on a more thorough study of the ABC trilayer model.

B. ABC trilayer: Instabilities with non-local interactions
From the results for pure onsite interactions and also inspired by recent experiments, we conclude that while ABC trilayers strongly support the formation of manybody states, the ABA * trilayer including remote hoppings most probably does not. We take this as a motivation to study further the instabilities of ABC trilayers for wider range of non-local interaction parameters, which brings us closer to the real materials. Ab initio values for the strength of the density-density interactions up to the third nearest neighbor were listed for single-layer graphene and graphite in Ref. 31. Most likely, one can safely interpolate the parameters for the bi-and trilayer case from this data.
Running the fRG for extended interactions, we find a number of different phases as previously described for the bilayer case in Ref. 32, namely a charge density wave (CDW), a quantum spin Hall state (QSH) and charge density wave with non-zero momentum transfer (CDW 3 ) The numbers on the axis specify the number of the patch as shown in Fig. 3. On the horizontal axis the wavevector k1 can be read off and on the vertical axis we enumerate k2. k3 is fixed on the first patch, k4 then follows from momentum conservation. Left Panel: Effective vertex function for the orbital combination, where o1 = o3, o2 = o4 = o1 and if o1 ∈ {a1, a2, a3} then o2 ∈ {b1, b2, b3}. If in the right plot both o1, o2 ∈ {a1, a2, a3} or ∈ {b1, b2, b3} the sign of the vertices changes.
alongside the AF-SDW. For the investigation of the phase diagram we take into account non-local interaction contributions, namely the nearest-neighbor in-plane interaction V 1 and the second-nearest neighbour in-plane interaction V 2 . More remote interaction contributions are neglected. In the bilayer case 32 we checked explicitly that a third-nearest neighbor repulsion V 3 does not change the picture. Also, we do not consider interlayer interactions. The ab-initio calculations in Ref. 31 showed that the nearest interlayer interactions in graphite (as well as in layered graphene) are of the order of the V 3 -term 58 .
The two leading non-local interaction terms V 1 and V 2 trigger the appearance of qualitatively different instabilities. For dominating V 1 we find an instability toward a charge density wave with a momentum signature of the effective interaction as shown in Fig. 7. This momentum structure can be written down as an effective interaction Hamiltonian of the form, This sign structure supports an enhanced occupancy of the A i sublattices as compared to the B i sublattices or vice versa. Furthermore, a mean-field decoupling of this effective interaction gives a gap in the single-particle spectrum.
A dominating V 2 yields an instability whose dominant interaction terms can be cast into an effective Hamiltonian of the following type, ) cos( The numbers on the axis specify the number of the patch as shown in Fig. 3. On the horizontal axis the wavevector k1 can be read off and on the vertical axis we enumerate k2. k3 is fixed on the first patch, k4 then follows from momentum conservation. Here, we see that the sharp vertical structure (k1 = k3) comes with the double magnitude as the horizontal structure (k2 = k3). Middle Panel: Effective vertex function for the orbital combination, where o1 = o3, o2 = o4 = o1 and if o1 ∈ {a1, a2, a3} then o2 ∈ {b1, b2, b3}. Right panel: Effective vertex function for the orbital combination, where o1 = o4, o2 = o3 = o1 and if o1 ∈ {a1, a2, a3} then o2 ∈ {b1, b2, b3}. If in the second and third plot both o1, o2 ∈ {a1, a2, a3} or ∈ {b1, b2, b3} the sign of the vertices changes .
be decoupled in a purely imaginary 'Kane-Mele' order parameter. This type of instability represents a manybody path to the quantum spin Hall state, where the mass term due to spin-orbit interaction in the original Kane-Mele proposal 46 is now provided by an interactioninduced mean-field. In the single-layer 48 and bilayer 32 honeycomb models, this instability was found in the same corner of interaction parameter space. For an odd number of layers, the mean-field Kane-Mele-ordered state will give rise to an odd number of helical edge modes, and thus represent a two-dimensional interaction-driven topological insulator with protected edge modes. Finally, in the niche of the parameter space for smaller U , we also recover the CDW 3 phase, that we already found in the bilayer system 32 , Fig. 9 for the characteristic momentum structure of the effective interaction. The order parameter due to the symmetry breaking N o Q = 0 is in complete analogy to the one in the honeycomb bilayer 32 , except for the adapted definition of the ǫ o . Within one layer, this order forms three inequivalent sites with different charge densities. The sign structure of the order parameter on different layers and sublattices is determined by the ǫ o -factors, so as to lower the energy contribution from (25). This leaves the total phase of the order parameter undetermined. Depending on this phase, the quadratic mean-field Hamiltonian for a single layer either exhibits a gapless spectrum with Dirac points shifted away from the K, K ′ points or a fully gapped state. Which case represents the mean-field The numbers on the axis specify the number of the patch as shown in Fig. 3. On the horizontal axis the wavevector k1 can be read off and on the vertical axis we enumerate k2. k3 is fixed on the first patch, k4 then follows from momentum conservation. Here, we see that the order parameter has non vanishing momentum transfer as the sharp feature is not at k1 = k3. Right Panel: Effective vertex function for the orbital combination, where o1 = o3, o2 = o4 = o1 and if o1 ∈ {a1, a2, a3} then o2 ∈ {b1, b2, b3}. If in the right plot both o1, o2 ∈ {a1, a2, a3} or ∈ {b1, b2, b3} the sign of the vertices changes.
ground state as function of the interaction parameters has yet to be determined variationally.

C. ABC trilayer phase diagram from fRG
For a systematic investigation of the ABC trilayer phase diagram, we scan a range of interaction parameters U , V 1 and V 2 , whose ab-initio values are listed in Ref. 31. As we expect the fRG to overestimate the critical scales, we take these ab-initio parameters as upper bounds for the range our phase diagrams. In Fig. 10, we show the fRG phase diagram obtained by identifying the leading tendencies in the effective interactions near the instability with an underlying contour plot of the critical scale Λ c in units of t.
We also mark the ab-initio values in the lower plots of Fig. 10, obtained by taking the values from Ref. 31 and scaling {U, V 1 , V 2 } → α{U, V 1 , V 2 } so as to hit the values U = 2t and U = 3t. In both cases, these choices place the system near the phase boundary between QSH and AF-SDW instability. For the more long-ranged single-layer graphene interactions, one finds a QSH state, while for the slightly shorter ranged graphite parameters, one gets a AF-SDW state. Hence, which order occurs might be a delicate issue that is decided by details. In our approximation, the critical scales interpolate smoothly across the phase borders, indicating a weaker competition between the different tendencies. Note that due to this borderline situation there is no true necessity for different layered graphene systems, e.g. with different environments, to exhibit the same ground state, and even the energy scales or gaps could come out similarly despite different states might be selected. Hence, regarding the leading instability, the situation in the ABC trilayer is very similar to the one found previously in the Bernal-stacked bilayer

system.
As mentioned before, we also find a rather exotic density wave phase CDW 3 with a tripling of the unit cell within the layers. This state is subject to further investigation. However it occurs only at quite unrealistic corners of the parameter space with dominant non-local terms. Hence, we do not discuss it further here.

D. Energy scales of the ordering phenomena
In our previous study of the bilayer system 32 , we encountered a problem of the current model studies that in principal also affects the present work. However, we will now also offer an explanation of what happens and how one should read the data in order to get reasonable agreement with and a more quantitative picture of the experiments.
As shown in this work and also in the previous paper on the bilayer system 32 , the simple models employed by virtually all many-body approaches to interaction effects in few-layer graphene can produce very large critical scales. This can be already seen in Fig. 10, where we indicate the values of the critical scales in units of the hopping t ≈ 3 eV. These scales are huge for most choices of the interaction parameters. This is not surprising. In our and other theoretical approaches, the large scales are simply caused by the perfect particle-hole nesting of the : FRG critical scale for ABC trilayer for the pure Hubbard model (blue solid) and rescaled ab initio parameters for graphite (red solid) and single-layer graphene (orange solid) as described in the text with t ⊥ = 0.1t. For comparison we also show the corresponding critical scale for the AB bilayer structure (dashed) and the singlelayer structure (dotted). For the Hubbard (blue curves) and the rescaled graphite parameters (red curves) we find the system to be in the AF-SDW phase, for the rescaled graphene parameters the system is in the QSH state (orange curves). The dashed horizontal line marks 10 meV, which is the order of magnitude, where the topology of the band structure in ABC trilayer graphene changes and remote hoppings become important.
band structure. Furthermore, from the comparison with QMC calculations in the case of pure onsite interactions 33 we know that the overestimate of the fRG is certainly not severe and, expressed conservatively, is less than an order of magnitude in the regime where also the QMC finds robust ordering. As the critical scale is an estimate for the energy scale of the spectral restructuring or gap opening in the ordered phase, a high critical scale would correspond to large energy gaps. If we took the ab-initio parameters of Ref. 31 literally, the theoretical gap estimates would exceed the experimentally observed gap scales ∼ 1 − 10 meV by orders of magnitude.
Let us now analyze the systematics of these critical scales a bit more deeply. In Fig. 11, we show the critical scales obtained from fRG for the single-, bi-and trilayer graphene for three cases, namely for onsite interaction U only, and for nonlocal interactions with the cRPA parameters V 1 and V 2 for graphite and graphene with repulsions up to the second nearest neighbor, for realistic interlayer hopping t ⊥ = 0.1t. The curves show the dependence on the overall magnitude of the onsite interaction, where the ratio between the local U and the non-local interaction parameters is kept fixed. Obviously there are two regimes: a high-scale regime with large critical scales that do not depend too strongly on the interaction strength and a low-scale regime with a very strong dependence on the interaction. For the single layer the second regime is very narrow and basically only contain the minimal critical interaction strength U c (V 1 , V 2 ) below which the semi- metal is stable. Also for bi-and trilayer, the high-scale regime starts above the single-layer U c (V 1 , V 2 ). Next let us consider the dependence of the critical scale Λ c on the perpendicular hopping t ⊥ in the ABC trilayer. Whereas the nature of the ground state qualitatively remains the same for all choices of t ⊥ = 0 the absolute value of this parameter has an impact on the size of the critical scale and shows different behaviors on the two different sides of the critical onsite interaction of single-layer graphene U c,SL , see Fig. 12, i.e. whether we are in the high-scale or in the low-scale regime. For large interactions U > U c,SL the size of t ⊥ is of no importance for the critical scale Λ c . This changes for smaller U . Here, the smaller t ⊥ is, the stronger is the Λ c -variation with the interaction strength.
The comparison with QMC in Ref. 33 was done at larger t ⊥ = t and U ≥ 2.8t, where the scales do not vary that strongly. For band structure parameters with t ⊥ 0.1t and small interactions U U c,SL , we also observe flows to strong coupling for the ABC and the ABA trilayer system, however, the critical scales turn out to be very small, an effect that reflects the behavior of the singlelayer system, where no instabilities occur for U < U c,SL . While this does not constitute a difficulty for the fRG method per se it makes the numerical evaluation very tedious and renders the comparison with QMC impossible.
We now argue that in order to account for the approximations made in the fRG scheme, mainly the neglect of self-energy corrections, and in order to obtain a realistic picture, we have to reduce the cRPA parameters with U ∼ 3t by hand. The reduction factor is chosen so as to obtain the experimentally verified semi-metallic solution for the single-layer case, i.e. roughly to U ∼ 1.5t (if we take the single-layer cRPA values) with an anal-ogous rescaling of the non-local couplings by the same factor. This shift on the interaction axis now takes us from the high scale regime into the regime of strongly varying scales, cf. Fig. 11. Now the critical scales for bi-and trilayer end up in the range 0.01t ∼ 30 meV which is much closer to the experimental values for gap sizes and already in the correct order of magnitude.
In the experiments on bi-and trilayer graphene of these types, the energy gap in the trilayer 19 came out a factor 2-3 higher than in the bilayer. In the fRG the ratio between the critical scales for these two systems depends on the parameter values for the interactions. However for intercations where single-layer graphene does not undergo a phase transition, the larger density of states near the Fermi level of the trilayer case leads to a larger critical scale in theory as well, see Fig. 11. Hence the different energy gaps in bi-and trilayers are qualitatively captured by this theory. Remarkably, with this choice of interaction parameters, the instabilities in ABC trilayer occur on an energy scale where remote hoppings start to be important (∼ 10 meV, see Fig. 11).

V. DISCUSSION
We have performed extensive fRG calculations on honeycomb trilayer systems with different stacking orders, as model systems for trilayer graphene. In doing this we have used as far as possible the available input parameters from ab-initio calculations. Moreover, we have taken into account the full 6-band band structure obtained within the model with one effective p z -type orbital per carbon site. The fRG approach we have used is certainly not exact, but it goes far beyond mean-field studies, random-phase approximation approaches and other perturbative calculations performed for these systems. Furthermore, in the bilayer case with onsite interactions only, the qualitative and quantitative information of the fRG compares quite favorably with the results of QMC calculations 33 . For non-local interactions, QMC encounters severe sign problems, and there is no possibility for benchmarks, but there is also no clear reason why the fRG method should get worse then. Hence, for the given theoretical model the fRG results should be rather robust both qualitatively and also quantitatively (with the adaptations discussed is Sec. IV D and mentioned below) at least as an order-of-magnitude estimate.
First let us mention the main results of the present trilayer study before we get to the connections with graphene systems and experiments. Comparing the different trilayer stackings with the Bernal-stacked bilayer and the monolayer results, we could identify the ABC trilayer as the system that is most prone toward instabilities, with larger energy scales for ordering than the AB-bilayer. The ABA trilayer without the additional interlayer 'remote' hoppings is comparable in its critical scales to the bilayer, but we showed the remote hopping terms with the suggested realistic parameter val-ues (called ABA * structure here) are likely to remove the instability at least for smaller, possibly realistic, interaction strengths. Interestingly, for the Hubbard onsite interaction case, the minimal U -value for obtaining a gapped ground state in the ABA * -structure is close to the one for the single layer. Taking the current experimental knowledge for bi-and trilayer graphene, these theoretical findings regarding the systematic differences are fully consistent with the observations.
Due to the uncertainties about the parameter values for the theoretical model and the approximations made in our approach, we cannot expect a fully quantitative description. However, the phenomenological input of requiring the single-layer to remain semi-metallic puts bounds on the bare interaction parameters that should be used in our model. The idea used in this paper is to scale down the ab-initio parameters for the interactions by an appropriate factor in order to compensate for the approximations, so as to render the fRG flow for the single layer free of divergences. With this remedy for the inexactness of our approach, the ABA * trilayer with remote hoppings remains semi-metallic and the energy scales of the unstable systems, namely the bilayer and the ABC trilayer come out in a quite realistic region below 30 meV.
The type of the leading instability in these trilayer systems can also be read off from the fRG effective interactions near the instability. Here, for the interaction parameters determined by ab-initio methods 31 , our calculations show a strong competition between SDW antiferromagnetic ordering and the interaction-driven 'Kane-Mele' quantum spin Hall state. Which tendency wins depends on the detailed spatial profile of the interactions. A more longer-ranged behavior favors the QSH instability, and for pure onsite repulsions, the AF-SDW state is the clear winner. Both states would open up a bulk gap. The SDW states should have an interesting modulation of the ordered moments depending on the number of nearest neighbors, with smaller moments for higher coordination number 33 . The QSH state should be a true two-dimensional topological insulator for the trilayer case, as time-reversal invariant edge defects will not suffice to gap the three pairs of helical edge states completely, in contrast with the bilayer case (for a discussion of the edge-state robustness, see the review in Ref. 50). Hence, at least one pair of counter moving, spin-resolved helical edge states should survive time-reversal-invariant edge disorder and could hence serve as smoking gun for such a correlated ground state. This perspective is rather exciting, as this would be the first realization of the concept of an interaction-driven topological insulator ('topological Mott insulator') 48 .
We acknowledge discussions with F. Zhang