Functional renormalization of spinless triangular-lattice fermions: $N$-patch vs. truncated-unity scheme

We study competing orders of spinless fermions in the triangular-lattice Hubbard model with nearest-neighbor interaction. We calculate the effective, momentum-resolved two-particle vertex in an unbiased way in terms of the functional renormalization group method and compare two different schemes for the momentum discretization, one based on dividing the Fermi surface into patches and one based on a channel decomposition. We study attractive and repulsive nearest-neighbor interaction and find a competition of pairing and charge instabilities. In the attractive case, a Pomeranchuk instability occurs at Van Hove filling and $f$-wave and $p$-wave pairing emerge when the filling is reduced. In the repulsive case, we obtain a charge density wave at Van Hove filling and extended $p$-wave pairing with reduced filling. The $p$-wave pairing solution is doubly degenerate and can realize chiral $p+ip$ superconductivity with different Chern numbers in the ground state. We discuss implications for strongly correlated spin-orbit coupled hexagonal electron systems such as moir\'e heterostructures.


Introduction
For decades the single-band Hubbard model has been the Standard Model of Correlated Electron Physics. Not only has it been thought of as capturing essential features of the phase diagram of high-temperature superconductors and related materials, it has also served as a reference model for the development of quantum many-body methods [1].
In terms of Fermi surface instabilities, the square lattice has been the dominating focus of theoretical research as it has been hosting the majority of quasi-two dimensional candidate materials for strongly correlated electron systems. More recently, however, the discovery of stronglycorrelated states in moiré materials, i.e. systems based on few-layer stacks of two-dimensional materials such as graphene or transition metal dichalcogenides (TMD) [2][3][4][5][6], has made a strong case for revisiting hexagonal lattice systems of correlated electrons from the viewpoint of stateof-the-art quantum many-body approaches [7][8][9][10][11][12][13][14].
Kagome, honeycomb, and triangular lattices all share the same hexagonal point group symmetry but differ in terms of Wyckoff positions taken by their respective lattice sites. The triangular lattice stands out as the local site symmetry matches that of the hexagonal point group symmetry. It has a high potential to offer exotic many-body states due an intricate interplay between frustration and correlations, see, e.g., [15][16][17][18][19][20][21][22][23] for a recent series of studies on that matter.
In systems such as TMDs, a sizable spin-orbit coupling breaks the spin-rotation invariance. As a consequence, effective models for moiré TMDs often involve several spinsplit bands [24][25][26] which need to be taken into account by adequate quantum many-body approaches. In an attempt to boil down moiré TMDs to its fermiological essence, this can thus lead to an effective model of spin-polarized interacting electrons (or spinless fermions). Note that in the absence of local Hubbard repulsion due to the removed spin degree of freedom, nearest-neighbor density-density interactions are the most elementary terms to consider, which we adopt for our paradigmatic toy model in the following.
A method that has been shown to be quite flexible when it comes to the description of competing instabilities of correlated-electron systems on various lattice geometries and for a broad range of fillings and interactions, is the functional renormalization group (FRG) [27][28][29]. The FRG has been used in numerous studies to identify the leading Fermi-surface instabilities with all competing interaction channels being treated on equal footing [30,31] Within the correlated-electron FRG different schemes have been employed for numerical implementations, most prominently the N -patch scheme, which divides the Brillouin zone into a number of N patches with the representative momenta lying on the Fermi surface [32]. The N -patch scheme allows for a relatively simple and straightforward numerical implementation, but becomes numerically expensive for high momentum resolution and also does not faithfully incorporate momentum conservation. More recently, an alternative scheme -the truncated-unity scheme [33] (TUFRG)based on a decomposition of the different interaction channels [34] has been devised, which separates stronger and weaker momentum dependencies and therefore allows for a more efficient numerical evaluation at high momentum resolution.
In this work, we establish the correlated phase diagram of of spinless electrons on the triangular lattice in the presence of competing interaction channels around Van Hove filling. To that end, we set up both, an N -Patchand a TUFRG approach for correlated fermions without spin-SU(2) invariance. We carefully study the convergence within both schemes and compare them to each other. The motivation of our work is twofold: 1. The FRG represents a very promising scheme for setting up sophisticated numerical implementations that can capture accurate multi-orbital/-band models for moiré TMDs. Our results can then be used for future reference of such implementations.
2. The systematic quantitative comparison between the two FRG schemes provides guidance to the choice of transfer-momentum resolution and form-factor expansions in future TUFRG studies, which are likely to be more appropriate for a faithful description of more involved models due to numerical efficiency.

Model
We consider a tight-binding model for spinless fermions on the triangular lattice where we add a nearest-neighbor density-density interaction, reading Here the operator c ( †) i annihilates (creates) a fermion on lattice site i, such that we allow for nearest-neighbor fermion hopping with rate t. The fermion density operator n i = c † i c i couples to the chemical potential µ to change the filling of the system and V 1 > 0 (< 0) is the strength of the repulsive (attractive) density interaction of neighboring fermions (see Fig. 1). We will study the effects of attractive and repulsive interactions for an extended range of fillings corresponding to µ. The energy band of this model is given via a Fourier transform, yielding with wavevector k = (k x , k y ). We note that at µ/t = 2 the band dispersion features saddle points at the three inequivalent M points of the Brillouin zone (BZ) (Fig. 1), giving rise to a Van Hove singularity (VHS). Our investigations of the emergent many-body instabilities of the system will be carried out in the vicinity of the VHS, but also beyond.

Fermionic functional renormalization group
The fermionic functional renormalization group (FRG) [27,28] has been established as a versatile approach to treat strongly-correlated electrons without bias towards a specific mean-field channel [30,31]. It is rooted in the functional integral description of quantum many-body systems and it allows for the investigation of a broad range of models without specific limitations for their kinetic or interaction parameters. Generally, the FRG acts as functional implementation of the Wilsonian renormalizationgroup (RG) idea, namely, one starts at an ultraviolet (UV) cutoff scale Λ UV and successively takes effects of fermionic fluctuations into account by approaching the infrared (IR) limit Λ IR = 0.
While the FRG description of a selected model is at a formal level exact, one needs to decide for truncations of the description to derive a feasible numerical application from the general principles. In the situation of competing interactions, this truncation will mostly concentrate on the evolution of the two-particle vertex as an indicator for emerging Fermi-surface instabilities. In the past, this has led to many successful applications of the method to strongly-correlated electron systems, for example, for models of spin-rotational invariant electrons on triangular and honeycomb lattices, see, e.g., [8][9][10][11][12][13][35][36][37][38][39]. In addition more specific models of these geometries have been investigated aiming at the description moiré materials [23,[40][41][42][43][44].
The FRG flow is realized by solving a system of coupled differential equations interpolating between the UV and the IR limit. In this work, we want to compare two specific computational schemes to track this FRG evolution of running couplings: (1) the N -patch scheme, which was one of the first well-established methods within the fermionic FRG framework, and (2) the truncated-unit FRG, a more recent approach which goes beyond the patching scheme and allows for a finer grained momentum resolution.

Flow equations
Our starting point is the action for a many-electron system whereψ, ψ are Grassmann-valued fields. Here, the quadratic term includes the free propagator G 0 (ω, k) = 1/(iω −ξ(k)) with Matsubara frequency ω and single-particle dispersion ξ(k), and the bracket (., .) denotes integrations over continuous and summations over discrete indices. The second term S int [ψ, ψ] in Eq. (3) is an interaction term, which can be read off directly from the interaction part of the microscopic Hamiltonian in Eq. (1). With the help of the action S, we can define the Schwinger functional and its Legendre transform -the effective action -Γ [ψ, ψ] = (η, ψ)+(ψ, η)+G[η, η] with ψ = −∂G/∂η andψ = ∂G/∂η, which generates the one-particle irreducible (1PI) correlation functions [45]. The central step for setting up the renormalization group scheme amounts to regularizing the free propagator by an infrared cutoff Λ, such that G 0 (ω, k) → G Λ 0 (ω, k). The cutoff implementation is, in some sense, arbitrary, as long as the ultraviolet (Λ → ∞) and infrared limit (Λ → 0) are smoothly connected. Here, we opt for implementing the temperature flow scheme introduced by Honerkamp and Salmhofer [46] which is employed for both FRG implementations. For now, however, we will keep the discussion general and refer the reader to App. A for details on the T -flow.
Having regularized the bare propagator, the effective action Γ [ψ, ψ] becomes scale dependent and its flow is governed by an exact differential equation [47], which reads where Γ (2)Λ = (∂ψ, ∂ ψ ) T (∂ ψ , ∂ψ)Γ Λ is the matrix of second derivatives of Γ Λ . Here, the appearance of the matrix of second functional derivatives of the effective action Γ (2) necessitates some truncation to derive a closed set of equations for the 1PI vertex functions. We employ a standard approximation scheme, which (1) neglects selfenergy insertions, such that undifferentiated fermion lines correspond to bare, unrenormalized propagators, (2) sets external Matsubara frequency arguments to zero and, simultaneously, does not account for the frequency dependence of the two-particle vertex and (3) truncates the three-particle vertex from the flow equations (an in-depth discussion of these approximations is reviewed in [30]). As a result, we obtain flow equations for the static two-particle vertex V (k 1 , k 2 , k 3 ) (the fourth momentum is fixed by momentum conservation), which allow us to determine Fermi liquid instabilities in an unbiased way. For spinless fermions, the flow equations read [48,49] d dΛ where denotes the pairing or particle-particle channel the crossed particle-hole channel and the direct particle-hole channel, respectively. Here the integral is defined as k = A −1 BZ T BZ dk iω and k = (k, ω) where A BZ is the area of the Brillouin zone.
Integrating these equations starting with the bare coupling in the Λ → ∞ limit, Fermi liquid instabilities are signified by singular contributions to V . We note that V is a function of three momenta and it is therefore costly to compute. For this reason, we rely on further approximations for its momentum dependence, two of which are presented in the following.

N -patch FRG
The first, well-established approximation of the momentum dependence assumes that the two-particle vertex is constant along elongated patches in momentum space [30].
To implement the patching scheme, we define a mapping π : 1.BZ → Z N FS , identifying momenta k in the first Brillouin zone with their nearest-neighbor π(k) in an angular discretization Z N FS of the Fermi surface, which consists of N points, see Fig. 2. This way, irrelevant couplings perpendicular to the Fermi surface are projected out and the vertex is fully determined by its value on the central patch points, which we place equidistantly. Note, that this treatment of the momentum dependence of the vertex spoils momentum conservation, since the fourth momentum k 4 = π(k 1 ) + π(k 2 ) − π(k 3 ) of the projected vertex V (π(k 1 ), π(k 2 ), π(k 3 )) will in general not align with a patch point and therefore require an additional transformation with π.
N -patch FRG calculations were successfully employed to track the flow of marginal couplings for prototypical model systems of high-T c superconductivity such as iron − 2π pnictides and cuprates, see, e.g., Refs. [29][30][31] and references therein. This legitimates the method as a valid starting point to determine the leading instabilities around the Fermi surface fixed point. In summary, the patching scheme describes the vertex with three projected momenta, i.e. V (π(k 1 ), π(k 2 ), π(k 3 )), such that for a selection of N patches, the numerical cost will scale with N 3 . In this work, we implemented a resolution of the Fermi surface using N = 192 patches.

Truncated-unity FRG
The truncated-unity FRG (TUFRG) [33] allows for a high resolution of the full Brillouin zone, i.e. in contrast to the N -patch scheme, it is not restricted to the Fermi surface. Instead, one can chose arbitrary points of momenta to evaluate the flow equations. The derivation of the TUFRG approach is based on the fact that the singular behaviour of instabilities are mainly depending on the transfer momenta inside the loops in Eqs. (6)-(8) connecting the two vertices [34]. Specifically, they are k 1 + k 2 in τ pp , k 1 − k 4 in τ ph,c and k 1 − k 3 in τ ph,d . Consequently, the interaction is reparametrized into different channels such that each object is accounting for one of the transfer momenta. In practice, V Λ is decomposed as where V Λ,0 (k 1 , k 2 , k 3 , k 4 ) accounts for the initial conditions of the model. The channels carry the important transfer momentum as first argument and each channel can be interpreted as representing a specific kind of interaction. The choice of these three channels was initially motivated by models of spinful fermions, where P will represent a pairing interaction, and depending on spin combinations, C and D represent magnetic and density-density interactions. Since our model Eq. (1) is spinless, both channel C and D will eventually represent density-density interactions and this choice is therefore redundant We keep this representation anyway such that a transfer of this method to a spinful model can be done in a transparent way.
To relate the channels to the diagrams with the same important momentum, we define the flow equations where Λ was dropped for brevity. Since the the last two momenta of the channels are deemed as less important, we will expand them in form-factors: with X ∈ {P, C, D}. This expansion can be imposed as long as the form-factors are forming a unity: indicates an instability in the particle-hole channels. Solving the respective gap equation using the renormalized vertex in (b), yields the largest eigenvalue for transfer momentum q = 0. The so-determined order parameter, indicated by light blue dots in (c), has an extended s-wave symmetry and transforms in the A1 irrep. with both first and second neighbor harmonics (a fit to the numerical data is plotted as a dark blue line).
The channel decomposition and the unity of the formfactors can now be used to reformulate the initial flow equations into a form which offers a computational advantage.
In the TUFRG approach we derive flow equations for P l,l (q), C l,l (q), D l,l (q) by taking the derivative d dΛ and inserting form-factor resolved unities on the right hand side of Eqs. (6)-(8) between the vertices and the loops, eventually leading to separating the three objects momentumwise while connecting them in terms of form-factors. The sum of the form factors introduced with the unity Eq. (15) can then be truncated safely to gain a numerical advantage. The final form of the TUFRG flow equations reads for details of the objects see App. B.1. The flow equations now scale with N q × N 2 f , where N q is the number of momenta q which discretize the Brillouin zone and N f is the number of chosen form-factors, see Fig. 3. In practice one has to choose much less form-factors than patches in the patching scheme. Therefore, we gain a numerical advantage over the scaling of the N -patch scheme (∼ N 3 ) and the freedom to choose a larger number of momenta N q in the Brillouin zone. In this work, we use N q = 180 and N f = 19 for comparison with 192 patches in the other approach. To discuss single points in the phase diagram we use N q = 540 and N f = 19. In the convergence checks we go up to N q = 792 and N f = 61. For details about the choice of momenta and form factors, see App. B.2.

Linearized gap equation
To obtain the gap function ∆(k) for the superconducting instabilities encountered during the FRG flow, we utilize standard BCS theory [50], that is, we perform a mean-field decoupling in the superconducting channel and derive a self consistent gap equation for ∆(k). Close to the critical temperature, where the gap is presumably small, the gap equation can be linearized and resembles an eigenvalue equation, which reads The only input required to solve Eq. (19) and determine the leading contributions to the gap function as the eigenvectors with the largest negative eigenvalues, is then given by the pairing potential For the patching approach, we rewrite the right hand side of Eq. (19) as an integral over a small energy shell − c ≤ ξ k ≤ c FS around the Fermi surface, where the most dominant contribution to the momentum sum stems from. The gap equation thus becomes where the integral evaluates to (21) allows to straightforwardly obtain ∆(k) on the Fermi surface within the patching approach.
If we work with the TUFRG approach instead, we can restore the pairing interaction straightforwardly by calculating the pairing interaction from the P channel, which is just given by the form-factor expansion in Eq. (13): Since the divergence has a sharp peak at q = 0, we set the superconducting pairing interaction as: and identify V BCS (k, k ) = Φ P (k, k ). Thereafter, the gap function is obtained by diagonalization of the N q × N q matrix Φ P (k, k ). Note, that while we have focused on pairing instabilities for the sake of brevity, one can generalize the discussion above directly to instabilities in the particle-hole channels by performing the respective mean-field decoupling and deriving a gap equation with an appropriate density instead of a pairing potential.
4 Attractive case V 1 < 0 We first investigate the case of attractive interactions V 1 < 0 at and away from Van Hove filling µ/t = 2. To that end, we apply both, the N -patch and the TUFRG scheme, and work out the qualitative and quantitative differences between these approaches. In order to generate a common starting point we initialize both methods as follows: the RG flow starts at where W = 9t is the bandwidth of the model. The respective flow equations are integrated down to the infrared, which we numerically define by T IR /t = 10 −5 . If one of the channels diverges, signified by its maximum exceeding 3W , the integration is terminated preemptively. As initial value for the vertex, we set V W (π(k 1 ), π(k 2 ), π(k 3 )) = V 1 in the patching scheme, and V W (k 1 , k 2 , k 3 , k 4 ) = V 1 , Φ W,X = 0 in the TUFRG.

Pomeranchuk instability at Van Hove filling
Tracking the evolution of the attractive case under the RG flow both employed approaches eventually detect a divergence of the particle-hole channels between T /t = 1 and T /t = 0.1, see Figs. 4 and 5. Due to crossing symmetry, which relates the direct and crossed particle-hole contributions, the flows of the respective maxima align and we, thus, reduce our discussion to τ ph,d for brevity.
In the patching scheme, we find the most singular eigenvalue to emerge from the linearized gap equation (see App. ?? for further details) with transfer momentum q = k 1 −k 3 = 0, corresponding to a Pomeranchuk instability [51]. The respective order parameter ψ k ψ k (see (c) in Fig. 4) is found to live in the A 1 irreducible representation (irrep.) of C 6v with an extended s-wave form factor including nearest and second-nearest neighbors. The momentum modulation, induced by the second neighbor harmonic, is, however, quite weak to the constant offset presented by the nearest-neighbor A 1 basis function.
In the TUFRG scheme, the instability almost exclusively affects the onsite-component D 1,1 (q) of the direct particle-hole channel, with a pronounced peak at the Γ point (see Fig. 5 (b)) and in agreement with the patching results. The reconstructed order parameter ∆ D (see (c) in Fig. 5) likewise transforms in the A 1 irrep., including a momentum modulation on the Fermi line. Note that, due to this modulation being weak compared to the nearestneighbor A 1 contribution, this is rather difficult to see from the colormap in Fig. 5.

Superconductivity below Van Hove filling
For fillings µ/t < 2, the Fermi surface is deformed and at some point the Pomeranchuck instability is overruled by a superconducting instability. We observe that, depending on the combination of chemical potential and interaction strength, both employed FRG schemes consistently predict two different kinds of superconductivity with q = 0 for an extended range of fillings.
More specifically, moving away from µ/t = 2 towards smaller values, the Pomeranchuck instability will at first be replaced by a region of f -wave superconductivity. Using the FRG data, we can reconstruct a gap function as detailed in Sec. 3.4. Indeed, we find that the gap function belongs to the one-dimensional B 1 irrep. of C 6v (see Figs. 6 and 7). The size of the filling range where the f -wave superconductivity instability occurs grows for decreasing interaction strength |V 1 |. In the case of V 1 = −0.4 this type of superconductivity is even the only one which persists. Notably, in the case of V 1 = −1.0, where the region is the smallest, the N -patch FRG scheme does not detect f -wave at all while TUFRG still resolves a small domain of this instability.  Lowering µ further, p-wave superconductivity becomes the leading instability, which is described by the twodimensional E 1 irrep. of the same point group (see Figs. 8 and 9). On a mean-field level, it is energetically beneficial for the superconducting order to open a full gap in the quasi-particle spectrum, which can be accomplished, for example, by constructing the superconducting gap ∆(k) as a complex superposition of the E 1 lattice harmonics. This leads to a p + ip superconducting state featuring a finite Chern number C = −1 which is thus topologically non-trivial (see App. D). Qualitatively, the two superconducting instabilities we find here are also consistent with the mean-field study presented in Ref. [52]. We note, however, that our FRG study includes additional fluctuations, which induce the Pomeranchuk instability when approaching Van Hove filling.

Phase diagram of the attractive case
In Fig. 10, we have mapped out the phase diagram for various V 1 /t < 0 using both, the N -patch FRG and the TUFRG. Generally, the phase boundaries, the respective ground state instabilities, and the critical scales are in reasonable agreement. Some deviations in the critical temperatures are visible, in particular, in the regions where the superconducting instabilities occur at very low scales. Notably, the transition from Pomeranchuk to the f -wave superconductivity is in good alignment in both methods while the second transition point towards p-wave superconductivity has a larger difference although deep into this particular phase the methods apparently converge.
To establish the reliability of our results, we have further studied the convergence of the TUFRG approach with respect to the momentum-and form-factor resolution (N q , N f ) in more detail, see App. B.2. 5 Repulsive case V 1 > 0 We now consider the repulsive case V 1 /t > 0. Here, we can expect that the occurring instabilities result from an interplay of the perfect nesting at the Van Hove point, whose effect can be mitigated by changing the filling, and Fig. 11. Patching results for V1/t = 1 at Van Hove filling.
(a) Flow of the channel maxima, indicating a simultaneous divergence in both particle-hole channels, consistent with the TUFRG result in Fig. 12. (b) Plot of the direct particle-hole channel right at the critical scale Tc. The corresponding plot for τ ph,c can be obtained via crossing symmetry, i.e. a permutation φ1 and φ2 and a flip of the overall sign.
a divergent susceptibility in the pairing channel, which eventually induces a superconducting instability.

CDW at Van Hove filling
Similar to the attractive case, both methods detect a divergence of the particle-hole channels for µ/t = 2. An analysis of the possible order parameters ψ k+q ψ k (see Figs. 11 and 12), however, reveals that the leading instability occurs for transfer momenta q, which coincide with the nesting vector M . The FRG results thus indicates the instability towards a charge density wave.

5.2p-wave superconductivity below Van Hove filling
As we have discussed for the attractive case, the Fermi surface loses its nesting property below Van Hove filling, and, thus, fluctuations in the particle-hole channels are weaker (but still finite). In contrast to our previous considerations, however, putative superconducting instabilities would now arise from a different mechanism. Since V 1 /t > 0, pairing is not directly encapsulated by the bare vertex and an attractive interaction in τ pp henceforth needs to be generated by inter-channel feedback during the RG flow. Indeed, both methods find an instability of the particleparticle channel for various fillings µ/t < 0 and, remarkably, the flows of the maxima in the different channels plotted in Figs. 13(a) and 14(a) underline the importance of particle-hole fluctuations for the emergence of superconductivity. While the pairing channel is negligible (in TUFRG) or at least smaller than the other contributions (in the patching scheme), the particle-hole channels first sharply increase and then converge to a constant value, which dominates the vertex. In the low temperature regime, however, an abrupt upturn in the τ pp flow can be observed, which ultimately results in a divergence of the RG flow. The respective gap function again transforms in the E 1 representation of C 6v , but requires both nearest-and second-nearest neighbor lattice harmonics, as indicated by an increased number of nodes on the Fermi surface (see Fig. 13(c) or Fig. 14(b)). We dub this instabilityp-wave to set it apart from its counterpart in the attractive case.
Notably, a complex order parameter constructed solely from the second neighbor E 1 basis functions likewise yields C = −1, whereas superpositions of both the first and second neighbor harmonics can generate an enhanced quantum Hall response due to Chern numbers |C| > 1 (see Fig. 19 for more details).

Phase diagram of the repulsive case
In Fig. 15 we finally show results for the phase diagram obtained from the patching scheme and TUFRG for various fillings and repulsive interactions. Interestingly, the temperature scales for thep-wave superconductor measured in the patching scheme are almost one order of magnitude higher than in TUFRG, though the nature of the instability remains the same. Moreover, the sharp drop in T c between the CDW and superconducting regime is absent in the patching results, where only a soft shoulder is indicative of the transition. Close to Van Hove filling on the other hand, the agreement is more reasonable. Since the central patch points coincide with the saddle points in the latter case, this generates the suspicion that the projection to the Fermi surface might be responsible for the observed discrepancy away from perfect nesting.

Discussion
We analyzed competing orders in a model of spinless electrons on the triangular lattice with nearest-neighbor interaction. Our study was motivated by the observation of correlated states in moiré bilayers of transition metal dichalcogenides. These systems are effectively described  Fig. 13. Patching results for V1/t = 1 and µ/t = 1.7. A superconducting instability, driven by strong particle-hole fluctuations, becomes visible as a divergence of the particle-particle channel (see (a) and (b)). The pairing potential, constructed from τpp at the critical scale Tc, has two degenerate gaps (light red and light blue dots in (c)), which can be fit by a linear combination of first and second neighbor lattice harmonics of the two-dimensional E1 representation of C6v (dark red/blue line).
by interacting electrons on a triangular lattice, although equipped with (pseudo)spin and/or orbital degrees of freedom. To distill out the minimal degrees of freedom, we considered the paradigmatic toy model of spinless electrons and showed that it still possesses a rich interplay of ordering tendencies in the vicinity of a Van Hove singularity. To resolve this interplay, we calculated the effective two-particle interaction vertex in an unbiased way with the functional renormalization group. It is crucial to accurately resolve the momentum dependence of the vertex and we used two different parameterizations -a patching scheme for the Fermi surface and a channel decomposition for the momentum transfers. Both of them give qualitatively consistent results.
With an attractive bare interaction, we find a Pomeranchuk instability in the s-wave channel directly around Van Hove filling and f -and p-wave pairing instabilities in its vicinity for smaller fillings. Within RPA, both the charge and pairing channel can develop an instability, although at weak coupling the pairing channel has a stronger divergence (logarithmic vs double logarithmic). Interestingly, in our calculations, the Pomeranchuk instability in the charge channel develops first due to non-universal effects (beyond the logarithmic scaling). The s-wave Pomeranchuk instability corresponds to a singular compressibility but is not associated with any symmetry-breaking order. This can signal the tendency to phase separation with domains of different density. Another possibility is that the divergence is cured by terms outside of our truncation, e.g., by self-energy terms, and makes room for a subleading instability. The p-wave pairing solution is two-fold degenerate and can form chiral p + ip superconductivity in the ground state. This topological triplet superconducting state breaks time-reversal symmetry and can host Majorana modes on its boundaries.
In the case of a repulsive bare interaction, we obtain a CDW instability closest to Van Hove filling, whose fluctuations mediate unconventional p-wave pairing at smaller fillings. The wave vectors of the CDW are the three nonequivalent M points of the Brillouin zone and the exact charge pattern of the associated order depends on their combination in the ground state. Due to the bare nearest-neighbor repulsion, we find the unconventional pwave pairing to be of extended size described by nearestand next-nearest-neighbor harmonics. This can yield topological p + ip states with higher Chern numbers, which increases, e.g., the number of chiral edge modes and the quantum Hall response.
where the integral includes the Brillouin zone area: k = A −1 BZ dk. These expressions can also be simplified by plugging in the plane wave form-factors exp(ikR l ) (see App. B.2) and expressing V Λ by the decomposition Eq. (9). Therefore the double integral over the Brillouin zone is exchanged by a simple sum over the selected form-factors L : The objects V X,0 l,l (q) encode the initial interaction of the model Eq. (1) by projecting it into the respective channels, see App. C.X l,l represents the Fourier-transformed channels, for example for the pairing channel P :

B.2 Choice of momenta and form-factors and convergence
One has the freedom to select different sets of form-factors as long as the unity condition Eqs. (14)-(15) are fulfilled. The simplest choice of form-factors have the form of plane waves: f l (k) = exp (ikR l ) where R l is a real space vector of the lattice of the investigated model, i.e. in our case the triangular lattice. This choice has the advantage, that the truncation of form-factors can be done within an interpretable reasoning: the inclusion of a form-factor f l (k) will correspond to taking effects of fermionic bilinears with distance R l into account [53]. Since we assume that the emerging physics in the RG flow will be predominately influenced by short-range effects, we will truncate all form-factors which exceed a chosen distance. In our calculations we mostly select N f = 19 form-factors, corresponding to on-site (i.e. R 1 = 0), first-, second-and third-nearest neighbors effects. For the convergence checks in Figs. 17,18 we will also use N f = 37 (i.e. up to 5th nearest-neighbors effects) and N f = 61 (i.e. up to 8th nearest-neighbors effects), see Fig.16. This specific choice of amount of form-factors is based on keeping a hexagonal-shell N s into account. This means, that we will include all plane waves with R l which are on or inside the N s − th hexagon of the real space lattice, cf. Fig.16.
Therefore the numbers N f = 19, 37, 61 correspond to the hexagon-shells N s = 2, 3, 4. For the momentum resolution, we choose evenly placed points in the Brillouin zone. Most of our calculations are done with N q = 180 momenta to compare it with the 192 patching points of the other approach, while for the convergence checks in Figs. 17,18 we also choose N q = 336, N q = 540 and N q = 792.
Actually, one does not have to calculate the RG flow for all momenta N q , but only for a fraction 1/12 × N q . The rest of the contributions can then be restored by symmetry relations since the symmetries of the initial model Eq. (1) are inherited by the flow equations, see [54] for details.

C Initial conditions
The initial condition for the FRG flow is given by the bare two-particle vertex V 0 , which can be directly read off the microscopic model in Eq. (1). For this purpose, one needs to identify the action S int with the vertex at the UV scale, i.e S int = V ΛUV = V 0 , and additionally account for crossing symmetries, such as V (k 1 , k 2 , k 3 , k 4 ) = −V (k 2 , k 1 , k 3 , k 4 ). The initial condition needs to be properly (anti-) symmetrized henceforth. On the level of the Hamiltonian, crossing symmetry can already be made explicit by reordering the Fock space operators as Transforming to momentum space, the initial condition for the FRG flow is thus where we sum over the nearest-neighbor displacement vectors δ. Projecting all momenta to the Fermi surface via π : 1.BZ → Z N FS , Eq. (37) directly serves as the initial condition for the patching scheme.  Fig. 18. Convergence of critical RG scales from TUFRG for the repulsive case V1/t > 0.(a) For the investigation of convergence in increasing momentum resolution Nq we find qualitatively the same phase diagram which quantitative deviations diminish for higher resolution. (b) The the investigation of including more form factors N f we still find qualitative alignment, while the critical temperature slightly grows for including more shells. Since we are primarily interested in the qualitative behaviour and the deviations are not too strong, we use Ns = 2 for the calculations in the sections 4 and 5.

D Chern numbers
To access possible topological properties of pairing instabilities, we consult a Skyrmion winding number formula [55,56] where . | . is the Euclidean scalar product. Here, m(k) denotes the pseudospin vector or Skyrmion magnetization, which follows the winding of the superconducting gap around the Fermi surface. In algebraic form, m(k) is given by where E(k) = |∆(k)| 2 + ξ(k) 2 is the Bogoliubov quasi-particle spectrum. It is immediately clear, that any real or purely imaginary gap function will result in a topologically trivial state with C = 0. In contrast, for a gap function corresponding to a two-dimensional irreducible representation, such as the p-wave instabilities we found in the main text, the possibility of non-trivial topology arises. In principle, one would need to minimize the mean-field free energy for a linear superposition of the respective lattice harmonics and determine whether or not a complex gap function prevails. Here, we resign from employing this variational approach and instead use a heuristic argument. Consider the ground state energy E 0 = − |∆(k)| FS for a gap function ∆(k) which we suppose to live in the complex two-dimensional space corresponding to a doubly degenerate eigenvalue of the linearized gap equation. If this linear combination is either real or purely imaginary, there will be momenta on the Fermi surface where |∆(k)| is gapless and no contribution to the ground state energy is obtained henceforth. If one assumes a complex linear combination instead, |∆(k)| will be fully gapped at the Fermi level and thus, a lower ground state energy is obtained. It is therefore natural to assume, that the energetically more beneficial superposition  Fig. 19. Example calculations for Chern numbers in the E1 representation for µ/t = 1.7. Motivated by our finding of a higher-harmonicp-wave instability for repulsive interactions V1/t > 0 (see Figs. 13 & 14), we perform exemplary computations of C for superconducting gaps of the form ∆(k) = [cos(α)δ E 1 1 (k) + sin(α)δ E 1 1 (k)] + i × [cos(α)δ E 1 2 (k) + sin(α)δ E 1 2 (k)], where δ E 1

1(2)
denotes the nearest-neighbor lattice harmonics of the E1 irrep. andδ E 1 1(2) the respective second neighbor functions. The model is chosen such that we recover the pure first (second) neighbor limit for α = 0 (π).
of lattice harmonics is a complex one. Computing C from the ansatz ∆(k) = δ E1 1 (k) + iδ E1 2 (k) for the nearest-neighbor or second neighbor lattice harmonics δ E1 of the E 1 irrep., for example, we find C = −1 over the entire range of fillings where the p-wave instability occurs. An admixture of both, the first and second neighbor functions may, however, yield a strongly enhanced Chern number, as exemplified in Fig. 19.