{'fraction_non_alphanumeric': 0.03259387327716359, 'fraction_numerical': 0.009082939658111333, 'mean_word_length': 4.463150192896102, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multiorbital effects in the functional renormalization group: A weak-coupling study of the Emery model

We perform an instability analysis of the Emery three-band model at hole doping and weak coupling within a channel-decomposed functional renormalization group flow proposed in Phys. Rev. B 79, 195125 (2009). In our approach, momentum dependences are taken into account with improved precision compared to previous studies of related models. Around a generic parameter set, we find a strong competition of antiferromagnetic and d-wave Cooper instabilities with a smooth behavior under a variation of doping and additional hopping parameters. For increasingly incommensurate ordering tendencies in the magnetic channel, the d-wave pairing gap is deformed at its maxima. Comparing our results for the Emery model to those obtained for the two-dimensional one-band Hubbard model with effective parameters, we find that, despite considerable qualitative agreement, multi-orbital effects have a significant impact on a quantitative level.


I. INTRODUCTION
Unconventional superconductivity is an extensively discussed topic of condensed matter physics. For cuprate high-T c materials, resonating valence-bond 1 and spinfluctuation 2 mechanisms have been proposed besides other approaches that require low-energy bosons in addition to an electron-electron interaction. While the former relies on a strong-coupling scenario, the latter may also apply to iron-based and other types of unconventional superconductors where the electrons interact more weakly.
A particular problem in the theoretical description of spin-fluctuation induced superconductivity is that often there is no clear separation of energy scales, i.e. the spinfluctuations are built up at least in parts by the same electronic degrees of freedom as the pairing tendencies. At weak coupling, renormalization group (RG) methods can be used to study both types of electronic fluctuations on equal footing. Both, for the two-dimensional one-band Hubbard model at weak coupling [3][4][5] and for iron-based superconductors, 6 the interplay of the spin-density wave (SDW) and superconducting (SC) channels has been investigated within RG approaches involving a small number of running couplings. Such studies have helped to understand the interacting ground states qualitatively. More generally, the RG interpolates between microscopic models at high scales and and effective models at lower scales. As they take into account the dominant fluctuations on equal footing, RG methods may provide a manybody framework that can become quantitatively precise at least for some unconventional superconductors whose parameters fall into the weak-coupling sector.
In the attempt to become more quantitative, functional renormalization group (fRG) techniques for fermions (for a recent review, see Ref. 7) shall be useful, as they allow for studying the flow of momentumdependent coupling functions around the Fermi surfaces which can resolve more details of the models under investigation. Using fRG, the ordering tendencies of the oneband Hubbard model at weak coupling have been classified in a number of works (see, e.g. Refs. [8][9][10]. Among others, also models for iron-based superconductors have been investigated within the fRG framework. [11][12][13] Here, relevant variations of the superconducting gap function due to material-dependent electrons structure differences have been found. 14 Generically, the multiple Fermi surfaces (hole and electron pockets) of iron-based superconductors with their varying orbital character have been proven to add new aspects such as competing pairing channels and gap anisotropies. 14,15 Triggered by these insights, multiband effects are now explored and revealed in a number of other systems (e.g. Refs. [16][17][18]. It therefore appears rewarding to study multiband effects in models for the cuprates. The main question is whether the underlying multiband character gives rise to deviations from the one-band picture, even though the single band at the Fermi level has overwhelming d x 2 −y 2 -character. Furthermore, some proposed ordering phenomena in Copper oxide planes such as ring currents 19,20 are not describable by models involving one d-orbital on the Copper atoms only. Hence, also from this perspective, an fRG study of the Emery model, 21,22 which includes the oxygen porbitals, appears to be worthwhile. Of course, within an fRG framework, cuprate models have to be considered at (possibly unrealistic) weak coupling strengths. However, as for the one-band Hubbard model, similarities between weak and strong coupling be-havior can be expected. Moreover, multi-orbital cuprate models seem very well suited as a testbed for methodological developments. For example, the impact of an interaction rendered nonlocal by the transformation to the band language 23 can be more easily studied than in the context of pnictide models. A recent fRG study of two-and three-orbital models for the cuprates with orbitals only residing on the Copper atoms has revealed that these so-called orbital-makeup effects may have a significant impact on the phase diagram. 16 In the strong-coupling case, other powerful many-body methods are applicable, such as dynamical mean-field theory [24][25][26][27][28][29][30][31] (DMFT), dynamical cluster quantum Monte Carlo (DCQMC) techniques 32,33 and the so-called variational cluster approach [34][35][36][37]. Also these methods are still being extended and developed further in order to hopefully promote a deeper understanding of the cuprate superconductors. Unfortunately, a practicable strong-coupling truncation of the fRG flow equations for fermions, which would allow for direct comparison between fRG and those other methods, is not known so far. As the fRG approach pursued in this work is hence confined to weak coupling, we do not primarily seek to compare the results of these strong-coupling approaches to our own findings for weak coupling.
In addition to a weak-coupling truncation which neglects three-particle and higher interaction terms, the fRG studies mentioned so far share a common feature: While the dependence of the interaction on the Matsubara frequencies of the fermionic fields is omitted by projecting to zero frequency, their momentum dependence is projected to the Fermi surface, which is divided into n patches. In order to keep the numerics tractable at increasing resolution, a channel decomposition has been proposed for the frequency dependence 38 and for the momentum dependence 39 of the coupling functions. In such an approach, the coupling function of the two-particle interaction is split into at least three contributions which depend strongly only on one frequency/momentum variable and weakly on the other two ones. For the one-band Hubbard model in two dimensions, the projection to zero frequency has been relaxed in a channel-decomposed approach. 40,41 In all those works, the weak momentum dependences have been accounted for by simple ansätze within an exchange parametrization. These ansätze have been found to describe the weak momentum dependences well in large parts of the phase diagram with the exception of a region around the transition from d-wave SC to ferromagnetism (FM). 42 Let us note in passing that also graphene and the kagome lattice models have recently been studied within a related channel-decomposed fRG treatment, dubbed singular-mode-fRG (SMFRG). 43,44 As for fRG studies of models for Copper oxide planes, the phase diagrams in Ref. 16 apparently show some discretization artifacts related to the projection onto the Fermi surface. Hence, employing a channel-decomposed fRG approach for multi-orbital models appears natural: Since the orbital makeup induces some non-locality in the interaction of multiband Hubbard models, it would be advantageous to resolve its momentum dependence away from the Fermi surface. Moreover, the material characteristics contained in such a multi-orbital model should appear in the high-rather than in the low-energy sector of these models, which makes some momentum resolution perpendicular to the Fermi surface even more desirable. Furthermore, the ansätze for the weak momentum dependences underlying the exchange-parametrization of Ref. 39 seem questionable in the presence of orbital makeup. Therefore, instead of restricting the study to a few form factors, we will use the channel-decomposed flow with momentum dependences discretized on chessboards in the present work. Note that this can be viewed as wavevector-based band-picture analog of the representation of the channels in terms of fermion bilinears on (short) bonds of the lattice in the real-space/orbital picture that is at the heart of the SMFRG approach by Q.-H. Wang et al. 43,44 . As we are primarily interested in a comparison to previous work in the band picture, the wavevector-based approach seems more adequate, but in principle the two setups can be transformed into each other. The SMFRG approach usually deals with all orbitals or bands of the effective model together, while our current study is reduced to the conduction band at the Fermi surface.
This work is dedicated to a channel-decomposed fRG instability analysis of the Emery model and to the comparison to the one-band Hubbard model with effective parameters. It is organized as follows: In Sec. II, we introduce the Emery model and the corresponding effective action for its conduction band. Then we give a prescription how the parameters of an effective one-band Hubbard model can be calculated and classify different kinds of multi-orbital effects. As we show in Appendix A, our coupling functions of the Emery model show the same trivial point-group behavior as those of the one-band Hubbard model and therefore this prescription is viable. Sec. III is devoted to the fRG formalism. After introducing the general form of the RG flow equations, we elaborate on a channel decomposition of those flow equations in Sec. III B. In Sec. III C, we further comment on the exchange parametrization of Refs. 39-41 and on how approximations made in these studies can be understood in a group-theoretic sense.
Our numerical results for the Emery model are discussed in Sec. IV. First, we comment on the strong competition of the antiferromagnetism (AFM) and dSC instabilities observed for most parameter sets considered in this work in Sec. IV A. The subsequent discussion of the dependence of the stopping scale on doping and oxygenoxygen hopping suggests that the system exhibits a firstorder transition between these to phases. In Sec. IV D, we compare the phase diagrams of the Emery model and of the one-band Hubbard model with the corresponding parameters. We finally conclude in Sec. V by discussing the importance of different kinds of orbital makeup effects.

II. MODEL
A. Three-orbital Emery model In this work, we study a three-orbital model introduced by Emery 21,22 for the description of the Cu-O planes of the high-T c compounds. Its Hamiltonian reads as where k denotes the 2D-momentum, µ the chemical potential and σ the spin index and where the interaction term H int enters in normally ordered form. The annihilation operators d σ (k) for the Cu 3d-orbitals and p x,σ (k), p y,σ (k) for the oxygen 2p-orbitals form the components of the orbital pseudo-spinor Here and throughout we set the lattice constant (i.e. the spacing between neighboring Copper atoms) to unity. The one-particle part of H is then determined by the matrix s x,y = sin (k x,y /2) , c x,y = cos (k x,y ) , with on-site energies ǫ d , ǫ p and hopping integrals t pp , t pd . We use LDA-values of these parameters 45 for La 2 CuO 4 as a starting point of our analysis. The importance of the oxygen-oxygen hopping t pp has been extensively discussed 22,26,[28][29][30] for the strong coupling case. At weak coupling, in its absence the leading instability of the system would correspond to commensurate AFM due to perfect nesting. In the basis chosen in Eq. (1), H 0 is not 2π-periodic due to the s x/y -entries. As pointed out in Ref. 46, all momenta still must be folded back to the first Brillouin zone.
In addition to this one-particle Hamiltonian, we consider a short-ranged interaction term where the brackets ij and jj ′ indicate that the sum only runs over neighboring orbitals of the respective types. We will restrict our study to weak interaction, i.e. the typical energy scales of the interaction are about one order of magnitude below the values given in Ref. 45 for La 2 CuO 4 . Interaction terms involving the oxygen porbitals are weak compared the the dominating U d -term Within the fermionic fRG approach pursued in this work, it turns out to be advantageous to write the quadratic part of the Hamiltonian in its diagonalized form, i.e. in terms of bands instead of orbitals. The field operators then correspond to Bloch states that do not get mixed by the one-particle part of the Hamiltonian. The band dispersion of the Emery model is depicted in Fig. 1 for typical parameter values. The chemical potential µ is chosen to values around van-Hove filling µ vH where the Fermi surface touches the saddle points at (0, π) and (π, 0) of the uppermost band. We then obtain one conduction band which is separated from two valence bands by an energy gap of about four times its width. Through the orbital weight imposed by the unitary transformation from the orbital to the band picture, the interaction acquires a nontrivial momentum dependence, dubbed orbital makeup 23 by some authors. From the form of the one-particle Hamiltonian Eq. (1), one finds that the hybridization of the d-and p-orbitals grows from the center to the boundary of the Brillouin zone.

B. Effective one-band action
In principle, an appropriate low-energy solver could be directly applied to the full three-band model. Such a solver would effectively resum diagrams up to infinite order, even if the underlying concept is of non-perturbative nature. The result of such a resummation at weak to moderate coupling will be dominated by diagrams with internal legs on the conduction band. It should therefore suffice to treat the valence bands perturbatively only up to a certain order.
In the fRG literature on multi-orbital models (see for example Ref. 47) high-energy modes above some ultraviolet cutoff or in bands that do not cross the Fermi level are usually neglected. In a recent publication, 48 we have discussed the impact of the most relevant perturbative corrections from the modes above this cutoff. In the present work, those modes correspond to the valence bands. As a starting point, we have considered an effective action for the conduction band represented by the Grassmann field χ − with the inverse propagator D − as a starting point. Here, the functional integral defining the effective interaction V runs over the Grassmann field χ + representing the valence bands with inverse propagator D + and S (4) denotes the bare interaction of the original model. Formally, all connected Greens functions of the conduction band can be reproduced exactly from S eff . In practice, however, the functional integral in Eq. (2) is evaluated perturbatively. In Ref. 48, we have argued that the most relevant corrections to the bare interaction of the conduction electrons (or holes) are given by a second-order tree-diagram with six external legs on the conduction band. This diagram is depicted in the second line of Fig. 2. In the following, we will also study the impact of the three-particle interaction on the fRG flow. If this three-particle term is neglected, the (t, t ′ ) oneband Hubbard model seems to be a good candidate for an effective low-energy model. Such a description can be considered valid if a low-energy solver yields similar results for S eff and for the action corresponding to the one-band Hamiltonian where n σ,i = c † σ,i c σ,i and where ij indicates that the sum runs over second neighbors.
Before we give a prescription how to calculate the effective parameters t, t ′ and U eff , let us elaborate on the FIG. 2. Two-and three-particle vertex of S eff for the effective one-band model. Small filled vertices correspond to the bare interaction S (4) , whereas the vertices of S eff on the left hand sides are denoted by empty circles. Solid lines represent the conduction band and dashed lines the valence bands. There are no propagators attached to the external legs. Self-energy effects will be neglected.
differences between the effective action Eq. (2) for the conduction band and the effective one-band Hamiltonian Eq. (3). The frequency dependence of the action corresponding to Eq. (3) will of course be trivial in contrast to the frequency dependence in S eff which is generated by integrating out the valence bands. Throughout this work, we will, however, project to zero frequency and hence we will not discuss such effects. Recently, frequency-dependent RG flows have been analyzed for the two-dimensional one-band Hubbard model. 40,41 In this context, also frequency-dependent multi-orbital effects appear to be worth further study.
As far as the momentum dependence is concerned, the Hamiltonian Eq. (3) is short-ranged in the sense that all terms are either on-site or describe hopping between first and second neighbors at most. In contrast, S eff may contain quite long-ranged terms both in the dispersion and in the two-particle interaction. At van-Hove filling, the restriction to only two hopping integrals in Eq. (3) can be justified in the spirit of a gradient expansion. Since we have a diverging density of states at the saddle points A = (0, π) and B = (π, 0), the integrals over momenta of internal lines in a perturbation expansion will be dominated by a small region around these saddle points. Consequently, in leading order in an expansion around the van-Hove points, the exact and an approximate dispersion should coincide. Since a k x k y -term is forbidden by symmetries, only two expansion coefficients remain in leading (second) order. They can be expressed in terms of effective first and second neighbor hoppings t and t ′ .
Away from van-Hove filling, the situation may be more involved and hopping terms between non-neighboring unit cells of the direct lattice may be needed for an effective model. Since this corresponds to Wannier functions that have support on more than one unit cell, such a description is not really of tight-binding type. The effective two-particle interaction will also bear traces of the multi-orbital character of the underlying model. More precisely, the orbital makeup renders the interaction nonlocal. Whether this feature plays a role for the low-energy physics remains a question to be answered by applying a low-energy solver.
Before doing so, we should however give a prescription how to obtain the effective Hubbard parameters t, t ′ and U eff . From the comments made above, a gradient expansion around the saddle points of the dispersion appears natural as a guiding principle. So the calculation of t and t ′ is straightforward.
For the interaction, we proceed as follows. First, let U (k 1 , k 2 , k 3 ) denote the coupling function that appears in the parametrization of the effective two-particle interaction which exploits the U(1), SU(2) and translation symmetries. 49 We further assume a basis of Bloch states that ensures the C 4v -symmetry of this coupling function.
In Appendix A, we will show that such a Bloch basis exists. This property is nontrivial, since the oxygen porbitals are mapped onto one another upon a rotation of π/2. As a second step, we also expand the two-particle interaction around the saddle-points. In leading order, we then have four two-particle couplings, namely These are the couplings of the so-called two-patch approximation. 5 For the Hubbard model, all four twopatch couplings g i are equal to U eff . For given V (4) , we therefore take the average of the two-patch couplings U eff = i g i /4 as the effective Hubbard interaction, while the hopping parameters t and t ′ are calculated from a gradient expansion. This means that we choose the parameters of the one-band Hubbard model such that it has a common two-patch approximation with S eff and that we further restrict the hopping terms to neighboring unit cells and the interaction to an on-site density-density term in the effective one-band Hamiltonian.
Note that, in contrast to the famous work by Zhang and Rice, 50 which has been tailored rather for the strongcoupling case, this method is non-perturbative in the hybridization between the orbitals.

C. Classification of multi-orbital effects
Before we discuss our method, a classification of multiorbital effects seems to be in order. Clearly, if we have a full model [in the present case S eff given in Eq. (2)] and an effective one-orbital Hamiltonian [e.g. the Hubbard Hamiltonian in Eq. (3)], effects contained in the full, but not in the effective model are of multi-orbital character. Of course, the multi-orbital nature cannot be attributed to a certain effect without referring to a prescription according to which the full model is mapped to its effective single-orbital counterpart. So, in the following, multiorbital effects will be classified with respect to the above prescription based on a gradient expansion around the saddle points of the dispersion.
In this picture, multi-orbital effects decay into three classes, namely i) effects resulting from the three-particle interaction, and other higher-order vertices generated by the high-energy bands, ii) orbital-makeup effects, in particular a detuning of the two-patch couplings g i , iii) hopping between non-neighboring unit cells.
The three-particle and higher-order vertices responsible for the effects listed as i) appear as (perturbative) corrections to the bare interaction in the effective interaction V in Eq. (2). Since these corrections also contain internal loops with valence-band propagators, they are in general frequency-dependent. Throughout this work, we will however neglect four-particle and higher-order terms as well as contributions with closed valence-band loops (see Fig. 2 for the remaining diagrams). Since we then include only tree-diagrams of bare vertices in the perturbative expansion of V (i.e. diagrams that are reducible to bare vertices by cutting one internal line), the twoand three particle interaction terms of S eff are frequencyindependent in this approximation. Hence, the frequency dependence of S eff is completely neglected in this work.
In the present work, the three-particle term of S eff is either neglected or fed back into the flow equation of the two-particle vertex using an extended truncation of the fRG flow equations. 48 It will turn out to play a minor role due to the large gap between the conduction band and the valence bands.
As far as orbital makeup effects are concerned, we thus restrict the discussion to the two-particle interaction although the three-particle term obviously bears signatures of orbital makeup. Deviations of the two-particle interaction from the on-site Hubbard type manifest themselves in a detuning of the two-patch couplings g i and in a nontrivial momentum dependence also away from the saddle points of the dispersion. These effects have been listed as ii) above. This implies that orbital makeup effects can be partly understood with the help of the flow equations in the two-patch approximation. 5 By looking at the bare one-particle part Eq. (1) of the Emery Hamiltonian, however, we observe that the hybridization between the d and p orbitals is stronger at the saddle points than in other parts of the Brillouin zone. For example, at the origin in k-space, there is no hybridization at all. Therefore, we have a pronounced momentum dependence which may lead to effects that cannot be captured in the two-patch approximation.
The effective action S eff and the effective single-orbital Hamiltonian H eff also differ in their quadratic parts. If the dispersion of the conduction band in the former is expanded around its saddle points, also hopping terms between non-neighboring unit cells appear in the coefficients in subleading orders. As already mentioned, these longer-range hopping integrals do not fit well into a tightbinding picture, as the corresponding Wannier functions would have long tails. Since the conduction band is predominantly of d-orbital character, the hybridization with oxygen p-orbitals can be said to create such tails. We expect hopping between non-neighboring unit cells to play a minor role for the fRG results at weak coupling.

A. General framework
In order to now extract physical properties at low temperatures for the effective one-band action, we have to resort to low-energy solvers. In the fRG approach used in this work, an exact flow equation interpolates between a microscopic action and the generating functional of the one-particle irreducible (1PI) vertices. More precisely, we are dealing with a first-order differential equation describing the renormalizations of this functional that occur when some infrared cutoff λ is lowered. The exactness of the flow equation allows for non-perturbative truncations, 7,51 where fermionic two-particle interaction is decoupled by a Hubbard-Stratonovich transform. At weak to moderate coupling, one may as well resort to a vertex expansion 7,49 in order to avoid some biases of the partially bosonized approach. The purely fermionic vertex-expansion approach is perturbative in the effective interaction (and hence of weak-coupling nature) and resums diagrams up to infinite order in the bare interaction in this way. If the two-particle interaction is parametrized in the same way as in Eq. (4), but with a renormalized coupling function V (k 1 , k 2 , k 3 ) instead of U (k 1 , k 2 , k 3 ), the righthand side of the resulting flow equation for the two-particle vertex comprises five diagrams. T pp denotes the particle-particle diagram and T cr ph the crossed particle-hole diagram while T d ph comprises three direct particle-hole diagrams including vertex corrections and screening. For their precise form, we refer to Appendix B 1. Three-particle-feedback corrections to these flow equations are given in Appendix B 2.
In this article, we employ the Ω-scheme 39 regularization, where the propagator G(k) is multiplied by a smooth frequency regulator according to This regularization scheme does not suppress a Stoner instability as a momentum-shell cutoff would.

B. Channel decomposition
In many previous works, frequencies have been projected to zero and the flow of the self-energy has been neglected. More recent studies 40,41,52 on one-band models taking into account parts of the self-energy and frequency-dependent vertices show that the flows to strong coupling are not changed in character if these two approximations are made simultaneously. They considerable facilitate the loop integrations in the flow of the two-particle interaction and reduce the number of running couplings. Since our focus rather lies on orbital makeup effects, we will hence resort to these approximations in the following.
Although this simplifies the flow equations to be solved, a direct and unbiased discretization of all external momenta in the remaining flow equation for the coupling function V (k 1 , k 2 , k 3 ) is still too costly from a numerical viewpoint. In an older approach to handle this problem, the momentum dependence of V was projected to a finite number of patches on the Fermi surface. 9 This Fermi-surface patching was designed to reproduce the low-energy physics properly, but renormalizations away from the Fermi surface are only crudely approximated. In multi-orbital problems, already the bare interactions, expressed in band representation, show a significant wavevector-dependence away from the Fermi surface. As a consequence, the Fermi-surface patching in Ref. 16 is plagued by discretization artifacts. Hence, an approach that is more suited to capture orbital-makeup effects is desirable. (Nevertheless, meaningful results for multiband models can by obtained from Fermi-surface patching fRG as in Refs. 13-18, 47, and 53.) As in Refs. [39][40][41]43, and 44, we therefore decompose the flow into three channels with the bare interaction U and the pairing, spin and charge coupling functions Φ SC , Φ M and Φ K . These coupling functions are generated during the flow according . The precise form of the one-loop terms on the right-hand side of these equations is given in Appendix B 1. The first argument of the Φs corresponds to the total or transfer momentum of the loops in the particle-particle and particle-hole channels, respectively. For weak coupling, the coupling functions therefore should depend strongly on their first argument (bosonic momentum) and only weakly on the other two (fermionic) momenta. The dependence on the fermionic momenta can be accounted for by an exchange parametrization based on a form-factor expansion (FFE). In Refs. 39-41, only s-and d x 2 −y 2 -wave form factors have been taken into account. It has however been shown 42 that the remainder terms neglected in that work can have a substantial impact on the stopping scale in some parts of the phase diagram of the one-band 2D Hubbard model.
In the present article, we therefore pursue a different approach. Instead of employing a FFE, we directly patch all three momenta of the Φs while we project to zero frequency. The fermionic momenta then live on a much coarser grid than the bosonic ones. The finest resolution is only used for the bosonic momenta around potential divergencies and for the internal loop momenta close to the Fermi surface. If not indicated otherwise, we use 6×6 fermionic quadratic patches and a bosonic resolution of 18 × 18 and of 126 × 126 patches away from and close to possible ordering vectors, respectively. The Pauli principle, point-group and particle-hole symmetries reduce the number of independent couplings further. Note that, close to the transition from d x 2 −y 2 -wave SC to ferromagnetism in the one-band Hubbard model, this resolution would be still too coarse in the fermionic momenta. 42 For the parameters considered in this work, however, we are far away from such a transition and a deformation of the form factors at low energies should at least qualitatively be captured within our approach.
Since we project to zero frequency, the Matsubara sums over internal frequencies can be performed analytically as in Ref. 39. For the Ω-scheme regularization employed in this work, the remaining two-dimensional loop integrals are numerically challenging. They can be efficiently calculated using an adaptive routine. 54 For different patches, these loop integrals are evaluated in parallel using OpenMP and then stored. When we subsequently assemble the diagrams, we again OpenMP parallelize the calculation for different external momenta.
We have sent external fields breaking the U(1), SU(2) and/or space group symmetries to zero right from the beginning. The resulting decrease of computational effort, however, comes with a price. The limit of vanishing external fields is only physical after the thermodynamic limit has been performed. Consequently, spontaneous breaking of a symmetry manifests itself in a flow to strong coupling. We therefore stop the flow when the maximum of the coupling functions reaches 7.7 t pd , and interpret the stopping scale λ c as an estimate for the critical scale of the respective instability. 55 Note that this stopping condition is rather weak, since 7.7 t pd corresponds to about thirteen times the bandwidth of the conduction band. In comparison, in Ref. 39, the flow is stopped once the interactions exceed 2.5 times the bandwidth. For the present work, we have not chosen a stronger stopping condition in order to keep track of the strong competition of AFM and dSC instabilities discussed in Sec. IV A.

C. Remarks on form-factor expansions
Now we comment on the form-factor expansion (FFE) put forward in Ref. 39. In particular, we shall analyze to what extend a channel-decomposed, renormalized twofermion interaction can be conveniently expressed as one resulting from a small number of bosonic channels. Regarding the classification of such order-parameter fields, we proceed similarly to Vojta et al. in Ref. 56, where different types of commensurate ordering within the d x 2 −y 2wave superconducting phase have been classified according to the irreducible representations (IRs) of the point group. The underlying group-theoretical lemmata will be laid out in Appendix C. For the SMFRG approach, similar considerations have been undertaken. 57 In the present case, the coupling functions Φ SC , Φ M and Φ K may be decomposed in the spirit of a Hubbard-Stratonovich transform using a set of orthonormal form factors f i . In the Cooper channel, for example, we have with bosonic propagators D ij (l). The form factors play the role of fermion-boson vertices, with indices i and j labeling different bosonic species (flavors). In the present from of the exchange parametrization, they are frequency-independent and can therefore be chosen real. They obey the orthonormality relation dq f i (q) f j (q) = δ i,j . Here and throughout, integrals dq run over the whole Brillouin zone and a normalization factor has been absorbed into the measure such that dq 1 = 1. For given Φ SC , the matrix elements D ij (l) are thus uniquely s-wave p-wave d x 2 −y 2 -wave dxy-wave g-wave defined. Since the two-particle coupling functions are C 4v -symmetric, it appears natural to choose basis functions of the IRs of C 4v as form factors, for example (Note that the IR corresponding to a p-wave is twodimensional, while the other ones are one-dimensional.) When one projects to zero frequency, the coupling function can be fully recovered by using a complete set of form factors. By (anti)symmetrizing the real Fourier basis functions on the first Brillouin zone with respect to the C 4v point-group operations, one can easily construct a complete basis set with elements that transform according to the IRs (cf. Tab. II). By equivalence transformations of the IR, these form factors can be rendered well-behaved under C 4v in the sense of Appendix C. The form factors mentioned above are the most slowly varying basis functions of the respective IRs, which corresponds to the formation of exchange bosons from constituents residing on the same site or on neighboring unit cells.
For the bosonic ordering vectors l = (0, 0) and (π, π), the little group L l equals the full point group C 4v . According to Corollary C.1, which is proven in Appendix C, matrix elements of D(l) mixing bosons of inequivalent IRs vanish at these momenta. Different form factors transforming according to equivalent IRs may however mix. In the following, we shall refer to this effect as to the admixture of higher harmonics. In flavor space, the nonvanishing matrix elements of D ij appear in d × d blocks corresponding to a d-dimensional IR. Note that Schur's first lemma 58 implies that all these remaining blocks are then a multiple of the unit matrix, if the form factors are well-behaved in the sense of Appendix C, where this statement is proven as Corollary C.3.
At l = (0, π) and (π, 0), the little group reduces to L l = C 2v and therefore, again by virtue of Corollary C.1, the five IRs of C 4v decay into three sets of form factors that do not mix with another. One contains s-and d x 2 −y 2wave, one d xy -and g-wave while the third one purely consists of p-wave form factors. This p-wave set splits into two, each transforming with a one-dimensional IR of C 2v . Altogether, this corresponds to the four IRs of the little group. For example, the most slowly varying pwave basis functions sin(q x ) and sin(q y ) then transform with two inequivalent one-dimensional IRs of C 2v , which may be referred to as p x -and p y -wave.
For bosonic momenta l on the boundary of the first Brillouin zone, i.e. for L l = C s , there are two such sets, one for s, p x/y , d x 2 −y 2 -and the other one for p y/x , d xy and g-wave. Again, those two sets correspond to the IRs of the little group. For bosonic momenta that do not lie on any of the symmetry axes, the little group just contains the identity element and all form factors may get mixed.
Let us now assume that the form-factors are wellbehaved in the sense of Appendix C. If the mixing between inequivalent IRs of C 4v is neglected, the bosonic propagators of the four one-dimensional IRs then inherit the full C 4v symmetry of the coupling function according to Corollary C.2. In contrast, the p-wave block still transforms with two-dimensional IR matrices. Let us note in passing that a mixing of different IRs has already been observed in Ref. 59 for the RPA pairing susceptibility at incommensurate Copper pair momenta.
So far, we have only considered a FFE in the Cooper channel. Of course, such an expansion may as well be performed in the other channels, which are then decomposed as If one wishes to simplify the RG flow equations by resorting to a FFE, the expansion has to be truncated behind a few terms in order not to exceed available computational resources. This may be conveniently done in the following way.
i.) Neglect the mixing between inequivalent IRs of C 4v (or the respective point group for other lattice geometries).
ii.) Only consider the most slowly varying form factor among equivalent IRs, i.e. neglect the admixture of higher harmonics.
(In addition to these approximations, the p-, d xy , and gwave sectors have not been taken into account in Ref. 39, as such form factors can be expected to play a minor role for the one-band Hubbard model at van-Hove filling.) If these approximations are adequate, the truncated FFE of the RG flow equations should in principle capture the important momentum dependences well. If, in contrast, the admixture of higher harmonics plays a role, a large number of bosonic channels might be needed. At least for the one-band Hubbard model, the above approximations seem to be fine in a large region of the parameter space. 42 The question now is, whether important orbital makeup effects are still captured within a viable truncation of an FFE.
Since the fermionic momenta are directly put on a grid in this work, we are in a position where we can easily keep track of mixing between inequivalent IRs. We expect this mixing to play a minor role, if the ordering vectors of leading and subleading instabilities are l = 0 or (π, π) or very close. By diagonalization of the coupling functions as matrices in q and q ′ with l fixed to one of these ordering vectors, optimized form factors can then be attributed to the respective instabilities. These optimized form factors will turn out to be close to the most slowly varying ones for most parameters considered in this paper, but in some cases also higher harmonics play a role. A sensible truncation of the FFE then consists in only retaining the terms corresponding to the most relevant eigenvalues. Clearly, the optimized form factors are scale-dependent in an fRG flow. In principle, it should be possible to parametrize this scale dependence. 60 Let us note in passing that similar effects have already been discussed within a Bethe-Salpether equation approach. 61 For incommensurate antiferromagnetism, however, the potentially non-zero mixing between inequivalent IRs of the point-group symmetry of the lattice may prohibit the calculation of an optimized form factor that is defined on the whole Brillouin zone. In such a case, a faithful truncation of the FFE would already contain too many terms to be numerically tractable. We will come back to the question of the applicability of a FFE when we discuss our numerical results in the next section.

IV. NUMERICAL RESULTS
We use the parameters displayed in Tab. I as a generic parameter set for which we run the RG flow in the conventional truncation, i.e. without the three-particle feedback. We then test the stability of our results against variation of some of these parameters and against threeparticle terms. We further compare the Emery model to the corresponding one-band Hubbard model. All results presented in this paper are for zero temperature and all parameter sets considered are on the hole-doped side.
Note that spectral functions are inaccessible in our approach, since we project to zero frequency and neglect the self-energy. But even if these approximations were relaxed within an instability analysis, the flow would be stopped at a finite scale, where contributions to the interaction start to diverge. Therefore, the spectral functions obtained from fRG and successive analytic continuation then would still not really comparable their DMFT counterparts in Refs. 27 and 30.
A. Nature of the leading instability For the generic parameters and at van-Hove filling, we observe a flow to strong coupling at about λ c = 2.9 · 10 −3 t pd , which roughly corresponds to 50 K for t pd = 1.3 eV. In order to determine the nature of this instability, we diagonalize the coupling functions Φ SC , Φ M and Φ K at the stopping scale. We attribute the largest of the eigenvalues of these three coupling functions to the leading instability, which is characterized by an optimized form factor given by the corresponding eigenvector.
For the parameter sets considered in the following, the most relevant eigenvalues in the pairing and the magnetic channels compete. Let us first describe our results for the generic parameter set of Tab. I. The optimized form factors for this parameter set are depicted in Fig. 3. In the pairing channel, contributions with total wavevector l = 0 dominate clearly. The optimized form factor corresponds to a d x 2 −y 2 -wave, with peaks that are a little broader than for f d = cos(q x ) − cos(q y ) (cf. Fig. 3). We will comment on the admixture of higher harmonics to this optimized form factor further below. In the magnetic channel, the optimized form factor corresponds to a deformed s-wave with small admixtures of higher harmonics. For the generic parameter set we find slightly incommensurate magnetic ordering vectors on the boundary of the Brillouin zone. Inequivalent IRs of C 4v may hence mix. In the case of our generic parameter set, we can indeed observe small admixtures of other IRs to the s-wave contributions of the optimized form factor. For example, a small d x 2 −y 2 -wave admixture is clearly visible.
Other types of instabilities such as a Pomeranchuk instability 8,62-65 and the formation of different types of loop currents 19,20,66 do not participate in the competition of the most relevant instabilities. For the former type of ordering, our results are in agreement with Ref. 65, where the experimentally observed nematic tendency [67][68][69][70] in cuprate materials appears to be of strong-coupling nature. Similarly, the absence of loop currents is not surprising in a weak-coupling scenario, since the critical interaction strengths are found to be quite large in meanfield calculations. 20 At van-Hove filling, one can expect that it is possible to tune the Emery model to a Stoner-like ferromagnetic instability by raising the value of t pp . For the parameter sets considered in this work, however, we do not find ferromagnetism to prevail over other ordering tendencies. We hence conclude that a Stoner instability only occurs for oxygen-oxygen hopping parameters that are far away from the generic value of t pp in Tab. I.
At the stopping scale, we also consider a FFE of the full coupling function V obtained within our new approach for comparison. We truncate this expansion behind the most slowly varying form factors. This is not to be confused with the results obtained from RG flow equations in an exchange parametrization, as we take mixing between different IRs of C 4v and of different harmonics within an IR into account in the integration of the flow equations. The propagator D d of d x 2 −y 2 -wave Copper pairs rescaled by the square strength of the fermion-boson interaction then reads as where f d (k) = cos(k x ) − cos(k y ). In the following, we will simply refer to D d as a bosonic propagator despite the rescaling by an energy-squared factor. Likewise, the propagator of the magnetic s-wave exchange-boson is obtained as with f s (k) = 1. Note that this FFE is only viable if the basis of Bloch states is properly chosen such that V is invariant under all point-group transformations of the lattice. In Fig. 4, D d (l) and M s (l) are depicted for the generic parameter set in Tab. I at the stopping scale. Both d x 2 −y 2 -wave Cooper-pair and magnetic s-wave propagators show peaks with values close to the corresponding eigenvalues of the coupling functions of the respective channel. Since the dSC peak is quite sharp while the incommensurate peaks of the magnetic propagator have a broader width, d x 2 −y 2 -wave superconductivity might prevail in a situation where the two most relevant instabilities are closely competing. In any case, the system is in a regime of two competing, mutually reinforcing instabilities. In this place, we would like to recapitulate that our stopping condition is quite weak. Therefore, if the magnetic and pairing channels are still of comparable strength at the stopping scale, these two channels are then indeed closely competing and this competition itself might have some physical content.
In Ref. 9, the parameter-space region of strong AFM-dSC competition was dubbed the saddle-point regime and interpreted as an insulating spin-liquid phase. It may also contain a region of homogeneous coexistence as described for the iron pnictides in Ref. 71. It is likely that a large part of this regime has a non-vanishing superconducting gap. Unfortunately, order parameters are not directly accessible within the present instability analysis. In a recent fRG approach to the one-band Hubbard model using rebosonization techniques 72 it has however been found that pairing is avoided inside the antiferromagnetic phase. This is not surprising, since at least parts of the Fermi surface are gapped away once spontaneous symmetry-breaking in one channel sets in, which hampers symmetry-breaking in the other channels. A very recent purely fermionic study on the fRG flow of the repulsive single-band Hubbard model into the superconducting phase is in full agreement with this picture. 73 In that work, a non-vanishing pairing gap is indeed found in a large subregion of the saddle-point regime. So a putative coexistence phase should be considerably smaller (µvH − µ) /t pd · 10 3 −15.4 −7.69 −3.85 0 3.85 7.69 n h · 10 2 6.5 8.7 9.9 11.5 13.1 14.3 than the saddle-point regime. In principle, a two-orderparameter mean-field approach is viable below the stopping scale, 10 but of course such a treatment is not free of bias.
As we shall find in the following, the character of the instability is quite robust against slight variations of the parameters. The stopping scale will turn out to be sensitive to the chemical potential µ and to the diagonal oxygen-oxygen hopping t pp further below. This latter dependence is already rather smooth, if the oxygen-oxygen hopping is not to far away from its generic value. So the behavior around the generic parameters suggests that the system is close to a first-order phase transition between AFM and dSC, as second-order transition would probably go along with a kink of the critical scale. Since the RG stopping scale is an upper estimate for the critical scale, such a kink may however be hidden. For the oneband Hubbard model close to van-Hove filling, a similar behavior has been observed. 9,10,39

B. Doping dependence
The system also stays in the saddle-point regime when the doping level is slightly varied. In Fig. 5, the stopping scale for the parameters in Tab. I is plotted as a function of µ (+-markers). In the following, we define the hole filling factor as the number of holes per unit cell and spin orientation that have been doped into the originally half-filled conduction band where ǫ(k) and Θ(x) denote the dispersion of the conduction band and the Heaviside step function, respectively. Note that n h may differ for the dispersions of the Emery model and the single-band models with effective parameters calculated according to Sec. II B. In the following, all values of n h will be for the Emery model. The doping level varies between n h = 6.5 · 10 −2 and 14.3 · 10 −2 and at van-Hove filling the filling factor is 0.115. We observe that between half-filling (n h = 0) and van-Hove filling, the stopping scale only varies slightly. The AFM and dSC instabilities are closely competing, except for n h < 0.09, where AFM clearly prevails. At hole doping beyond van-Hove filling, the stopping scale decreases rapidly and the tendency to d x 2 −y 2 -wave pairing gets a little stronger. Qualitatively, this behavior is analogous to the hole-doped one-band Hubbard model. 9 The enhanced d x 2 −y 2 -wave pairing at hole doping beyond van-Hove filling, however, comes with a broadening of the d x 2 −y 2 -wave form factor at the van-Hove points as can be seen from Fig. 6. Moreover, the magnetic s-wave propagator is of interest. In Figs. 7-9, we observe that for n h ≤ 0.09 the ordering vector l = (π, π) is commensurate. At about van-Hove filling, the commensurate peak of M s is split into four peaks at incommensurate ordering vectors. The deviation of the ordering vector from (π, π) then increases with hole doping, corresponding to the shift of Fermi surface segments at high density of states.
The highly incommensurate peaks at n h = 0.14 (i.e. at hole doping significantly beyond van-Hove filling) still allow the Kohn-Luttinger effect to generate an attractive d-wave pairing component. This however goes along with a deformation of the form factor f d for d-wave pairing (see Fig. 6). Such a simultaneous occurrence of incommensurability in the magnetic channel and a deformation of the fermion-boson vertex in the Cooper channel has already been observed in the one-band (t, t ′ ) Hubbard model 42 and may be explained as follows. Consider a singlet Copper pair with momenta (k, −k) scattered to (k ′ , −k ′ ) by the interaction in the Cooper channel which shall be mimicked by a one-loop particle-particle diagram with two spin-channel vertices. If these vertices have their peaks at transfer momentum l = Q = (π, π), the main contribution to the Cooper channel comes from k ′ = k. For incommensurate ordering vectors, the important contributions come from k ′ = k as well as from k ′ = k + Q i + Q j , where i and j may correspond to all possible combinations of the ordering vectors. The dependence of fermion-boson vertex in the Cooper channel on the fermionic momenta is hence smeared out around (0, π) and (π, 0) resulting in a broadening of the form factor. A shoulder-like broadening of the peaks of M s would give rise to the same effect in a similar way.
While higher harmonics do not contribute at the so-called anti-nodal points (0, π) and (π, 0), they may change the slope of the gap at the nodal points. Such an occurrence of multiple energy scales for the gap has been observed in Raman spectroscopy 74 and angle-resolved photoemission spectroscopy 75 experiments. In contrast to our results for weak coupling, Ref. 74 suggests decreasing contributions of higher harmonics with hole doping. This may be due to the strong-coupling nature of real cuprate materials.
As the Fourier transform of the form factor in the pairing channel corresponds to the distribution of the distance between the electrons forming a Cooper pair, deviations from the cos(q x ) − cos(q y )-form may also be analyzed in real space (see lower row in Fig. 6). One should, however, be aware that an interpretation in real space requires some care, since basis sets of Wannier functions may strongly differ in their localization properties (cf. Ref. 46). The position argument in the real-space form factors then corresponds to the relative distance between the constituents of such a pair, i.e. the two electrons or holes involved. Note that this distance can be resolved up to only n sites in all directions for 2n × 2n fermionic patches. We have therefore studied the flow for some parameters with a resolution of 8×8 fermionic and 24×24 or 120×120 bosonic patches away from and close to possible ordering vectors, respectively. The results are displayed in Fig. 6. We observe that the most important contribution corresponds to a cos(q x )− cos(q y )-form. But already at van-Hove filling, an admixture of higher harmonics is visible, which partly get shifted further away from the origin at hole doping beyond van-Hove filling. A thorough discussion of the minor contributions corresponding to Cooper pairing beyond nearest neighbors may require a resolution higher than 8 × 8 fermionic patches.
Before we analyze the impact of the t pp -hopping parameter, a remark on the effects of the coupling between the different channels seems to be in order. These effects go far beyond the spin-fluctuation induced generation of an attractive pairing interaction. In particular, if the magnetic propagator was calculated within RPA, i.e. if the Cooper and forward scattering channel were neglected in the flow, the stopping scale would be about one decade higher. Moreover, the magnetic propagator would be less sharply peaked. This behavior can be attributed to the feedback of the Cooper on the magnetic channel, which hampers antiferromagnetism before the d x 2 −y 2 -wave pairing interaction gets attractive. On a qualitative level, this effect is already captured in the two-patch approximation. 5 C. Dependence on hopping-between the p-orbitals So far, we have only investigated the impact of doping away from van-Hove filling, but not the interplay between µ and t pp . If, in this spirit, the hopping between the p-orbitals is now changed to t pp = 0.53 t pd , the effective second-neighbor hopping in the conduction band gets stronger and the tendency to d x 2 −y 2 -wave pairing should be enhanced. Indeed, the flow can now be attributed to the saddle point regime for all filling factors considered. The corresponding curve in Fig. 5 (x-markers), however, looks similar to the one for t pp = 0.50 t pd except for the insignificantly lower stopping scale. Moreover, at hole doping beyond van-Hove filling and at t pp = 0.53 t pd , the tendency to dSC is only slightly enhanced compared to t pp = 0.5 t pd . It therefore seems that a considerable region of the parameter space has to be attributed to the saddle-point regime as for the one-band Hubbard model in Ref. 9.
In Fig. 10, the dependence of the stopping scale on t pp is depicted for van-Hove filling both for 6 × 6 and for 8×8 fermionic patches. The curves for the two resolutions almost coincide, indicating that 6 × 6 fermionic patches are sufficient. The stopping scale behaves as follows: For the Emery model, we find a decrease of λ c with increasing t pp of less than one decade. Such a behavior is quite generic as, in the absence of orbital makeup, a more rounded Fermi surface depresses the stopping scale in other models. 16 Once t pp exceeds 0.54 t pd , the flow to strong coupling bears rather the signatures of pure d x 2 −y 2 -wave pairing than of the saddle-point regime. For the one-band Hubbard model with effective parameters, however, the situation is different: The stopping scales are much lower and drop abruptly as soon as the system enters the pure dSC regime at about t pp = 0.54 t pd . This is apparently caused by an abrupt growth of the hybridization at the saddle points which reduces the effective interaction strength U eff . Comparing the Emery model to its effective one-band Hubbard counterpart, we find that orbital-makeup effects in the Emery model partly counteract the decrease of λ c with a more rounded Fermi surface and that they prevent the stopping scale from dropping abruptly.
In Ref. 16, a similar behavior as been found for multiband models only involving orbitals on the Copper atoms. Also the DCQMC results of Ref. 33 support this conjecture. Namely, the critical temperature is found to increase with the value of t pp in that work. At first glance, this seems to contradict our fRG results. How-    not fully stringent, since, at strong coupling, mechanisms may be at work that do not occur in a perturbative picture.

D. Three-band vs. single-band model
Let us now look at multiband effects more systematically. So far, we have discussed results for the full Emery model as given by the parameters in Tab. I and variations of t pp and µ in a conventional truncation of the flow equations, i.e. without a three-particle term. If the feedback of this term is taken into account within the truncation proposed in Ref. 48, we do only find insignificant changes of the stopping scale. Moreover, the tendency to dSC is slightly enhanced between half-filling and van-Hove filling, so that the system stays in the saddle point regime down to n h = 6.5 · 10 −2 , as can be seen from the phase diagrams depicted in Fig. 11. Such minor modifications of the phase diagram seem quite surprising as those feedback terms had a great impact on the phase diagram of a two-band model in a two-patch approximation. Such a behavior of the three-particle feedback may, however, may be an artifact of the two-patch approach pursued in that older study. Namely, the two-patch approximation only allows for a small number of strong-coupling fixed points which results in a mutual exclusion of the Cooper and Stoner instabilities. Moreover, in the Emery model, the large gap between the conduction band and the valence bands results in a flat momentum structure of the diagrams corresponding to three-particle feedback terms. In a frequency-resolved study, the three-particle feedback may however play a more important role.
Before we elaborate on other multiband effects, the role of the oxygen p-orbitals in the two-particle term shall be discussed. First, we turn our attention to the impact of a U pp -term, which should be absent according to Ref. 45. We have varied U pp from zero to 0.1 U d . The stopping scale then only changes insignificantly and the systems stays in the saddle-point regime.
In a second step, the role of the interaction terms involving legs on the p-orbitals is now discussed, i.e. we compare the result for the full Emery model to what is obtained if all interaction terms except U d are ignored. Note that this interaction is still non-local in the band language, and that it is therefore different from an on-site Hubbard term. Data points for this level of approximation with dispersion parameters according to Tab. I are represented by stars in Fig. 5. These points almost coincide with the ones for the full interaction (+-symbols). Therefore interaction terms involving the p-orbitals only slightly affect the stopping scale.
However, this does not mean that these terms do not have any impact at all. Actually, the tendency to AFM is enhanced if the p-orbital interaction terms are neglected.
In particular, between half filling and van-Hove filling, the flow to strong coupling shows features of an AFMinstability rather than of the saddle-point regime. A form factor deformation above van-Hove filling still occurs, but this effect is a little weaker without U p and U pd . One may now try to understand the enhancement of the AFM tendency by considering the two-patch couplings g i in the ultraviolet. For the full interaction g 1 , g 2 and g 4 have the same value of about 0.269 t pd while g 3 is lowered by a few percent to 0.264 t pd . For an interaction involving the d-orbital only, in contrast, g 1 = g 4 = g 2 = 0.255 t pd and g 3 is now enhanced to 0.257 t pd . So we have an overall decrease of the two-patch couplings and their detuning differs from the case of full interaction, resulting in a relative increase of the d-wave pairing and the AFM components of the interaction in the two-patch approximation, g 3 − g 4 and g 1 + g 3 , respectively. Apparently, a lowering of the stopping scale induced by the former is compensated by the latter and, as a net effect, the AFM tendency gets stronger.
We now continue with the discussion of the other multiband effects that have been listed in Sec. II B. First, we replace the two-particle interaction for the conduction band by an on-site Hubbard interaction while we leave the dispersion unchanged. Its strength U eff is chosen to be the average of the two-patch couplings that correspond to the full interaction. From the open squares in Fig. 5, we infer that the stopping scale is significantly lowered in this approximation. As the detuning of the two-patch couplings is rather small on this level of approximation, we conclude that phase-space regions away from the saddle points play some (minor) role even at van-Hove filling. Since the hybridization of p-and dorbitals is strongest at the saddle points, the contributions of those regions to the diagrams on the right-hand side of the flow equations are underestimated by an onsite interaction with strength U eff . Therefore a flow to strong coupling occurs at lower scales.
One may therefore wonder, whether another prescription for choosing U eff might give results that are more close to those for the three-band model. Since the hybridization between d and p-orbitals is strongest at the saddle points and since the U p and U pd interactions are of minor importance, the choice U eff = U d seems appealing as well. If one were to follow this alternative prescription, the stopping scale of the effective model would overshoot the value for the original model by roughly a factor of two. This is not surprising, since the averaged interaction strength of the two models then already differs at the saddle points. Moreover, the prescription U eff = U d can be regarded as the leading-order result of a gradient expansion around k = (0, 0). However, this is inconsistent with the expansion around the saddle points of the conduction band underlying the calculation of the effective hopping parameters t and t ′ . We therefore continue to use the prescription given in Sec. II B, since the approach seems to be the most systematic one.
The suppression of orbital makeup then generically lowers the stopping scale. Away from van-Hove filling (in particular at hole doping beyond), this lowering is more pronounced as the gradient expansion gets worse. Moreover, the system stays in the saddle-point regime for all filling factors considered. There are now several possible mechanisms giving rise to the enhanced tendency to dSC between half filling and van-Hove filling. As the attractive d x 2 −y 2 -wave pairing component is generated by fluctuations in the magnetic channel, this enhancement may simply be caused by the lowering of the stopping scale. Furthermore, a detuning of the two-patch couplings that hampers the dSC instability is now absent. Let us note in passing that the broadening of the d x 2 −y 2 -wave form factors around the saddle points is still restricted to hole doping beyond van-Hove filling. Finally, we consider the t, t ′ , U eff one-band Hubbard model (filled squares in Fig. 5), i.e. we now approximate also the dispersion by the leading-order result of a gradient expansion around the saddle points. Compared to the previous data set, the stopping scale is again lowered by almost a factor of two at van-Hove filling. Away from van-Hove filling, this depression of λ c again grows. At hole doping beyond van-Hove filling, this effect is more pronounced between half filling and van-Hove filling and at about n h = 0.14 the system enters the pure dSC regime. Again, a form factor deformation occurs at van-Hove filling and larger hole doping. Compared to the original model, the stopping scale is a factor between five and ten too low. We therefore conclude that in the effective action for the conduction band long-range hopping terms play a significant role, since they enhance the stopping scale.

V. DISCUSSION AND OUTLOOK
In this article, we have studied the hole-doped threeband Emery model at weak coupling by considering the RG flow of an effective interaction for the conduction band. This has been done within an improved channeldecomposed approach, which considerably improves the momentum dependence of fermion-boson vertices of previous treatments of the Hubbard model within an exchange parametrization. In Sec. III C, we have classified the contributions that are present in our approach, but potentially neglected in Refs. 39-41, albeit without qualitative changes of the leading instabilities. As new results we present tentative phase diagrams for the Emery model, the wavevector structure of the effective interactions, in particular in the spin channel, and the structure of the (deformed) d-wave pairing gap. We also provide a detailed comparison between the RG flow in the threeband case in one-band models with comparable Fermi surfaces.
In the Emery model for the parameters considered here, the leading instabilities are d x 2 −y 2 -wave pairing and antiferromagnetism. The leading instability often can hardly be distinguished from the subleading one, as also found in the one-band case 9 . For most of the parameter sets considered, the system appears to be in some intermediate region between the regimes with clear AFM and dSC instabilities, where both tendencies are comparably strong even very close to the critical scale. This so-called saddle point regime might contain a coexistence phase and should be superconducting to a large part. Signatures of this saddle-point behavior survive at doping levels a few percent away from van-Hove filling, again similarly to Ref. 9. Between half filling and van-Hove filling, the stopping scale shows a plateau. Once the system is hole-doped beyond van-Hove filling, the stopping scale decreases rapidly. This goes along with an increasing incommensurability of the AFM ordering vector giving rise to a broadening of the superconducting gap at the saddlepoints of the dispersion. In real space, this corresponds to Cooper pairing between electrons in Wannier states that are centered around points in non-neighboring unit cells. Such effects could not be explored with previous Npatch fRG schemes, but are now available by virtue of the refined wavevector resolution of the channel-decomposed vertices. They are not intrinsically of multiband nature, and can also be observed in the single-band model. All these properties are quite robust against the omission of interaction terms involving the oxygen p-orbitals.
The effective second-nearest-neighbor hopping t ′ of the conduction band can be tuned by changing the oxygenoxygen hopping t pp of the Emery model. An increased value of this parameter therefore leads to a more rounded Fermi surface at and close to van-Hove filling. We observe that the stopping scale then decreases only slowly with increasing t pp . For the corresponding single-band Hubbard model, we find a decrease of the stopping scale with t pp as well, but now with a peculiar feature: Between t pp = 0.53 t pd and 0.54 t pd , where the copper-oxygen hybridization in the full Emery model rises abruptly at the van-Hove points, the stopping scale suddenly drops by a factor of two. Therefore, orbital-makeup effects included in the three-band model apparently counteract and almost compensate the effect of a more rounded Fermi surface, similarly to what was found in Ref. 16. In other words, not only the Fermi surface shape and the density of states matters for the energy scale of those weakcoupling instabilities.
Close to van-Hove filling, it seems reasonable to determine the parameters of the effective one-band Hubbard Hamiltonian in Eq. (3) from a gradient expansion around the saddle points of the dispersion, as the C 4v symmetry shared by both models in a suitable Bloch basis allows for such a procedure. In order to obtain a single parameter for the interaction, the effective on-site interaction is chosen as the average of the four two-patch couplings, which appear in leading order in a gradient expansion and would constitute the running couplings in the two-patch approximation. Except for a detuning of the two-patch couplings, the single-band and the three-band model coincide at the saddle points of the dispersion.
We find that the differences away from the saddle points are crucial in the sense that the Emery model formulated in the band language contains hopping and interaction terms connecting non-neighboring unit cells, which turn out to enhance the stopping scale. Yet on a qualitative level, the one-band Hubbard model still has a similar phase diagram. On a more quantitative level, these longer-ranged terms play a role both in the dispersion, where hopping between non-neighboring unit cells has some impact, and in the interaction, which is decorated by orbital makeup. For the effective Hubbard model, the lowering of the stopping scale enhances the tendency to d x 2 −y 2 -wave pairing, whereas the detuning of the two-patch couplings in the original model is of minor importance. This is in analogy to recent VCA results 37 for the strong coupling case, where the effective one-band Hubbard model also accounts for the more universal features of the phase diagram. Of course, a different prescription for calculating the parameters of the effective one-band Hubbard model could have been chosen in the present work. The one presented, however, appears to be rather systematic. It allows one to determine the origin of the deviations from the original model which should be related to k-space regions away from the saddle points. More precisely, by using the interaction values a the saddle points for the effective one-band model, the interaction away from the Fermi surface and away from the saddle points is underestimated. This happens because the hybridization of the d x 2 −y 2 -orbital with the p-orbitals is largest at the saddle points. The relatively strong impact of the other regions on our results therefore suggests that those degrees of freedom are important in the sense that their presence does not only allow for a rich fixed point structure of the RG flow, but that they also directly influence the phase diagram.
A significantly large region of the parameter space constitutes the saddle point regime, where the dSC and AFM instabilities are closely competing and mutually reinforcing. The present instability analysis lacks directly accessible measurable quantities in potentially symmetrybroken phases which would facilitate the interpretation of such a behavior. In order to find a sharp phase boundary or a coexistence phase between AFM and dSC, it would be advantageous to enter the symmetry-broken phases within a purely fermionic approach. Furthermore, other, more exotic possibilities such as a Fermi surface truncation without long-range order 5,9 are also inaccessible without proper flows for the self-energy. Including these effects together is a formidable task beyond the current frontier in RG methods. We therefore have to refrain from explicitly breaking the U(1) and/or SU(2) symmetries, mentioning that fermionic flows have been continued into the superfluid phase 73,76-79 and that we have addressed parametrization questions and mean-field models for the AFM phase in a recent publication, 80 while we are currently studying the AFM phase beyond mean-field models within a channel-decomposed approach. Moreover, as recently suggested by Giering and Salmhofer,41 the parameters of an effective partially bosonized theory may be derived within a purely fermionic RG flow, whereas symmetry-broken phases are entered within a mixed flow.
Such issues left aside, the strong-coupling nature of abinitio parameter sets still prevents us from thoroughly discussing the applicability of our weak coupling approach (with interaction parameters one decade smaller than typical literature values) to real cuprate materials. We therefore refrain a discussion of some parameter trends observed in the strong-coupling literature (see, for examples, Refs. 31 and 33). However, the Emery model at weak coupling is shown to have the same leading instabilities than the one-band Hubbard model, complementing a strong-coupling VCA study by Kiesel et al. 37 and a DMFT study by de' Medici et al. 27 Moreover, it has turned out to form a good, rather simple testbed for the new discretization scheme presented here. This scheme can of course be carried over to other, more complicated multi-band systems where a weak-coupling approach is indeed realistic. Candidates are the iron superconductors or strontium ruthenates (for a recent SMFRG study, see Ref. 81), both with three and more Fermi surfaces. We are also glad to acknowledge the friendly support provided by the highperformance computing team at RWTH. This work was supported by the DFG priority program SPP1458 on iron pnictide superconductors and by the DFG research unit FOR 723 on functional renormalization group methods.
Appendix A: Point-group symmetries of the effective one-band action Although this might be hard to see at first glance, the Emery model possesses a hidden C 4v -symmetry. In Ref. 46, two of us have discussed how point-group symmetries manifest themselves in a large class of multiband models. Also the Emery model belongs to that class. Namely, there exists a three-dimensional representation {MÔ}Ô ∈C4v such that the Emery Hamiltonian (including the interaction terms) is invariant under the substitution Here, the RÔ denote the two-dimensional rotation matrices, which form a faithful representation of C 4v . For a coordinate exchangeÎ ′ and a reflectionÎ with respect to the x-axis, the respective three-dimensional represen-tation matrices read as Since all point-group operationsÔ ∈ C 4v can be written as products of the identity operation,Î andÎ ′ , the other MÔ follow from the group law MÂ = MĈ MB forÂ = BĈ. Note that the representation matrices MÔ decay into irreducible blocks -a one-dimensional A 1 block for the Copper d-orbitals and a two-dimensional E block for the oxygen p-orbitals.
In our fRG approach, the model is expressed in a band language, i.e. in terms of new pseudo-spinors with a unitary, wavevector-dependent u(k), where u(k) H 0 (k) u † (k) is diagonal. Of course, we would like to exploit the hidden C 4v -symmetry of the model in the numerical integration of the RG flow equations. As shown in Ref. 46, this can be accomplished as follows. There is some freedom in the transformation from orbitals to bands, corresponding to different choices of the phase of the eigenvectors in the lines of u(k). These phases can be chosen individually for each wavevector k and for each band. Since in the second-quantized language kdependent field operators create Bloch states from the vacuum state, this phase plays the role of a global phase of these Bloch states. In a so-called natural Bloch basis, one has u(RÔk) = u(k) MÔ ∀ k = RÔk .

Consequently, expressed in such a basis, the Hamiltonian is invariant under
and hence the coupling function U (k 1 , k 2 , k 3 ) in Eq. (4) shows a trivial point-group behavior, i.e.
This property is as well inherited by the coupling functions V (k 1 , k 2 , k 3 ) of the renormalized interaction in Eq. (B1) and U β (k 1 , k 2 , k 3 ) which enters the threeparticle feedback in Eq. (B 2).
In this paper, we therefore work in a natural Bloch basis. According to Ref. 46, discontinuities in the interaction then only occur in the terms involving the p-orbitals, which will turn out to be of minor importance. exploits U(1), SU(2) and translational symmetries. The quartic part Γ (4) of the (scale-dependent) oneparticle irreducible functional in the conduction band is parametrized in terms of a coupling function V (k 1 , k 2 , k 3 ) according to

Conventional truncation
In the conventional truncation, three-particle and higher vertices are neglected in the flow. (For a derivation see Refs. 7 and 49.) The scale derivative of the coupling function V consists of three parts The particle-particle contribution and the crossed particle-particle part can each be represented by one diagram with one loop containing (scale-dependent) propagators G(k) on the conduction band. Since self-energy effects are neglected, we have where ǫ 1 denotes the energy of the conduction band and where the Ω-scheme regulator R λ is given in Eq. (5).
Note that an RPA resummation in the Cooper or particle-hole channel is equivalent to an RG flow in which all terms except T pp or T cr ph , respectively, are neglected. (The corresponding Bethe-Salpether equation is then equivalent to the flow equation for the two-particle vertex. 39 ) Vertex corrections and screening, however, are accounted for by the direct particle-hole diagrams T d ph = dp [∂ λ G(p) G(p + k 3 − k 1 )] × [2V (k 1 , p + k 2 − k 3 , p) V (p, k 2 , k 3 ) − V (k 1 , p + k 2 − k 3 , k 1 + k 2 − k 3 ) V (p, k 2 , k 3 ) −V (k 1 , p + k 2 − k 3 , p) V (p, k 2 , p + k 2 − k 3 )] .
If all external frequencies are projected to zero, the Matsubara sums in the loops can straightforwardly be evaluated analytically within the Ω-scheme using contour techniques. The initial condition V = U is given by the quartic part of the bare action for all three orbitals, where X i = (k i , σ i , α i ) is a short-hand notation for the quantum numbers of the fields and where the α i denote orbital indices. The coupling function F may be parametrized as in Eq. (B1) F (X 1 , X 2 , X 3 , X 4 ) = [W α (k 1 , k 2 , k 3 ) δ σ1,σ4 δ σ2,σ3 − Wα(k 2 , k 1 , k 3 ) δ σ1,σ3 δ σ2,σ4 .
Terms that have the valence-band index on another leg than the first can be reconstructed from this coupling function by exploiting the Pauli principle and the particle-hole symmetry of the interaction. When the valence bands are integrated out, a three-particle term is generated. In this Appendix, we will only give the resulting flow equation for the three-particle feedback and omit intermediate steps.
(For a more detailed discussion of the case of one valence band, we refer to Ref. 48, as the generalization to a larger number of valence bands is straightforward.) These three-particle feedback terms are now given separately as corrections R pp , R cr ph and R d ph to the particleparticle and the crossed and direct particle-hole terms T pp , T cr ph and T d ph , respectively.
In the band-mixed loops, G β (q) = [iq 0 − ǫ β (q)] −1 denotes the propagator of for valence-band hole with band index β = 2, 3 and the single-scale propagator S(q) is given by the scale-derivate of the propagator G(q) for the conduction band. Zero-frequency projection allows again for an analytical evaluation of the Matsubara sums.
Corollary C.1 (No mixing) Let P (l) be a bosonic propagator that has been projected out of a G-symmetric coupling function Φ according to Eq. (C1). Suppose that, for fixed l, P (RÔl) = P (l) ∀Ô ∈ K, where K is a subgroup of G. If the basis set of form factors is then organized in blocks corresponding to IRs of K, blocks that mix form factors of inequivalent IRs of K must vanish.
Proof Consider the projection rule (C1) for P ij (RÔl) and matrices M α O of the αth IR of K which transform the form factors according to We further label the IR of K that transforms f i with α and the one transforming f j with β. After substituting the integration variables p and q by RÔp and RÔq, respectively, we exploit the point-group symmetry of Φ. For the block P αβ (l) relating the αth and βth IR, we then find Since the left-hand side equals P αβ (l) according to the premise, Schur's second lemma applies. Therefore, if the αth and βth IR of K cannot be related by an equivalence transformation, the block P αβ must vanish.
In Sec. III C, we have used Corollary C.1 for K = L l , L l being the little group of the bosonic momentum l. If the premise P (RÔl) = P (l) ∀Ô ∈ K were satisfied as a consequence of some approximation, this approximation neglects the mixing of inequivalent IRs of K.
For the next two corollaries, the following definition appears useful.
Definition We call a set of form factors well behaved under a point-group K if its elements transform according to either identical or inequivalent IRs of K consisting of unitary matrices.
The point-group behavior of the bosonic propagator is governed by the following law.
Corollary C.2 Consider again a bosonic propagator P (l) that has been projected out of a G-symmetric coupling function Φ with form factors which are well behaved under G. Blocks P αβ relating equivalent onedimensional IRs α and β of G then are fully G-symmetric, i.e. P αβ (RÔl) = P αβ (l) ∀Ô ∈ G for arbitrary l.
Proof We observe that, also in the present case, Eq. (C2) holds, with α and β now labeling IRs of the full point group G. This equation is trivially fulfilled if a block vanishes according to Corollary C.1. For non-vanishing blocks, well-behaved form factors give rise to M α O = M β O . If α and β then label one-dimensional IRs of G, these matrices are just phase factors, which cancel.
If the mixing of inequivalent IRs of G is neglected, the remaining one-dimensional irreducible blocks of P (l) are hence G-symmetric.
Finally, Schur's first lemma directly gives rise to the following corollary.
Corollary C.3 Suppose that, for fixed l, P (RÔl) = P (l) ∀Ô ∈ K, where P (l) is a bosonic propagator obtained from a G-symmetric coupling function and where K is a subgroup of G. For a well-behaved set of form factors under K, the non-vanishing irreducible blocks P αβ (l) of P (l) then are a multiple of a unit matrix, with α and β labeling IRs of K.
Proof Again, Eq. (C2) holds. For well-behaved form factors, the representation matrices M α O and M β O witĥ O ∈ K are equal to one another, if P αβ does not vanish according to Corollary C.1. Since P (RÔl) = P (l) ∀Ô ∈ K, Eq. (C2) simply states that M α O and P αβ commute ∀Ô ∈ K. According to Schur's first lemma, P αβ then must be a multiple of a unit matrix.