Instabilities of quantum critical metals in the limit $N_f\rightarrow0$

We study a model in 2+1 dimensions composed of a Fermi surface of $N_f$ flavors of fermions coupled to scalar fluctuations near quantum critical points (QCPs). The $N_f\rightarrow0$ limit allows us to non-perturbatively calculate the long-range behavior of fermion correlation functions. We use this to calculate charge, spin and pair susceptibilities near different QCPs at zero and finite temperatures, with zero and finite order parameter gaps. While fluctuations smear out the fermionic quasiparticles, we find QCPs where the overall effect of fluctuations leads to enhanced pairing. We also find QCPs where the fluctuations induce spin and charge density wave instabilities for a finite interval of order parameter fluctuation gaps at $T=0$. We restore a subset of the diagrams suppressed in the $N_f\rightarrow0$ limit, all diagrams with internal fermion loops with at most 2 vertices, and find that this does not change the long-range behavior of correlators except right at the QCPs.


Introduction
Interacting fermions at finite density form a rich physical system that can be used to describe many condensed matter phenomena. Such systems are generally governed by the Fermi liquid IR fixed point but are susceptible to several symmetry breaking instabilities resulting in states with superconductivity (SC), ferromagnetism (FM), a Fermi surface breaking rotational symmetry, or a state with charge or spin density waves (CDW, SDW) breaking translational symmetry.
Transitions to these states can sometimes be tuned to zero temperature in condensed matter systems by changing e.g. pressure or the concentration of different dopants. When such a transition is second order, it is associated with strong fluctuations of the order parameter field near the zero temperature quantum critical point (QCP). The interaction between these order parameter fluctuations and the fermions is relevant in two dimensions and may thus alter the Fermi liquid IR fixed point. More dramatically, the interaction with near-critical fluctuations may themselves induce an instability to a symmetry broken state. This makes interacting fermions in two dimensions near a QCP an even richer system to study, we refer to them as quantum critical metals. The modified fixed points are called non-Fermi liquids (NFLs) or strange metals. The sharp quasiparticles of the Fermi liquid are destroyed but the Fermi surface remains, they are still metallic but their scaling laws are modified.
Many experimental systems of interacting fermions that do not fit into a Fermi liquid description are in fact (quasi) two-dimensional and close to QCPs. Examples are cuprate [1,2] and iron pnictide [3,4,5] superconductors showing NFL physics in the normal phase and critical temperatures higher than what can be explained within Fermi-liquid theory. Critical fluctuations can lead to NFL physics but the destruction of the quasiparticles is thought to limit pairing. At the same time, certain critical fluctuations can work as a pairing glue that enhances superconductivity. This makes it important to study quantum critical metals to understand which of these competing effects dominates and whether the vicinity of the QCPs can explain the enhanced superconductivity.
A fermion-boson model [6,7,8] is typically employed to study quantum critical metals. The fermions represent low-energy excitations at a Fermi surface and the boson takes the role of the critical fluctuations of the order parameter. The symmetry broken at the QCP is reflected in the symmetry of the boson. There is an important distinction between transitions breaking translational symmetry where the order parameter fluctuations have a finite momentum Q and those that do not break translational symmetry and the order parameter fluctuations are centered around Q = 0. We only consider the Q = 0 case in this work. As mentioned, these are strongly coupled systems in two dimensions and the fermion-boson model is thus not amenable to perturbation theory. As a consequence, several approximation methods have been employed. A common approach has been to extend the theory to get a new expansion parameter that is taken to a limit where we can do calculations. The resulting system is not the one we are ultimately interested in, but the hope is that some of the interesting physics remains. One example of this is to not study the model in 2 dimensions, but in 3 − dimensions [9,10,11]. is treated as a small perturbation away from the upper critical dimension where the theory can be studied perturbatively.
Another approach is to change the field content of the model. This can be done in several different ways. Instead of considering just a single fermionic field we may consider N f flavors of fermions all coupling to the same bosonic field with a rescaled coupling constant g = λ/ N f . This extended model has been studied in the limit of large N f [12,13]. It was later found that the scaling laws presented in these works do not survive when adding higher loop corrections [14,15].
A related approach is the matrix large-N limit in which the boson is an N × N component matrix transforming in the adjoint representation of a global SU (N ) symmetry group and the fermion an N component vector transforming in the fundamental. The coupling is scaled as g = λ/ √ N and the theory is studied in the limit of large N . All fermion loops are suppressed in this limit as well as all non-planar diagrams. This means that the boson receives no corrections from the fermion and only a subset of the boson corrections to the fermion are kept. Polynomially many diagrams contribute at each order as opposed to the factorially many in the full theory. Both non-Fermi liquid physics [16,17,18,19] and the superconducting instability [20,21,22] of quantum critical metals have been studied in this limit.
Finally we arrive at the small N f limit which is the approach considered in this paper. This is set up similar to the large N f limit, but we do not rescale the coupling constant and we bring N f to 0 instead of infinity. It may seem awkward to set N f to zero, in the large N f and N cases above we have a sequence of physical theories with integer numbers of fields that approach the limit under study. Here there is no such sequence of physical theories but in practice there is no difference. All diagrams come with an order of N f and we keep the lowest order instead of the highest. It turns out that the N f → 0 limit keeps all of the diagrams kept in the matrix large-N limit and it additionally keeps the crossed diagrams and thus contain factorially many diagrams at each order. Only the fermionic loops are suppressed by taking N f → 0. The momentum space retarded two-point function was calculated in the small N f limit [23]. Surprisingly, it was possible to obtain an (almost) closed form expression, and it showed that the fermion dispersion becomes non-monotonic due to corrections from the critical fluctuations. This was then extended in [24] where a framework was developed to calculate general fermion n-point correlation functions in the small N f limit of a fermion-boson model.
In addition to keeping a superset of all diagrams of the matrix large N limit, the small N f limit allows us to calculate explicit correlation functions, albeit in real space and for distances longer than 1/k F . This provides a useful complement to the results that have been obtained using other limits that are mostly in the form of scaling laws and renormalization group flows.
The small N f calculations are done in a long distance limit where we consider correlators between local operators separated far compared to the inverse Fermi momentum. This lowenergy limit does not commute with either the matrix large-N limit nor the small N f limit. An important effect that is missed by taking these limits before the low-energy limit is called Landau-damping. The boson gets dressed by fermion bubbles that change the IR scaling of the boson. There are earlier works in the small N f limit where the authors incorporated these effects from the start by using an explicitly Landau-damped boson when calculating the real space fermion two-point function [25,26,8,27,28]. This takes into account some of the low-energy effects missed in taking the N f → 0 limit first. However, this was not done in a systematic way. Landau-damping effects were included systematically (with further constraints on the order of other limits involved in the definition of the theory) in the particular limit N f → 0, k F → ∞, N f k F = constant in [29]. The momentum space fermion two-point function was calculated there as well and it was found that the non-monotonicity of the fermion dispersion disappeared for a non-zero N f k F .
The above works in the small N f limit studied non-Fermi liquid physics of quantum critical metals in the normal phase. The cuprates and iron pnictides additionally show enhanced superconductivity compared to what would be expected from Fermi liquid theory and this has not yet been studied in the small N f limit. In this paper we explore how critical fluctuations can lead to instabilities of finite density fermions in the small N f limit. We do this by using the framework developed in [24] to calculate correlation functions of several fermion bilinears. As opposed to the fermion self-energy studied in earlier works, these higher order n-point functions are sensitive to the symmetry of the fermion-boson coupling so we consider some different types of QCPs and additionally fermions coupled to a U (1) gauge field which can be described by the same fermion-boson model. Experimental systems are always at a finite temperature, and it has been pointed out that non-Fermi liquids may behave quite differently at finite temperatures compared to at T = 0 [30,31,32]. Because of this, we go beyond earlier small N f works by working at both T = 0 and finite T . Additionally, we depart from the critical point and also consider gapped order parameter fluctuations to obtain a full phase diagram on the disordered side of the QCPs.
The paper is organized as follows. In Section 2 we present the model we have studied and the particular limit this is studied in. We start out with a very general model to first explore the possible phases without restrictions. In this section we also extend the framework in [24] to account for momentum dependent couplings. In Section 3 we calculate pair, charge and spin correlation functions of the general model in search of unstable modes. We consider the effect of Landau-damping corrections by redoing this in the double limit of [29] in Section 4. In Section 5 we apply these results to specific physical systems that can be described by the model that we have been exploring. We consider charge and spin nematic transitions, ferromagnetic transition, circulating current transition and coupling to an emergent U (1) gauge field, in different subsections. We also consider the case of fermions coupled to multiple near-critical fluctuations. Each of these subsections are ended with a discussion where we compare our results to earlier findings on similar models. Finally in Section 6 we summarize our work and pose some open questions for the future.
2 Setup and calculation of fermion n-point functions in the Similarly to [13,21] we consider a general theory that has applications to several systems composed of a Fermi surface coupled to Q = 0 gapless bosonic excitations. The theory contains free parameters whose values depend on the particular application we consider. For now, we keep the values of these parameters unspecified. In Section 5 we consider physical realizations of this model with more concrete parameters. We consider a Fermi surface of N f flavors of spin-1/2 fermions ψ i,σ and order parameter fluctuation fields φ a : i is the fermion flavor index and takes values i = 1, ..., N f , σ is the fermion spin σ =↑, ↓. The fermions transform in the fundamental representation of a global U (N f ) flavor symmetry group. a = 1, ..., N b enumerates the order parameters/order parameter components. r a is a tuning parameter across the QCP where φ a = 0 for r a < 0 and φ a = 0 for r a > 0. For simplicity, we limit ourselves to the case of a circular Fermi surface (k) = k 2 /2m. We define the bare Fermi momentum k F = √ 2mµ and bare Fermi velocity v F = 2µ/m. The fields φ a couple to the operator: where we use the notation kx = −ωτ + k · x and the Fourier transformed fields: The coupling function λ a,σ (k) characterizes how the fermion couples to the order parameter fluctuations. Only the dependence on the direction of k is relevant for low energy excitations close to the Fermi surface. A φ 4 term can be considered irrelevant and is left out of the action 1 . We consider dynamic bosons but the static case of our results is obtained by the and the limit c a → ∞.
With an eye towards finding instabilities to pairing and charge/spin order we consider correlation functions of pair creation/annihilation operators and charge/spin density operators: where ρ σ is invariant under the global U (N f ) whereas b is not. We consider real space correlators of these operators: Writing the charge and spin correlators in terms of ρ ↑ , ρ ↓ correlation functions we have: The above operators all contain two fermionic fields. This, together with the N f → 0 limit, means that order by order the contributing diagrams only contain two continuous fermionic lines, each starting and ending at the above bilinears. There is only one way of connecting these lines in the case of the pair operators and the ρ ↑ (0)ρ ↓ (x) correlator but for the ρ α ρ α correlator these lines can be connected in two different ways. The fermionic lines attached to the density operators can either connect back to the same operator or connect the other insertion, see Fig. 1. We get an additional sum over fermionic flavors when the propagators connect back to the same insertion as in Fig. 1(b). This means that those contractions are subleading in the small N f limit. The ρ ↑ ρ ↓ correlator only has that type of contraction and will always be subleading: We can thus not differentiate between charge and spin density fluctuations at leading order in the N f → 0 limit. Let us now see how we can calculate these correlation functions. We start off by writing a generating functional with sources J σ,i , J † σ,i for the fermionic fields: Next we integrate out the fermionic fields. To do that we write out the interaction in terms of real space fermionic fields: The fermionic integral is Gaussian and the generating functional is obtained as The background field Green's function is diagonal in spin and flavor indices and independent of flavor so we can write it as G σσ ,ij [φ](z, z ) = δ σσ δ ij G σ [φ](z, z ). The background field Green's function is the solution to and In momentum space: we have Following the procedure of [24] we expand around a pointnk F on the Fermi surface. Now we additionally expand the functions λ a,σ (k) = λ a,σ,n + O(k −nk F ), and Fourier transform back to real space: The Euclidean time variables τ 1 and τ 2 are constrained to the interval (0, β) when the theory is considered at a finite temperature. The fermionic Green's function G σ [φ] has antiperiodic boundary conditions and the bosonic fields φ a have periodic boundary conditions in the τ direction. The solution to Eq. (24) with antiperiodic boundary conditions is given by where and f ± n (x) are the β-periodic 2 (+) and β-antiperiodic (−) (in τ ) solutions of These functions can be expressed as: We define the cross product as a × b = a x b y − a y b x . The found solution G σ,n [φ](x 1 , x 2 ) is only valid when acting on momentum modes close to k Fn . As in [24] we project onto these modes and integrate over different directionsn(θ) to obtain an operator G IR [φ](x 1 , x 2 ) that is valid as long as We consider k F v F T, |x 2 − x 1 | −1 and a field φ without momentum components of order k F . The η integral can be viewed as a Fourier transform to the variable k − k F . If we assume v F |τ 2 − τ 1 | k −1 F then the η dependence of the rest of the integrand has no momentum components of order k F and this Fourier transform will be small unless |k − k F | k F . We later comment on what happens for v F |τ 2 − τ 1 | ∼ k −1 F . For the leading behavior we can thus assume k is of order k F . We then see that the exponent ikn(θ) · (x 1 − x 2 ) makes the θ integral oscillate rapidly except at the two points wheren is parallel or anti-parallel to x 12 = x 2 − x 1 . We perform saddle-point approximations around these points Following the reasoning in [24] we perform the k and η integrals to obtain: We would also like to consider equal-time correlation functions so now we consider the case where v F |τ 2 − τ 1 | ∼ k −1 F . Consider Eq. (30). The η integral now gets a contribution also for |k − k F | ∼ k F since the fraction has frequencies of order k F when η ∼ v F |τ 2 − τ 1 |. To calculate this contribution we make use of Since v F |τ 2 −τ 1 | ∼ k −1 F and φ contains no scales of order k F the first term can be neglected. The second term can also be neglected since η ∼ v F |τ 2 − τ 1 | for the contribution we are interested in. We can set φ = 0 to leading order. Note that this is only true for the contribution missed in the saddle-point approximation. We thus have Finding the free propagator in real space we find that the last two terms above cancel and the saddle point solution actually works for small τ 2 − τ 1 as well. Now we consider two-point functions of the composite operators. Differentiating the generating functional with respect to the sources and only keeping the leading contribution for small N f we have: The determinant action has now been omitted since it is subleading in N f . Note that the density correlator contains one more contraction but it is also subleading in N f and can be ommited, see Fig. 1. Note that from here on, in "subleading" we include terms that are either subleading in the small N f limit or in the x ≡ |x| k −1 F limit. As in [24] we now integrate out the fields φ a to obtain: Here the functions h ± a (τ, x) show up when contracting the sources for the fields φ a . The sources comes from G σ,IR [φ a ] in Eq. (35,36). h + a (τ, x) (h − a (τ, x)) comes from contractions with the same (opposite) signs of s in Eq. (32). h ± a (τ, x) is defined as The denominators come from writing (28) in momentum space. The sum is over bosonic Matsubara frequencies ω n = 2πnT . D a (ω n , k x , k y ) are the φ a Green's functions which for the action in Eq. (1) are given by:

Divergences in long distance correlation functions
At the critical points T = 0, r a = 0 there is no internal scale in (39) so by dimensional analysis we have For τ = 0, the h ± a functions are linear in the separation x and the prefactors of x in the exponents of Eq. (37) and Eq. (38) are quadratic forms of the different coupling constants λ ± ↑,↓ . It is possible to find couplings such that correlations grow exponentially in the spatial separation x unless these quadratic forms are negative semi-definite. An unbounded growth is unphysical so it is interesting to investigate whether it appears for any setup.
In the introduction we mentioned that quantum critical fluctuations might induce instabilities. A treatment of a theory with an instability by expansion around the naive vacuum where the expectation values of all fields of Eq. (1) are 0 will then be incorrect. Excitations should instead be considered from the true vacuum, which breaks some of the symmetries of the action. We do not attempt to describe the true ground state but instead study expectation values obtained by expansion around the naive vacuum. Although these expectation values can not be trusted if an instability is present, they can show inconsistencies, like the above unbounded growth, that indicate the presence of an instability. Specifically, we study how the correlation functions we calculated in the previous section behaves at large spatial separations. The correlators are expected to decay at a certain rate with separation and when this is not the case, we attribute this to an instability.
With that as motivation, we study the behaviors of these exponents in the long distance limit in this section. We do that at the QCP and away from it, both at T = 0 and T > 0. We start out by only considering fermions coupled to one order parameter fluctuation field so we drop the index a for now.
Before calculating the h ± functions in different limits we give a physical interpretation of them and find some of their general properties. Looking back at where these functions show up we see that h + captures processes where the fermion exchanges bosons with fermions within the same patch of the Fermi surface. The fermion self-energy is obtained from this type of process so quasi-particle smearing will be governed by the behavior of h + . h − captures processes where a fermion exchanges bosons with a fermion in an opposite patch of the Fermi surface. This can give rise to the attractive glue needed for e.g. pairing.
Using that D(ω, k x , k y ) is positive and even in energy and momentum we have where we used the triangle inequality to move the absolute value inside the integral of Eq. (39) for the last inequality. This last identity will turn out to be important, it means that the interactions between opposing patches in a sense is stronger than the effects of quasi-particle smearing.

The quantum critical point
We now consider the theory at the quantum critical point where T = r = 0. We can write the correlation functions as where We diagonalize A, B to find out whether they are positive-/negative-or indefinite. The eigenvalues of the above matrices are both are given by By using Eq. (46) we find that The exponents are thus indefinite quadratic forms of the couplings, regardless of the form of D(ω, k). This means that it is, for both correlators, possible to find couplings λ ± σ such that the exponent is positive, regardless of the form of D(ω, k). For a scale-free D(ω, k) at the critical point we can find the dependence of the separation by a simple scaling argument. Instead of considering the uncorrected Green's function of Eq. (40), we might as well consider a generalization with critical exponent η: By rescaling the momentum integrals we have that We find that there are couplings such that the exponents of Eq. (37) and Eq. (38) grow unboundedly in the separation for all critical exponents η < 2. This leads to unphysical correlation functions that similarly grow indefinitely in the separation. In the case of the uncorrected Green's function for φ in our action we have η = 1. This unbounded growth of correlations is interpreted as an instability of the theory. We discuss this in more detail in subsection 3.3 and in section 5 where we consider a type of QCP at a time. Now we investigate for what couplings we find instabilities. Considering only a spatial separation, we have Note that R is independent of x at the QCP. We have not yet considered the oscillating part of the pair-pair correlation function. To do this we need to calculate h + σ to find its sign. We do this for a spatial separation x and find that it is positive everywhere, see Fig. 2 (result presented in Appendix A). Since it is positive we find that the last two terms of Eq. (48) decay for all couplings. We only considered spatial separations here but for τ = 0 it is possible to obtain an h + (τ, x) that grows in τ . However, this is harder to interpret since τ is Euclidean time. It seems likely that when continued to real time this once again is negative but that is something left for the future to be studied.

Zero temperature -away from criticality
We now continue to work at T = 0 but away from the critical point. We consider D(ω, k) of Eq. (40) with a finite gap r > 0 and once again calculate the above considered correlation functions at large spatial separations to see whether they show signs of instabilities. We have additional symmetry at the point where c = v F (see Chapter 5 of [33]) and we can calculate h ± explicitly there. We find where Ei is the exponential integral and γ is Euler's constant. Although not immediately evident from these expressions, for small x this behaves as in the critical case and h ± are linear in |x|. This linear growth slows down at a separation |x| ∼ v F / √ r and in the case of h + it approaches a constant while for h − it crosses over to a logarithmic growth. These asymptotic behaviors have been calculated for a general v F /c (see Appendix B): The h − function will thus dominate the exponents of Eq. (37) and (38) for large separations. This gives corrections to the power-law of the free theory, x −3 , at long distances: The prefactors C i are given by the asymptotic value of h + and the finite part of h − and depend on λ ±,2 σ , v F , c, r. Note that these expressions apply to any gapped boson propagator under the constraints mentioned in Appendix B. Interestingly, the corrected exponents may become larger than 1/2 meaning that the correlation functions grow unboundedly for large separations, just as at the QCP. The growth is now a power-law with an exponent that approaches ∞ as we approach the QCP where r = 0. The conditions for unbounded growth are: These are not the same conditions as what we had at the QCP in the limit r → 0, see Figure  3. The difference can be understood as the |x| → ∞ and r → 0 limits not commuting. The h + function becomes important when taking r → 0 first and affects the large |x| behavior. However, it is neglected when considering |x| → ∞ at any finite r since h + approaches a constant (that diverges as r → 0). The correlation functions do not have to grow unboundedly to indicate instabilities. The static pair and density susceptibilities can be calculated through a Fourier transform of the real space correlation functions where we neglect small separations and only consider x > Λ, Λ 1/k F . The static susceptibilities give the linear respons to sources for a pair or density operator. A diverging static susceptibility indicates a finite response without a source and thus an instability. First let us consider the T = 0 static susceptibilites of the free theory. Performing the τ and angular integral we have where J 0 is the zeroth Bessel function of the first kind. The density susceptibility is finite for all q whereas the pair susceptibility diverges at q = 0. There is an instability at T = 0 in the free theory leading to BCS superconductivity when we add any finite attractive four-Fermi interaction. Let us now see how interactions in the N f → 0 limit change this. We only Figure 3: The colored regions indicate couplings for which the corresponding correlator grows unboundedly in the separation at T = 0. The plots show this both at the QCP (blue) and for √ r = 0.1 (red). Note that the red regions move closer to the origin as r → 0 and thus cover the whole quadrants in the critical limit. Here we have used c = v F and thus R = 1/ √ 2.
consider the large x asymptotics of the integral to see whether it diverges in the IR so we can use Eq. (62)-(63). We find that the density susceptibility diverges at q = 2k F when This means that whenever λ − σ λ + σ < 0, there is a critical gap 0 < r < r c for the field φ where an instability is seen in the density correlator at wavevector Q = 2k F . Considering the pairsusceptibility, we find that the instability of the free theory is cured when both of these two conditions are satisfied: We have a faster divergence in the opposite case where the left hand sides are positive and the integral diverges algebraically instead of logarithmically in the IR.

Spontaneous symmetry breaking
In the previous subsections we found combinations of parameters where correlations of operators at different points grow unboundedly in their separation, and in a broader region, diverging static susceptibilities. We interpret the divergences as instabilities towards states where the corresponding operators have infinite expectation values. In practice, this is cut off by higher order terms that we have omitted in the action, but which become important for large values of the fields. The true ground states of the physical theories we are interested in are expected to instead be states where these operators have finite expectation values. We can not find the exact ground states without considering these higher order terms. The typical approach is to use mean-field theory to study the ground state. We leave that for future work and instead make the following natural guesses: where 0 < A, 0 < θ < 2π. The first is a pairing state that spontaneously breaks the U (N f ) symmetry whereas the second is a charge or spin order that breaks translation symmetry. The form of the density fluctuation is only symbolic, we expect oscillations around the momentum ∼ 2k F but how this happens in practice depends on the angular dependence of the coupling and what higher order terms stop the infinite growth of the correlator. We consider the expected ground states in more detail for specific systems in Section 5.

Finite temperature
The unbounded growths of correlations at T = 0 are interpreted as instabilities towards symmetry-breaking ground state. However, spontaneous symmetry breaking of continuous symmetries is not allowed at finite temperature in two dimensions for short range interactions by the Hohenberg-Mermin-Wagner [34] theorem. The interactions mediated by φ are short range for any finite r so we do not expect correlations to grow in the separation at finite temperatures away from r = 0. Fluctuations in the parameter θ in Eq. (71) and Eq. (72) are the Goldstone modes and these proliferate at finite temperature destroying any long-range order. However, it could still be possible to find quasi-long-range order such as in the twodimensional XY model below the Berezinskii-Kosterlitz-Thouless transition [35]. A phase with quasi-long-range order has correlators that decay with a power-law of the separation at finite temperature. Now we consider the long distance behavior of the pair and density correlation functions at a finite temperature T to verify lack of spontaneous symmetry breaking and to see whether they exhibit quasi-long-range order. The h ± functions are calculated through a Matsubara sum at finite temperature. It is useful to consider the zeroth Matsubara frequency (n = 0) and the rest of the sum (n = 0) separately. We define: Note that the functions h ± n=0 and h ± n =0 individually satisfy the properties in Eq. (42) to (46) and that Additionally, using we find that Figure 4: The single diagram contributing to the fermion self-energy at second order in the coupling.
The n = 0 contributions will thus not be very important since they are finite. We can express h ± n=0 (x) in terms of a Meijer G-function, see Appendix C. They diverge logarithmically for small r: This is an IR divergence. The same divergence shows up at each order in perturbation theory and can be seen by calculating the fermion self-energy perturbatively (see Figure 4): The ω m = 0 contribution diverges logarithmically as r → 0. Despite these IR divergences at each order, once all perturbative corrections are summed up into the exponentials of Eq. (37) and (38), we find that the correlation functions do not diverge in the r → 0 limit. This can be seen by consider the matrices A, B as in Eq. (47) and (48). Using Eq. (74) and Eq. (76) we can write Now both A and B are negative semidefinite (in the limit r → 0) so, regardless of the couplings, the IR divergence does not make these correlators diverge but may suppresses components of them to 0. We have thus found that correlation functions diverge as we take r → 0 at any finite order in perturbation theory. However, as we sum up all diagrams in the N f → 0 limit, we can safely take r → 0 and the considered correlation functions remain finite. This is similar to what was found in [30]. They study the same model (albeit at matrix large N ) but in d = 3 − dimensions. Similarly they find an IR divergence in the form of a diverging decay rate of electrons, in a perturbative treatment. The matrix large N limit allows them to sum contributions at all orders and they find that this cures the IR divergence and they obtain a non-diverging fermion self energy. A recent work has considered the same limit in d = 2 [36]. The IR divergence is then more severe and they find that it is not cured by summation of rainbow diagrams. Instead they argue that a φ 4 interaction that is irrelevant at T = 0 becomes relevant at finite T and there contributes a mass to the boson that cures the IR divergences. It is interesting that summation of rainbow diagrams did not cure the IR divergence while here the N f → 0 sum of diagrams cured it without any corrections to the boson.
Let us now depart from the critical, but finite temperature, theory and consider our correlators at a finite r and T and a large separation (x c/ √ r and x v F /T ). The denominators in Eq. (37) and (38) are given by hyperbolic functions at finite T . This gives a decay of exp(−2πT x/v F ) at large separations but interactions may give corrections to this. The h ± functions behave as for large |x|. We now have These are once again negative semidefinite and proportional to x. This means that interactions can only contribute to a shorter decay length at finite T except for the exceptional cases of interactions in the null spaces of A and B where it is unchanged. The long distance behavior of these correlators does thus always decay exponentially in the spatial separation. The static susceptibilities are then also finite since the τ integral is over a finite interval at finite T . No instabilities are present at finite T and since correlators decay exponentially in the separation, they also do not show quasi-long-range order. Let us consider the static susceptibilities at finite r in the T → 0 limit. The h ± functions behave as at T = 0 for separations x min(v F , c)/T, τ 1/T since the temperature only affects the boundary conditions of the theory (however, we show more rigorously that the asymptotic behavior is the same as at T = 0 up to |x| ∼ v F /T in Appendix C). The correlators are exponentially suppressed for large spatial separations as exp(−2πxT /v F ) or faster and the time direction is periodic with periodicity 1/T . This means that the finite temperature effectively introduces cut-offs in the Fourier transforms of the T = 0 static susceptibilities. We consider an explicit cut-off at x 2 /v 2 F + τ 2 < Λ ≡ /T where is small such that the T = 0 correlators in Eq. (62) and Eq. (63) can be used up to this cut-off. The Fourier transform is largest at momenta where oscillations cancel, this happens at q = 2k F in the density susceptibility case and for q = 0 in the pair susceptibility case meaning that any potential divergences show up first at these momenta. Performing the cut-off Fourier transform at these momenta we find: By finite we mean terms that are bounded as T → 0. We expect the contribution from outside the cut-off to behave the same way since the integral is bounded in the τ direction  Figure 5: These figures show the static charge (or magnetic) susceptibility against momentum for r in the region with an instability at T = 0 (left) and outside this region (right). The susceptibility has been calculated by a Fourier transform of (37) where the h ± functions have been obtained numerically. Here we consider the somewhat artificial case of λ + σ = −λ − σ ≡ λ in all directions such that we have a rotationally symmetric density correlator. We further use k F = 5λ 2 /v 3 F and c = v F . We see that the susceptibility diverges upon lowering temperature for r < r c and approaches a constant for r > r c . and exponentially decays in the spatial directions. The full susceptibilities then behave like this: for low temperatures. The susceptibilities can easily be calculated numerically as well, see Figure 5 for the charge/spin susceptibility at different momenta.

Landau-damping corrections
The N f → 0 limit removes the corrections from the fermions onto the field φ such that it is still Gaussian and can be integrated out. However, we can add some corrections to the φ two-point function without making it non-Gaussian. These are important corrections to the IR at finite N f so being able to add these is interesting because it gives us an indication as to how N f corrections modify the results away from N f = 0. This was done systematically in [29] by considering a double limit of N f → 0 and k F → ∞ keeping N f k F constant. This limit was shown (with some caveats concerning the order of other limits) to suppress all symmetrized Figure 6: A two-vertex fermion loop giving corrections to the field φ.
fermionic loops with more than two vertices, order by order in perturbation theory, and thus allowing us to still integrate out the order parameter fluctuations. Fermion loops are also suppressed in the matrix large N limit studied in several works mentioned in the introduction. A similar double limit reintroducing some of the effects from fermion loops was studied in [19].
We would like to consider the same N f → 0, k F → ∞, N f k F = constant limit as in [29] but for the systems and correlation functions we are studying here. In this work we allow the coupling function to depend on momentum and we consider finite temperatures so we need to verify that the same suppression of fermionic loops with more than two vertices remains. This can actually already been seen concerning corrections to the φ two-point function from our above calculated density correlator. The field φ couples to the fermion densities so the amputated φφ correlator is given by ρρ . From Eq. (37) we see that ρρ is composed of a term with small momenta (compared to k F ) and a term with momenta of order k F . Only the former term is relevant since the low-energy φ modes have Q = 0. This term is independent of the couplings and is thus simply what is obtained from the 2-vertex fermion loop. To verify suppression of higher n-point φ corrections in this limit we explicitly calculate such fermion loops in Appendix D.
Having verified this cancellation we may consider what this limit entails for the correlation functions we are interested in. The earlier strict N f → 0 limit lets us use the free boson correlation function when evaluating the h ± -functions. Now we should calculate the h ± functions with a φ two-point function that is corrected be 2-vertex fermion loops, see Fig. 6. The N f → 0 limit was also used to suppress diagrams of the type in Fig. 1(b) contributing to the density-density correlator. Those diagrams come with an extra factor of k F as opposed to Fig. 1(a) and can seemingly not be neglected in the double limit we now consider 3 . The diagrams of the type in Fig. 1(b) can not simply be found by using the background field Green's function since we can only calculate it for insertions separated a distance x 1/k F and here we would need to have two insertions at the same point. We thus need to evaluate these diagrams some other way. Luckily, this class of diagrams turns out to be manageable. Let us first consider such diagrams in momentum space, with a large momentum incoming to the density vertex q k F . The momenta going through the bosons connecting the two fermion loops are now of order 4 k F . The propagators are of the form 1/(ω 2 + c 2 q 2 + r 2 − Π) so these diagrams get a further k F suppression and need not be considered.
Next consider the case of a small momentum incoming to the density vertex, q k F . In the above we found that all loops with more than 2 vertices cancel to leading order in large Figure 7: This figure shows a class of bubble chain diagrams that contribute in the combined limit N f → 0, k F → ∞ that is not captured by the above framework and has to be calculated separately. The ellipsis indicate that this class of diagrams has any number of bubbles, larger than 1.
k F . This assumed all incoming momenta were small compared to k F and it is now applicable. We thus have that the only diagrams contributing to leading order are the ones with a single boson connecting the left and right fermion loop in Fig. 1(b) such that both loops each have two vertices. The boson being exchanged however receives corrections from further 2-vertex fermion loops so when q k F we must sum all diagrams of the form shown in Fig. 7. The first and last vertices are momentum independent for the density-density correlator. All other vertices come with a momentum dependent coupling constant λ(θ). We write the two-vertex fermion loop as Π 2 [f ](ω, q) where f σ (θ) is the combined momentum and spin dependence of the two vertices. We have that the summed contributions from diagrams of the type in Fig. 7 is given by We now need to calculate the two-vertex fermion loop Π 2 [f ](ω, q) both for this additional contribution to the density correlation function and to be able to calculate the h ± functions. Π 2 [f ](ω, q) was calculated at a finite temperature in d = 3 in [30]. They found that the result is identical to the zero temperature result and this is true also in d = 2. Only fermions close to the Fermi surface contribute for φ momenta small compared to k F and we can approximate (see Appendix E) this loop by the integral where θ is an angle measured from the external momentum q and f σ (θ) is the momentum dependence of the vertices evaluated at the Fermi surface: f σ (θ) = f σ (k Fn (θ)). To solve this we expand the vertex momentum dependence in Fourier series Calculating the sums and integrals (see Appendix E) we havẽ Now we consider the boson self-energy and replace f σ (θ) = λ 2 σ (θ). The self-energy dominates over the ω 2 in the (critical) bosonic propagator and the low-energy scaling is The fraction in the power can then be replaced by −i sgn(ω n ) in this limit. We then find Assuming λ σ (θ) is parity symmetric or antisymmetric we find exactly the result of [13]: to leading order the boson only couples to fermions with perpendicular momenta. Interestingly, for a non-centrosymmetric theory where λ σ (k) 2 = λ σ (−k) 2 , we find that this is not the case and there is an extra contribution from the boson coupling to fermions all around the Fermi surface given by the principal value integral in Eq. (94). This extra contribution is odd in both frequency and momentum (the angle θ was defined w.r.t. the incoming momentum q).
Analytically continuing this to Minkowski signature we find that it is not a damping term but instead it contributes a singular real part to the polarization. We are not aware if this extra term appearing in non-centrosymmetric systems has been studied before but it is outside the scope of this paper so it is left for future studies. We only consider centrosymmetric systems henceforth. Now that we have calculated the boson self-energy we first use this to calculate the ρρ(q) bubbles contribution (see Eq. (88)) to the density-density correlation function. Eq. (88) is written in momentum space but we are interested in correlators in real space at large separations and thus Fourier transform this. We stay away from criticality and consider distances longer than c/ √ r. We do this in Appendix F (for a constant λ(θ) = λ) and find: for large x and τ and this contribution does thus not affect any conclusions regarding instabilities so we need not consider this more. Now we use a corrected D when calculating h ± LD : We use the subscript LD to indicate that D and h ± have been calculated with the above Landau-damping correction. This calculation is performed in Appendix G. We calculate the long distance behavior of h ± LD (τ = 0, x) to see if the instabilities found in the previous section have changed. Away from criticality (r > 0), at both T = 0 and T > 0 we find that only the finite part of h ± LD (τ = 0, x) has changed compared to h ± (τ = 0, x). The growth at long distances is the same and the conclusions about instabilities are thus unchanged by the inclusion of resummed 1-loop Landau-damping corrections. However, right at the QCP where r = T = 0 we find that the h ± LD (τ = 0, x) functions no longer grow linearly in the separation but now grow with a power-law: Note that there are now two scales at the QCP: λ 2 and M D so the h ± LD (τ = 0, x) functions only have this behavior asymptotically. Since the h ± LD functions both have the same powerlaw growth we have to consider them both when considering instabilities. Analogously to the non-Landau damped case we find that unbounded growth of correllators is governed by Eq. (56) and R is now independent of all parameters and given by R = 1/ √ 2.

Results for specific systems
In the following subsections we use the above results to study specific fermionic systems near QCPs to obtain their N f → 0 phase diagrams. We start each subsection with a brief description of the quantum phase transition and how its fluctuations couple to the fermions, i.e. what the structure of λ σ (θ) is. We finally also consider another system that can be described by the same model, fermions coupled to an emergent gauge field. While the study of these models is motivated by experimental systems, it is not very useful to make a direct comparison with our results since we consider very simplified models of them. Instead we focus on some earlier theoretical results on the same and similar models when making comparisons.

Charge nematic transition
In this subsection, we consider an l = 2 charge nematic quantum critical point such as the ones that have been observed in the Cu- [1,2] and Fe- [3,4,5] based SCs. These QCPs are due to electronic instabilities towards a state where the Fermi surface shape breaks the C 4 rotational symmetry of the lattice. The resulting nematic state has a remaining C 2 rotational symmetry and symmetry under spin inversion. The order parameter is a single real scalar field that is odd under π/2 rotations. This means that λ σ (k) is independent of σ and has d-wave symmetry. In the analysis of the previous sections we considered correlators of fermion bilinears separated spatially in a chosen directionn. For a particularn we then have λ + ↑ = λ + ↓ = λ − ↑ = λ − ↓ ≡ λn where λn has d-wave symmetry and thus is 0 in four different directions, the so-called cold spots. For T = 0 and r > 0, we find from Eq. (62) that nematic fluctuations suppress the charge and spin correlators at long distance: This is identical to what happens near a ferromagnetic QCP so we defer the discussion of this to the next subsection. From Eq. (87) we find that the nematic fluctuations lead to enhanced pairing fluctuations as T = 0 is approached, most strongly near the QCP: Instead of a logarithmic growth of fluctuations, now fluctuations grow with a power-law as T → 0 and furthermore, the exponent diverges as the QCP is approached. As noted in the previous section, these results are unchanged by the inclusion of Landau-damping corrections. While we have considered s-wave pairs, this result does not necessarily mean that the ground state is an s-wave SC. We could similarly consider e.g. p and d-wave pairs that contain spatial derivatives and different spin structures. This can be investigated within the framework we have used here but due to spatial derivatives breaking rotational symmetry it will be a longer calculation that we leave for the future. Instead, we make an educated guess of the result. The spin is unimportant in this case since λ ↑ (k) = λ ↓ (k). The conclusions regarding instabilities come from the asymptotic behavior of the correlation functions. We have derivatives acting on the background field Green's function of Eq. (36) for non-s-wave pair correlators. The derivative may hit any of the factors in Eq. (32). Hitting the exponent −isk F |x 12 | only brings down a constant and the resulting term does thus have the same asymptotics as for the s-wave case. While it is possible some cancellation occurs for e.g. pwave, this is not expected to be universal and happen for states of all symmetries apart from s-wave. This indicates that whether the pairing instability we found here results in s-wave or a gap of some other symmetry is something we cannot tell using this analysis. It seems likely several of these correlators diverge at the same points in the phase diagram and the groundstate symmetry is determined by the next to leading order N f and 1/k F corrections and higher order terms in the action. Since the coupling λn has a d-wave form factor and thus is 0 on four nodes at the Fermi surface, there will be directions in which the instabilities are not seen. This allows for nodes or lukewarm regions along the Fermi surface where the gap is 0 or small and thus compatible with a d-wave state.
We have found that nematic fluctuations lead to enhanced superconductivity in the N f → 0 limit but the symmetry of the gap is not definite. This is consistent with studies of the charge nematic transition in two dimensions in the matrix large-N limit [21] where nematic fluctuations enhance SC but also there the symmetry of the SC groundstate is non-universal. The matrix large-N limit has also been employed on the charge nematic transition in 3 − dimensions [20,22] where is small. There the order of the → 0 and N → ∞ limits decide whether the system is a SC or not [22]. [37] studies the same model of a nematic QCP (at N f = 1) as us and additionally include four-fermion interactions. By considering small interactions and staying a small distance away from the QCP they can integrate out the nematic fluctuations and get a correction to the four-fermion interaction, essentially staying within a Fermi liquid framework. They solve the gap equation and find that this correction gives an increase in T c compared to the BCS case as criticality is approached. Like us, they find that nematic fluctuations increase susceptibility towards pairing of both s-wave and dwave symmetry and find that the bare four-fermion interaction is what decides the symmetry of the resulting superconducting state and not the near-critical nematic fluctuations.
The authors of [38] include one-loop self-energy corrections to both the boson and the fermions and argue that this allows them to access the strongly coupled regime at the QCP. The order parameter field is static in their case but they are otherwise in the same limit as we have assumed; they consider l = λ 2 /v F k F a small parameter. They find that nematic fluctuations lead to superconductivity, with a finite T c , and they also find that the SC symmetry is non-universal.
Finally we mention that the two-dimensional charge nematic transition has been studied using sign-problem-free determinant quantum Monte-Carlo (DQMC). This was done by putting fermions on a finite lattice and coupling them to pseudospin-1/2 degrees of freedom [39,40]. The spins have a nematic QCP whose fluctuations give strong corrections to the fermions. The first work does not find any pairing instability however they find large pairing fluctuations. In the second study the authors increase the coupling between the spins and the fermions beyond what is physical in a microscopic interpretation of the model and they instead regard it as an effective model. They then find a superconducting dome covering the QCP with s-wave symmetry in the symmetric phase.
As opposed to these earlier studies, we only find a pairing instability at T = 0. However, we believe this is expected due to the Mermin-Wagner theorem. By departing from the strictly two-dimensional case and including four-fermion interactions at finite N f we expect an increase in T c compared to that of BCS theory since we find an increase in pairing susceptibility due to nematic fluctuations. This is something we leave for future works to study.

Ferromagnetic and spin nematic transition
Ferromagnetic and spin nematic ordering both break spin-inversion symmetry. Electrons spontaneously order their spins in the same direction in a ferromagnetic phase.
The spin-nematic phase is similar to the charge-nematic case, the Fermi surface spontaneously breaks rotational symmetry, but in perpendicular directions for spin up and down. In fact, all of the orders considered here can be seen as coming from Pomeranchuk instabilities. Ferromagnetic ordering is simply the Fermi surface of one spin spontaneously becoming larger than for the opposite spin.
There will be soft fluctuations in the order parameter in the case where these transitions are continuous. The fluctuations couple with opposite sign to fermions of opposite spin, λ ↑ (k) = −λ ↓ (k), for both the ferromagnetic and the spin-nematic transition. However, λ ↑ (k) has s-wave symmetry in the former case and d-wave symmetry in the latter. For a chosen spatial directionn we then have λ + ↑ = λ − ↑ = −λ + ↓ = −λ − ↓ ≡ λn for the coupling of both ferromagnetic and spin-nematic fluctuations. Their effects on the fermions will thus be the same when studied using the framework of this paper, with the caveat that in the spin-nematic case there are four coldspots in λn and none in the ferromagnetic case.
Plugging in the coupling function in Eq. (56), we do not find any instabilities at the QCP. In fact, we find that components of the leading small N f CDW/SDW fluctuations and pairing fluctuations receive faster power-law fall-offs at T = 0 due to interactions: This means that the pair susceptibility is now finite at T = 0 (see Eq. (87)) and the divergence in the free theory is cured by interactions with the ferromagnetic or spin-nematic fluctuations.
We thus have a naked QCP with strong NFL behavior but suppressed pairing and suppressed Friedel oscillations. At finite temperatures we find that the correlations lengths ln of the above components receiving corrected power-laws also decrease due to interactions (see Eq. 37-38): We have been considering the leading contribution to the correlators at x 1/k F . With these corrections due to interactions we find a faster decay in separation x for some of the terms above so it is possible that the subleading behaviors that we have not considered are actually dominating these. To be perfectly consistent we should then remove these terms above and regard them as part of the neglected terms that are subleading at long distances. The leading long-range behavior does thus not necesarily behave as described by Eq. (101)-(103) at large separations. The above asymptotic result should be taken with caution and conservatively viewed as simply a cancellation of the diagrams at leading order for large x.
Fermions in 2D coupled to ferromagnetic fluctuations were studied in [41] using DQMC. The authors consider a system with two flavors (orbitals) of spinful fermions. They find spin triplet SC tendencies but at much lower temperatures than where NFL effects set in. This is opposed to what is found in the charge nematic case [40] where superconductivity sets in at a similar temperature as NFL effects. By adding an orbital index to the fermions we may also consider orbit-singlet spin-triplet pairing instead of the spin-singlet pairing studied so far in this work. We define the spin-triplet pair operators We can calculate the pair correlator the same way as in the case of the spin-singlet pair, the only difference is that we never encounter λ σ (k) with mixed spins, the result is the same as what is found for the spin-singlet pair correlator near the nematic QCP studied in the previous subsection. This means that in a two-orbital system coupled to ferromagnetic fluctuations, we also find a pairing instability in the spin-triplet channel in the N f → 0 limit. However, this behaves similarly to in the charge-nematic case whereas [41] found that there is a qualitative difference between ferromagnetically induced spin-triplet SC and charge-nematic induced spin-singlet SC.

Circulating current transition
The circulating current phase proposed by Simon and Varma [42] breaks lattice C 4 rotational symmetry and time-reversal symmetry. It requires a two component order parameter φ a = (φ x , φ y ) that couples to the fermions through a coupling function with p-wave symmetry independently of spin: λ a,σ (k) ∝ k a [13]. As we only consider two antipodal patches in the directionn at a time we may simply consider a single φ corresponding to φ a projected in thê n direction. We then have λ + ↑ = λ + ↓ = −λ − ↑ = −λ − ↓ ≡ λn and no cold-spots despite the pwave symmetry since we actually have a two-component order parameter. Again considering Eq. (56), we find that fluctuations at the circulating current QCP lead to an instability. This time in the charge/spin density channel at momentum Q = 2k f . Departing from the QCP but still at T = 0, we find that the divergence in the charge and spin susceptibilities persists up to a critical order parameter fluctuation gap r c , see Eq. (68): Density fluctuations diverge upon approaching the T = 0, 0 < r < r c interval either from finite temperatures or larger gaps. A phase diagram is shown in Fig. 8. The interval terminates at a new QCP where there now are large fluctuations at wave-vector Q = 2k F . These fluctuations may give rise to important finite N f corrections, something we leave for the future to consider. As in the above ferromagnetic case, we find that s-wave pairing is suppressed and there is no superconducting instability at T = 0. In the ferromagnetic case this was due to the coupling acting with different signs on the electrons with opposite spins making up a Cooper pair and we saw that it is still possible to have triplet pairing in the ferrromagnetic case. In this case the coupling is independent of spin, but it couples with opposite sign to fermions on opposite sides of the Fermi surface and the coupling does not induce triplet pairing in this case.

Spin liquids: Emergent gauge field
As mentioned in the introduction, Eq. (1) can also describe fermions coupled to a U (1) gauge field A which appears in models of spin liquids. The electric potential A 0 is screened and can be neglected [12]. By defining φ a = A a we can capture the interaction of the magnetic potential and the fermions through our action Eq. (1) with the following choice of coupling function where e is the coupling strength. The longitudinal part of A a is set to 0 in the Coulomb gauge (k · A = 0). Once we pick a directionn and only consider patches k = ±k Fn we may thus neglect An and simply call the transverse component φ. We thus find that we have a similar coupling function as in the case of a circulating current phase transition: Also in this case, there are no cold-spots since we picked the directionn arbitrarily. At T = 0 we thus find that this system also has an instability towards forming charge/spin density waves around Q = 2k F . Now r = 0 since φ is massless by gauge invariance and this instability is present for arbitrarily weak but non-zero interaction strength e. This result is opposed to what is found in the (vector) large N f limit [12] where no antiferromagnetic instability is found at T = 0. This model at T = 0 has also been studied at matrix large N in [16]. The authors find that an instability at 2k F in the charge density is present at large N , but suspect that the point of interest N = 2 may be stable. They proceed by considering a modified boson dynamical critical exponent z b = 2 + (in our model z b = 3) and a combined limit where both N → ∞ and → 0. N tunes the relative strength of the quasiparticle smearing that suppresses instabilities and what they call "Amperean attraction" which enhances it. Taking these limits, they find that whether the theory is unstable depends on the order of the two limits. Whether the physically interesting case of N = 2, = 1 is unstable then depends on what the -N phase diagram looks like. Now we consider finite temperatures. The apparent IR divergence discussed in the paragraph after Eq. (80) becomes relevant for this system since r = 0. We find that the N f → 0 spin liquid has the leading large distance pair correlator completely suppressed at finite temperatures because of the divergence in the h ± functions and it is entirely composed of the subleading terms we have not considered in this work. Since the divergence is due to the n = 0 Matsubara sum term, it is unaffected by the inclusion of Landau-damping.
The n = 0 terms in the h ± functions cancel out in the charge (and spin) density correlator exponent and it decays exponentially at long distances for all finite temperatures. Its precise form is found by calculating the finite piece of the h ± functions. However, since the correlator grows exponentially in separation at T = 0 (or as ρρ ∼ exp(Ax 1/3 ) when including Landaudamping, see Eq. (97), (98)), we expect the charge and spin susceptibilities to diverge faster than a power-law as T = 0 is approached and the τ, x cut-offs move to infinity as v F /T , see end of Section 3.4.

Multiple critical points
We can easily consider a system close to multiple QCPs of different natures. Several types of order parameter fluctuations give corrections to the fermions that have to be treated simultaneously. As we see in Eq. (37) and Eq. (38), these corrections simply sum in the exponents. We have seen interaction effects both suppress and enhance the long distance correlation functions so it is possible that the instability induced at T = 0 near one QCP is cured by the vicinity to another. Whether an instability shows up or not depends on the relative strengths of the interactions in the long distance limit. This depends on the coupling strengths that are direction dependent so we consider a directionn at a time. We find that the relevant strength of a QCP is of the form where a refers to the index on the corresponding order parameter field φ a and λ a,n are the relevant components of the interaction strength as defined in the above subsections. These strengths are unaffected by the inclusion of Landau-dampng corrections away from the QCP.
We find the following conditions for CDW/SDW and pairing instabilities: where CN, FM, SN, CC refers to charge nematic, ferromagnetic, spin nematic and circulating current transitions, respectively. The effects of the different fluctuations are as expected from the above subsections. α −1 a can be viewed as the effective distance to the QCP labeled a. A three-dimensional phase diagram is shown in Fig. 9 with distances to charge nematic, ferrromagnet and circulating current transitions on the axes.
Noting that the strengths α a , as defined, are always non-negative, we find that there is no phase competition between the pairing state and the CDW/SDW states, only one of the conditions above are satisfied at a time. However, we have only considered a single direction n. The nematic fluctuations have cold spots and this means that we may have systems with phase competition. I.e. they may show a CDW/SDW type instability when considering one part of the Fermi surface and a pairing instability when considering another part and we cannot resolve within the current framework which one is winning.
Phase non-coexistence/competition of antiferromagnetism (AF) and SC is also seen in unconventional SCs with AF T c going to 0 at the point where SC T c is optimal [5,43,44]. While this may be indicative of AF fluctuations near an AF QCP enhancing superconductivity, we find that another possible mechanism giving a similar result may be seen for an AF phase in the vicinity of a nematic QCP. Pairing is most strongly enhanced right at the nematic QCP and at the same point we find that the exponent in the power-law decay of the normal-state AF correlator goes to −∞. It is not entirely clear how the interactions giving rise to the AF phase would interact with the nematic fluctuations. However, if those interactions can be written in the form of gapped CC fluctuations, then we find from the above that they are completely outdone by the fluctuations at a nematic QCP and the AF state should disappear at or before the nematic QCP where α CN = ∞.

Discussion
In this paper we have investigated instabilities of fermions in two dimensions coupled to order parameter fluctuations near different QCPs and emergent gauge fields, in the N f → 0 limit. Taking this limit allowed us to calculate charge, spin and pair correlation functions at long distances, x 1/k F . First we considered a general model without specializing to a particular system. While we were agnostic to the momentum dependence of the coupling functions, we assumed an otherwise rotationally symmetric Fermi surface for convenience.
We found expressions for correlators in terms of functions we label h ± . They are responsible for interactions between fermions within a patch, and for interactions of fermions in opposite patches, respectively. This means that h + is solely responsible for quasi-particle smearing effects while the attractive glue between fermions giving rise to pairing, CDW or SDW is generated by h − . We find an inequality, Eq. (46), that settles that the former effect is stronger.
We considered both gapped and critical fluctuations and zero and finite temperatures. We found that an IR divergence in perturbation theory at finite temperatures for critical fluctuations was cured by summing all diagrams in the small N f limit. All considered correlation functions show exponential decays at large separations at finite T and power-law behavior at T = 0, for finite gaps. We found that there are cases where the powers are negative such that correlation functions show unphysical asymptotic growths at long distances, but also milder cases with power-law decay where nonetheless static susceptibilities diverge. We interpreted this as instabilities towards states where the corresponding operator condenses at T = 0. We proceeded by calculating the behavior of susceptibilities as T = 0 is approached from finite T .
Next we considered a modified boson propagator that includes the one-loop Landaudamping correction. We show that also in this more general setup, as compared to [29], we can interpret this as the same controlled approximation in the double limit N f → 0, N f k F = const. We find the inclusion of fermion loops with two vertices does not change the asymptotic correlator decay lengths at finite T nor the exponent in the power-law at T = 0. However, the long distance behavior right at the QCP was shown to be affected by the inclusion of Landau-damping corrections and the constant prefactor and subleading terms are affected by Landau damping.
The diverging susceptibilities are all proportional to N f which is taken to 0 here so one may question whether the divergences indicate instabilities since we find them in the limit N f → 0. We argue that for the purpose of studying the N f → 0 limit to get insights about the finite N f case, these divergences should be taken seriously. Instead of defining the operators ρ α and b as sums over fermionic bilinear, we could have defined them as ψ † 0,α ψ 0,α and ψ 0,↑ ψ 0,↓ , respectively, instead. Perhaps a bit awkward, but we could then consider susceptibilities of such an operator as a function of N f and take N f → 0 to find that at least there, the susceptibility diverges. This corresponds to the susceptibility in the N f = 1 case diverging, when only summing diagrams without internal fermionic loops. We found that after adding all diagrams with fermionic loops with only 2 vertices to this, the divergence of the susceptibilities around (but not at) the QCP are unchanged. If that is any indication as to what happens when adding further finite N f corrections in terms of diagrams with fermionic loops with more than 2 vertices, we would expect to find the same divergences of susceptibilities at finite N f . However, this indication may be deceiving and these further corrections may cancel out or dominate the divergence found here.
Finally, in Section 5 we applied our findings to specific systems of fermions near QCPs to see whether fluctuations induce instabilities in the N f → 0 limit. We considered charge nematic, ferromagnetic, spin nematic and circulating current transitions, systems close to multiple of these transitions and finally we also considered fermions coupled to an emergent gauge field.
One of the main reasons for studying this kind of model is the unexplained SC dome in the vicinity of nematic transitions in the cuprates and iron pnictides. As mentioned in the introduction, nematic fluctuations may both lead to a breakdown of quasiparticles that limit pairing and they act as an attractive force inducing pairing. Here we found that the latter effect wins and superconductivity is enhanced by nematic fluctuations in the small N f limit. Furthermore we argued that the gap symmetry is non-universal, it depends on terms neglected in the low-energy limit we are considering.
In the case of a continuous ferromagnetic or spin-nematic transition we found a naked metallic QCP in the case with one fermion orbital. This is interesting since it admits the study of strong NFL physics all the way to T = 0. The naked QCP is crucially single orbital and it is thus not amenable to DQMC relying on two-orbitals to cancel the fermion sign-problem. This makes an alternative approach like the small N f limit currently necessary to study this naked QCP. In the case of two orbitals we find a spin-triplet SC groundstate as also develops in DQMC as temperature is lowered [41]. The results near the rotation and timereversal symmetry breaking circulating current transition was maybe the most intriguing.
Here we find a phase diagram with a SDW/CDW state induced in a finite gap interval at T = 0 near the original QCP and suppressed spin-singlet SC.
We made some comparisons to earlier studies on the same systems. Most interesting are perhaps comparisons to DQMC results. While our results are largely consistent in what instabilities we find, it is hard to make sharp comparisons at this stage. The results found here are valid at distances 1/k F while the sizes of systems amenable to DQMC are still quite limited. Additionally, the instabilities we find are only at T = 0 which is inaccessible to DQMC. One could approach T = 0 and compare the growth of correlators but to do this we would need to match the couplings of our field theory to the interactions in the lattice models studied using DQMC. Admittedly, the N f → 0 limit might be too far from the finite N f case to warrant this type of comparison. However, DQMC provides an easy way to study the N f → 0 limit at finite temperature. The acceptance-rejection probability within DQMC depends on the ratio of fermion Green's function determinants and by raising it to the power of N f , one can easily account for different values of N f . That way it should be possible to make contact with the type of results we provide here and DQMC simulations. The benefit is two-fold, firstly it provides a non-trivial way to benchmark DQMC simulations with analytical results. While the analytic results are at N f = 0, it is still at a strongly interacting point not previously studied. Secondly, it gives us a way to measure the effect of tuning N f that can shed light on whether the phenomena found here are specific to N f = 0 or if they in some form extend to finite N f .
The fact that the k k F part of the ρρ correlator receives no corrections from interactions is an interesting result. This was also found in [45] and [46] due to cancellations upon symmetrizing fermion loops with density vertices. We have generalized that result by allowing for momentum dependent vertices and finite temperatures. This has prompted further investigations that will be presented in [47].
When adding Landau-damping corrections to the boson with momentum dependent coupling, we found that corrections from non-centrosymmmtric systems seem to disobey the earlier belief [12,13] that only two patches are important. In this case we in fact find that fermions receive relevant corrections from bosons that have relevant corrections from fermions in all patches. This also warrants further investigation.
A surprising find is the curing of the finite temperature IR divergence in perturbation theory since this was not found to be cured in the recent matrix large N treatment of the same model in [36]. It would be interesting to understand whether the additional crossed diagrams cured the divergence, or if there is some other difference between our works.
The framework developed in [24] and used here has proved to be quite convenient for studying quantum critical metals, albeit at N f → 0. Many results can be found analytically and the calculations are not more difficult than those in one-loop perturbation theory. In particular we find the calculations easier than the much employed matrix large N limit that contains a subset of the diagrams in the small N f limit. We believe it would be useful to consider this framework further to develop an understanding similar to that of perturbation theory. Specifically it would be interesting to systematically understand how properties of the boson Green's function affects the h ± functions and what combinations of h ± functions show up in the exponents. For all the coupling functions we considered, and for all the fermion bilinear correlators we considered, we did not find any quasi-long range order and we found that the finite temperature IR divergence was cured. While not manifest, it seems likely that these findings are more generally true, for all couplings and correlators. It would be interesting to investigate whether this is the case, and what happens for other boson Green's functions.
A h ± at the QCP Calculating h ± for τ = 0 at the QCP T = r = 0 we find for v F = c and for v F = c we have Despite not manifestly so, these functions are continuous in v F at v F = c.
B Calculating h ± with a gapped boson at zero temperature Consider a gapped boson propagator D that after integrating over k y is continuous at ω = 0, k x = 0: and further decays algebraically at large ω and k x . We want to calculate the large τ , x behavior of the corresponding h ± -functions at T = 0. The h ± integrals can be viewed as regularized, symmetrized Fourier transforms of the h ± integrands without the numerator, where is a cutoff that regularizes the divergence at ω = k x = 0 and The unregularized Fourier transform diverges for all τ , x due to the singularity at ω = k x = 0 but the divergence is uniform so the subtraction of the τ = x = 0 component allows us to remove the regulator . For any > 0, we have that d ± is absolute square integrable when using the same cutoff and so is its Fourier transform by the Plancherel theorem. This means that the Fourier transform necessarily decays at large τ , x for any finite . This means that if the h ± functions were to not approach 0 at large τ , x, then the behavior at infinity can be obtained from d ± (ω, k x ) for ω 2 + k 2 x < 2 for any > 0. We can thus use to calculate the large τ , x limit for a gapped D at T = 0. We then find: Here the finite part depends on how D(ω, k x , k y ) behaves away from ω = k x = 0.
C Calculating h ± at finite temperature We calculate the n = 0 contribution to the h ± (τ, x) functions at finite temperature for φ propagator We find where G is the Meijer G-function. In the large |x| limit we have and taking r → 0 we find Now we consider the double limit c/ √ r x v F /T . Using that the φ Green's function decays monotonically in |ω n | we can bound the h − Matsubara sum with T = 0 integrals on both sides shifted by either an extra or a removed h − n=0 to obtain Further we use that to find We thus find that

D Symmetrized loop cancellation
Here we show that to leading order in large k F , fermion loops with 3 or more vertices cancel upon symmetrizing the vertices. We do this at a finite temperature T and allow for momentum dependent vertices. We start off with two lemmas: Lemma D.1. For n ∈ N, and n ≥ 2 we have n j=1 n i=1,i =j Proof. Consider a contour C at ∞ that surrounds all poles: The contour integral is 0 for 2 ≤ n.
Definition D.1. We call a set of complex numbers Z nonexceptional if for every non-empty proper subset Z 0 ⊂ Z we have where w i , z i ∈ C, {z i } is nonexceptional and S n is the symmetric group over the first n positive integers.
See [48] for a proof of this identity. A simplified proof together with a generalization of the following loop cancellation is to be published [47]. Now consider the finite termperature symmetrized fermion loop with n vertices with external momenta q i = ( i , q i ): where ω m = (2m + 1)πT are fermionic Matsubara frequencies, k = (ω m , k), and By energy and momentum conservation we have Q σ n = 0. The function f is the momentum dependent coupling function and it varies between momenta of order k F . The momentum integral is dominated by the region close to the Fermi surface. We restrict the momentum integral to a patch P of the Fermi surface small compared to k F . The coupling function can then be treated as a constant within this patch, with the corrections being next-to-leading order in large k F . We then need to prove that the cancellation of the large k F leading order term happens after symmetrizing, but before summing up all different patches. We use patch coordinates: k x across the Fermi surface and k y parallel to it and we linearize the Green's function: First we integrate out k x with a contour around the upper half-plane.
where θ is the Heaviside step functions. Now the ω m sum is trivial. The summand is clearly 0 for ω m < min j (−Ω j ). For ω m > max j (−Ω σ j ) we find that the ω m -summand is of the form n j=1 n i=1,i =j and we can apply Lemma D.1 to find that it is zero there as well. We can thus restrict the ω m sum to some upper cutoff ω N > max j (−Ω σ j ). The benefit is that each of the terms in the j sum now converge and we can exchange the orders, perform the now restricted ω m sum to obtain Here we have used Lemma D.1 again to get rid of the upper limit of the primitive function. We note that this is now independent of temperature. Finally, we are interested in the large k F limit. We cannot exchange the limit with the integral right away since the k y that are relevant will grow with k F in the limit. We therefore make the change of variables k y → sk F and we have Now it is safe to take the large k F limit of the integrand and we obtain an expression that it is of the form of the sum in Lemma D.2. It thus cancels out to leading order in large k F upon symmetrizing. Note that this happens before performing the k y integral. This argument is made more rigorously and is further generalized in [47].

E Calculating 1-loop boson polarization
For brevity, we consider spin-and flavorless fermions here. The one-loop boson polarization is given by: Performing the Matsubara sum we have 2(iω n + (k) − (k + q)) .
The denominator vanishes at 2k = −q if ω n = 0 but the numerator also vanishes there so the above integrand is finite and bounded as long as the coupling function is. We go to polar coordinates k, θ where θ is the angle between k and q. The difference of the hyperbolic tangent functions decays exponentially away from the Fermi surface with decay length T /v F . We can thus linearize the dispersion, coupling function and integral measure around k = k F for q, T /v F k F and additionally we can extend the k integral to all of R: − tanh v F (k−k F +q cos(θ)) 2T 2(iω n − v F q cos(θ)) + subleading.
Performing the k integral we find We find that this is q independent and generally non-zero for ω n = 0. This means that interactions generate a boson gap but we have already considered a gap r so we can simply absorb this in a redefinition of r. We continue to work with the polarization with this constant subtracted and define:Π We now expand the function f in Fourier series according to Eq. (90) and discard the sin(θ) terms since they are odd under the θ integral: We write the fraction as a geometric series to obtaiñ dθ v F q cos(θ) iω n l f m cos(mθ) Here we strictly have to assume v F q < |ω n |, but we later see that we can extend the result through analyticity in ω n . We expand the cosine power using the binomial theorem and finally the θ integral is trivial. Finally we perform the l sum to obtain: This is manifestly analytic in ω n away from the imaginary axis so we may use this also for v F q > |ω n |.

F Real space resummed 1-loop boson polarization
Here we calculate the Fourier transform of the resummed bubbles of Eq. (88).
ρρ(x) bubbles = 1 λ 2 Here Π 2 (ω, k) is the 2-vertex fermion loop with density vertices, see Eq. (91). Instead of performing the geometric sum we do the Fourier transform term by term. Going to polar coordinates for the spatial integral and performing the angular integral we have where J 0 is the zeroth Bessel function of the first kind. Considering |τ | r −1 and x v F r −1 we are only interested in energies and momenta k rv −1 F and |ω| r and we can approximate the boson propagator as 1/r 2 . We can then perform the k integral where Finally we perform the ω integral and the sum in n ρρ(x) bubbles, where This decays as (x 2 + v 2 F τ 2 ) −3/2 at large x and τ .

G Calculating h ± with Landau-damping corrections
We want to calculate the h ± functions with a corrected boson propagator D LD (q) = 1 ω 2 + c 2 (q 2 x + q 2 y ) + r − Π 2 [λ(θ) 2 ](q) with Π 2 [λ(θ) 2 ](q) given by Eq. (91). For T = 0, 0 < r we can use the result of Appendix B since D is gapped and still decays at least algebraically for large energies and momenta. We proceed to finite temperatures. Eq. (76) still applies and we need only consider the n = 0 contribution. However, the boson self-energy is 0 for ω n = 0 so there is no difference from the non-Landau damped case apart from in the finite part. Finally we consider the QCP. Here we use a simplified propagator D IR (q) = 1 motivated by the low-energy scaling in Eq. (92) to find the large x limit (y = 0 so we neglect q x above) and verify the result using numerics and the full propagator. We find: at the QCP.