Momentum dependence of quantum critical Dirac systems

We analyze fermionic criticality in relativistic 2+1 dimensional fermion systems using the functional renormalization group (FRG), concentrating on the Gross-Neveu (chiral Ising) and the Thirring model. While a variety of methods, including the FRG, appear to reach quantitative consensus for the critical regime of the Gross-Neveu model, the situation seems more diverse for the Thirring model with different methods yielding vastly different results. We present a first exploratory FRG study of such fermion systems including momentum-dependent couplings using pseudo-spectral methods. Our results corroborate the stability of results in Gross-Neveu-type universality classes, but indicate that momentum dependencies become more important in Thirring-type models for small flavor numbers. For larger flavor numbers, we confirm the existence of a non-Gaussian fixed point and thus a physical continuum limit. In the large-$N$ limit, we obtain an analytic solution for the momentum dependence of the fixed-point vertex.


I. INTRODUCTION
Surprisingly many condensed-matter systems support gapless Dirac-fermionic modes as emergent low-energy degrees of freedom, exhibiting a quasi-relativistic dispersion relation. This plethora of systems includes, for example, materials of graphene type, 3 He, surface states of topological insulators, or d-wave cuprate superconductors. Despite being rather different microscopically, these systems are summarized as the class of Dirac materials, see e.g. Refs. 1 and 2 for reviews.
The corresponding investigation of low-energy effective theories for such systems and associated phase transitions, e.g., semimetal-insulator transitions, has been a focus of theoretical research in recent years  . These theories are considered to transcend the standard Landau-Ginzburg-Wilson paradigm for critical phenomena, as gapless fermions remain present at the phase transition 36,40,41 possibly being subject to further spectator symmetries 42 .
The present work is devoted to a study of two wellknown models of relativistic Dirac fermions in 2+1 dimensional spacetime: the (irreducible) Gross-Neveu model 43 and the Thirring model 44 . Both models are classically defined in terms of a local four-fermion interaction, and exhibit a U(2N f ) symmetry, where N f denotes the number of massless (reducible 4-component) Dirac flavors. Both models are naively "perturbatively nonrenormalizable", but have shown to be nonperturbatively renormalizable, e.g. to any order in a 1/N f expansion [45][46][47][48][49] , and thus are simple examples 50 of the concept of asymptotic safety 51 . Despite their similarity, the long-range behavior can be rather different.
By contrast, the Thirring model does not exhibit a phase transition at large N f in agreement with results in the strict N f → ∞ limit 72,73 . Also, the order parameter is a bifundamental scalar in flavor space 74 , describing the breaking to a U(N f )⊗U(N f ) subgroup in case of gap formation. Early studies using a variety of methods suggested the existence of a critical flavor number N f,crit , above which no gapped phase exists, and below which a transition to a gapped phase occurs at sufficiently strong coupling. Quantitative estimates of N f,crit since then have spanned a wide range of values 72,[75][76][77][78][79][80][81][82][83] . As analyzed in Ref. 74, various reasons may have led to this diversity, e.g., also the presence/absence of exact chiral symmetry in the different methods. Recent lattice simulations with exact chiral symmetry 69,[84][85][86] indicate that there may be no symmetry breaking for any integer N f ≥ 2 at all. Whereas simulations using domain wall fermions provide evidence for 0 < N f,crit < 2 with some indications for N f,crit > 1 84 , the latest results of Ref. 86 using SLAC fermions and an analytic continuation of the parity-even theory to arbitrary real N f show clear evi-dence that the critical flavor is below unity, N f,crit < 1.
In view of these similarities and differences between the two models, we wish to address three immediate questions in this work: (i) Is there a way to quantify the higher degree of complexity of the Thirring model compared with the Gross-Neveu model? (ii) If so, can this information be used to construct better approximation schemes for the Thirring model? (iii) Is the Thirring model a UV complete (asymptotically safe) quantum field theory, despite the absence of a phase transition for a wide range of N f values? While the second question aims at an improvement of analytical methods, the answer to the third question is also relevant for the existence of a continuum limit of discrete lattice realizations.
For addressing these questions, we need a theoretical tool that can simultaneously deal with both models in the continuum, preferably for all values of N f and directly in d = 2 + 1 dimensions. These criteria are met by the functional renormalization group (FRG) which offers a substantial degree of flexibility of devising systematic nonperturbative approximation schemes. In fact, the FRG as formulated in terms of the Wetterich equation 87 has already widely been used for the two models under consideration. Common approximation schemes focus on local fermionic interactions 31,42,82,88 , or use partial-bosonization techniques together with a derivative expansion to also account for the emerging composite bosonic degrees of freedom 14,17,24,40,50,65,68,74,[89][90][91] . Whereas these methods work very well for the Gross-Neveu model, exhibiting apparent convergence for the quantitative estimates of critical exponents 68 , the Thirring model already on this level is more involved, requiring composite scalar and vector boson fields 74 as well as dynamical bosonization techniques [92][93][94] .
In the present work, we go beyond the limitations of these FRG approximation schemes by investigating the momentum dependence of the fermionic self-interactions. In this manner, we can specifically address the question as to whether the assumptions underlying the use of the derivative expansion are justified. For the present exploratory study, we concentrate specifically on the momentum dependence arising in the s channel, using modern pseudo-spectral methods as developed in Refs. 90 and 95 to resolve this dependence to a high accuracy.
Our results provide at least partial affirmative answers to all three questions raised above: (i) The quantitative relevance of momentum dependence provides for a measure to assess the higher degree of complexity of the Thirring model. By contrast, our results for the Gross-Neveu model lend strong support to the quantitative use of the derivative expansion. (ii) More adapted approximation schemes for the Thirring model hence have to resolve the more complex momentum dependence of the fermionic self-interactions in particular for smaller flavor numbers; for future work, a combination with dynamical bosonization techniques absorbing momentum dependencies in the composite channels appears recommended.
(iii) Our results provide further evidence for the existence of a non-Gaussian Thirring fixed point for a wide range of flavor numbers. The 2+1 dimensional Thirring model hence can be a UV-complete quantum field theory, also implying the existence of a continuum limit of lattice formulations independently of the observation of a specific chiral-symmetry breaking phase transition.
This work is structured as follows: in Section II we shortly discuss the relevant aspects of the FRG. Section III discusses the general structure of the renormalization group equations, whereas in Section IV we present the explicit flow equations. After reviewing the complementary bosonization approach and considering some limiting cases in Sections V and VI, we solve the flow equations for particular values of the flavor number in Section VII. We finish with a discussion of the results and a conclusion in Section VIII.

II. FUNCTIONAL RENORMALIZATION GROUP FOR FERMION SYSTEMS
We consider theories with N f flavors of relativistic fermions described by the Grassmann-valued spinor fields ψ i ,ψ i (i = 1, . . . , N f ) in d = 3 Euclidean dimensions. We also define a collective field variable Ψ := (ψ, ψ T ) T in the spirit of Nambu-Gorkov spinors. For the spin base, we use γ matrices, γ µ ∈ C dγ ×dγ , generating the Clifford algebra {γ µ , γ ν } = 2δ µν . We begin with a reducible representation with d γ = 4. A transition to an irreducible representation is described below.
For the following nonperturbative investigation, we use the functional renormalization group (FRG) which facilitates to systematically integrate out quantum fluctuations at all energy scales. We specifically make use of the FRG to compute the effective action Γ of a theory starting from its classical, microscopic action S. The key ingredient is a scale-dependent effective average action Γ k that only includes fluctuations with momenta p 2 k 2 . Thereby, Γ k interpolates between the microscopic action Γ Λ S for some high-energy scale k = Λ and the full effective action Γ k=0 = Γ. For this, we add a scale-dependent regulator term in momentum space to the classical action, which effectively decouples modes with momenta p 2 k 2 from the theory. We use symmetry preserving regulators of the form R k (p) = − / p r(p 2 /k 2 ) with a dimensionless regulator shape function r that diverges sufficiently fast for small argument to provide for a k-dependent gap of the fermionic propagator at low momenta. Note that representing an exact functional differential equation for Γ k that describes the flow from the microscopic action to the full effective action. Here, Γ k denotes the second functional derivative with respect to the fields Ψ, and the supertrace STr is a shorthand for a summation over discrete and integration over continuous variables.
While Eq. (2) is an exact equation, exact solutions are generally difficult to obtain. In practice, it is useful to work with systematic truncation schemes, concentrating on an accessible subset of theory space. There are two common choices to truncate the theory space. A derivative expansion spans the action in terms of operators sorted according to the number of derivatives acting on the quantum fields; the advantage is that full field dependencies can be retained. This approximation scheme is especially useful for investigating the action of bosonic order-parameter fields in the study of spontaneous symmetry breaking and the description of condensation effects. By contrast, a vertex expansion spans the action in terms of correlation functions with an increasing number of external legs. This allows to address the full momentum dependence of the correlations, having become a successful strategy to explore gauge theories and gravity where the couplings are inherently momentum dependent.
For fermionic systems as considered here, previous studies have concentrated on either lowestorder derivative expansions in a purely fermionic description 42,82,88 , or derivative expansions in a partially bosonized description within the so-called LPA' approximation 14,17,24,40,50,74 , including multi-bosonic interactions 65 ; for some models, even complete next-toleading order results are available 68,91 .
In this work, we use a vertex expansion to study the momentum dependence of four-fermion couplings for the first time for the present set of models. More precisely, we account for the full momentum dependence of the propagator, a projected choice of momentum dependencies of the four-fermion vertices, and neglect possible contributions from higher vertices.
Resolving momentum dependent propagators and the flow of higher-order vertex functions has become a central tool for analyzing flow equations for the theory of the strong interactions 108-111 as well as quantum gravity [112][113][114] . Within the FRG, also approximation strategies to resolve both momentum and field dependence are available 115,116 . The present work is among the few examples, where not only the full propagator, but also the functional dependence of vertex functions on a momentum channel are resolved near an interacting fixed point.

III. FLOW OF MOMENTUM-DEPENDENT COUPLINGS
We consider models defined on the basis of the flavor and Dirac structure of four-fermion interactions, i.e. models with local interactions of the form where p 1 ≡ −p 2 − p 2 − p 4 in the Fourier transformed version in the second line. Our convention assumes all momenta as in-coming into the vertex. The couplingḡ parametrizes the strength of the interaction and O 1 and O 2 are flavor and Dirac space matrices characterizing the vertex. In the following, we mostly consider the Gross-Neveu 43 and Thirring 44 models defined in terms of O GN = 1 ⊗ 1 and O Th = 1 ⊗ γ µ , respectively, for both O 1 and O 2 . Here, the first factor of the tensorial structure refers to the flavor index and the second to Dirac space.
Upon integrating out fluctuations, the corresponding 4-point correlators will acquire nonlocal contributions. These can be parametrized by promoting the coupling in (3) to a function of the external momenta, g =ḡ(p 1 , p 2 , p 3 , p 4 ). Because of momentum conservation, only three of these vectors are independent. Lorentz symmetry reduces the number of variables to a maximum of six linearly independent invariants that can be constructed from three independent momentum vectors. It is still useful to stick to the redundant notation above, as it allows us to keep track of the flow of loop momenta in the Wetterich equation (2) and goes along well with the pictorial language of Feynman diagrams.
Since we are eventually interested in fixed points and critical behavior of specific models, we switch to dimensionless quantities by measuring all dimensionful quantities in units of the RG scale k. E.g. we use dimensionless momentap = p/k, and also introduce renormalized fieldsψ(p) := Z 1/2 (p 2 )ψ(p),ψ(p) := Z 1/2 (p 2 )ψ(p), where Z(p 2 ) denotes a momentum-dependent wave function renormalization.
Similarly, the dimensionless, renormalized four-fermion couplings parametrizing the 4-point correlator are defined asg( This rescaling leads to additional contributions to the flow equation of the couplings,  where the ellipsis stands for the quantum fluctuations, i.e. the loop-induced flow ofḡ. Here, we have also introduced a generalized, momentum-dependent anomalous dimension η(p 2 ) := −∂ t ln Z(p 2 ). In the following, the dimensionless notation is implicitly understood, so that we drop all tildes and use the bar notation to distinguish dimensionful couplings from their dimensionless counterparts.
The Wetterich equation (2) has a one-loop structure. A general analysis of generic four-fermion vertices shows that there are four different processes or diagrams contributing to the RG flow of the coupling functions. These are depicted in Fig. 1, where the process in panel (a) exists in two variants, see below. The diagrams are characterized by different momentum transfers across the loop and are conveniently labeled in terms of the Mandelstam variables 117 s = ( In this spirit, we refer to the diagrams in Fig. 1 (a)-(c) as s, t, and u channel processes. We emphasize that the diagrams denote processes involving full propagators and fully-momentum dependent vertices and hence are nonperturbative.
In order to cast them into a more explicit form, we introduce the threshold kernel a1 a2 K(p, q; η) := which encodes the propagation of the loop fermions. Here q denotes the loop momentum, p the momentum transfer, η is the anomalous dimension function, and the regulator shape function r was introduced below Eq. (1). By convention, we always choose q to be the momentum flowing through the regulator insertion. The prefactors a 1 and a 2 parametrize different momentum modes of the propagating loop fermion. The threshold kernel is related to the threshold function  Fig. 1 then are in correspondence with the diagram functionals a2 K(p 3 +p 4 , q; η) g 1 (p 1 , p 2 , q , −q)g 2 (p 3 , −q , q, p 4 ) + (q ←→ −q ) q =q+p3+p4 , (7b) a1 a2 B t [g 1 , g 2 ; η](p i ) := 1 N q a1 a2 K(p 2 +p 3 , q; η) g 1 (p 1 , q , −q, p 4 )g 2 (−q , p 2 , p 3 , q) + (q ←→ −q ) q =q+p2+p3 , (7c) a1 a2 B u [g 1 , g 2 ; η](p i ) := 1 N q a1 a2 K(p 2 +p 4 , q; η) g 1 (p 1 , q , p 3 , −q)g 2 (−q , p 2 , q, p 4 ) + (q ←→ −q ) q =q+p2+p4 , whereÑ = N f d γ = 4N f counts the independent spinor components. The notation "+ (q ←→ −q )" here implies that all preceding terms within the same bracket level have to be added with q and −q interchanged. In relation to the diagrams, we identify the left/top vertices with the coupling g 1 and the right/bottom vertices with the coupling g 2 . The two distinct versions of the s-channel diagram mentioned above differ by an exchange of the 1and 3-legs in the second vertex.
In a similar fashion, we can analyze the self-energy contributions from the four-fermion couplings entering the flow of the wave function renormalization. The corresponding diagram, Fig. 2, exists in two variants again. The associated diagram functionals are These definitions will allow us to express the flows of wave function renormalizations and coupling functions of arbitrary four-fermion models to fourth order in the vertex expansion.

IV. FLOW EQUATIONS OF THE GROSS-NEVEU-THIRRING MODEL
After these general considerations, we now study a class of models defined by the symmetries of the Thirring model. The Thirring model is defined by the microscopic action 44 with the conventionÑ = N f d γ = 4N f , as introduced below Eq. (7). The vector-type Thirring interaction (ψγ µ ψ) 2 exhibits obvious invariances under U(N f ) flavor rotations as well as continuous transformations ψ → e iαM± ψ,ψ →ψe ±iαM± generated by the Dirac algebra elements M + ∈ {γ 4 , γ 5 } and M − ∈ {1, iγ 4 γ 5 } (forming a U(2) group). Because of the partial similarity to chiral or axial transformations in d = 4, the latter symmetry is often called "chiral" in the literature. In fact, these two symmetries are actually subgroups of a larger U(2N f ) symmetry, which will become more obvious in the irreducible spinor formulation introduced below. In addition, the model also exhibits discrete C, P , and T symmetries which can be formulated in various ways 74 . Up to Fierz transformations, there is only one more linearly independent local four-fermion interaction sharing the symmetries of the Thirring model 42,74,82,89 . In our conventions, this can be written as (ψγ 45 ψ) 2 with γ 45 := iγ 4 γ 5 .
The U(2N f ) symmetry of the model becomes explicit by decomposing the N f four-component Dirac spinors ψ,ψ into 2N f two-component Weyl spinors χ,χ transforming in the irreducible representation of the Clifford algebra 74,82,89 . Extracting χ,χ as chiral projections of ψ,ψ by means of the projectors P ± = 1 2 (1 ± γ 45 ), the fermion bilinears are mapped intō with σ µ denoting the Pauli matrices, The interactions discussed above hence correspond to the Thirring interaction and the Gross-Neveu interaction 43 of the irreducible representation, featuring an explicit U(2N f ) symmetry. Even if the bare action was defined by the Thirring interaction (9) at some high scale, fluctuations would generically generate also the Gross-Neveu interactions. Therefore, we include both interaction structures in our ansatz for the effective average action, which reads in momentum space: (12) Inserting Γ k into the Wetterich equation (2) and ignoring contributions from n > 4-point correlators, we can extract the flow equations or β functions for the dimensionless couplings, (13b) j=1 η(p 2 j ). In addition, we suppressed the dependence on the anomalous dimension function η in the diagram functionals. We observe that all four types of diagrams worked out in Sect. III enter the flow equations.
However, we expect that the s-channel dominates for several reasons: For the flow of g GN , there are more schannel diagrams than t-or u-channel ones. Also, the Th ] come with opposite signs, such that mutual cancellations can occur. For the flow of g Th , this is the case for the diagrams of order g 2 GN as well as the ones of order g 2 Th . Moreover, the s2-diagrams are the only ones remaining in the limit N f → ∞. Our expectation of s-channel dominance is self-consistently verified below.
The flow of the wave function renormalization is governed by the anomalous dimension function which can similarly be extracted from the Wetterich equation, Eqs. (13) and (14) form a system of coupled integrodifferential equations for the coupling and anomalous dimension functions to which, unfortunately, no exact so-lution is known to us. In view of the expected s-channel dominance, we can reduce the complexity substantially, if we assume the coupling functions to depend on s exclusively, g(p i ) = g(s = [p 3 + p 4 ] 2 ) for both g GN and g Th . This s-channel approximation allows us to compute fixed points and critical exponents in Sect. VII.

V. PARTIAL BOSONIZATION
Let us first establish the connection between our momentum-dependent fermionic formulation and partial bosonization techniques used in previous studies 17,24,40,50,65,68,101 . Partial bosonization, in fact, allows one to recover parts of the momentum dependence of four-fermion vertices in a setting with pointlike Yukawa interactions 92,101,118 . Moreover, it gives direct access to the chiral condensate ψ ψ , making it particularly useful to address questions of dynamical symmetry breaking 74 . Along the lines of a Hubbard-Stratonovich transformation 119,120 , we rewrite the Gross-Neveu-Thirring action (12) with the aid of scalar and vector boson fields φ and V µ . An ansatz for the action with local interactions in coordinate space reads in momentum space: Indeed, for vanishing kinetic terms of the boson fields (Z φ = Z V =Ā V = 0), φ and V µ become purely auxiliary Hubbard-Stratonovich fields: they can be integrated out immediately, and give back the pointlike limit of Eq.
As the kinetic terms are generated inevitably by the RG flow, the action (15) describes momentum transfer between pairs of fermions as mediated by the bosons. Thus, the mixed boson-fermion theory parametrizes an effective momentum-dependent four-fermion interaction. In the above form, our ansatz for the bosonized effective action corresponds to an improved local potential approximation (LPA'), where only the lowest order of the effective potential, namely the mass termsm 2 φ φ 2 /2 and m 2 V V 2 /2, is kept. Of course, it is straightforward to include higher truncations of the effective potential 50,74 or higher-order derivatives 68 , but the present simple form suffices for the following argument.
In order to assess the momentum dependence encoded in the bosonized action, we can map the action (15) back onto a purely fermionic action. As Eq. (15) is quadratic in the boson fields, the use of the equations of motion for the boson fields is exact also on the level of the functional integral. Solving the classical equations of motion δΓ k /δφ = 0 and δΓ k /δV µ = 0, we find The first line illustrates, for instance, that the vacuum expectation value of φ is linked to the condensate χχ = ψ γ 45 ψ .
Substituting Eqs. (16) back into the action (15), we obtain the corresponding purely fermionic action, where p s = p 3 + p 4 . A comparison to the action (12) of the purely fermionic model reveals that the bosonization procedure induces the Gross-Neveu and Thirring interactions with an effective s-channel momentum dependence. Moreover, the vector boson parametrizes a further interaction vertex resembling the square of the kinetic term, (χ / p s χ) 2 . This type of interaction is present ifĀ V = Z V and is not resolved in our ansatz (12) for the momentumdependent Gross-Neveu-Thirring model. We come back to this observation below.

VI. LIMITING CASES
Let us go back to the purely fermionic case described by the action (12) and analyze the simplifying limits of pointlike interactions and infinite flavor number, respectively.

A. Pointlike limit
The RG flow of the Gross-Neveu-Thirring model in the limit of momentum-independent (pointlike) couplings has been investigated in detail 31,42,74,82,88,89 . This limit forms the basis of derivative expansions using bosonization techniques to restore parts of the momentum dependence 17,24,40,50,68,74,89,101 . In this limit, the self-energy diagrams (8) do not acquire a nontrivial momentum dependence, and hence the anomalous dimension vanishes. Using (6) the coupling flows reduce to in accordance, e.g., with Refs. 42 and 82. We show the resulting phase diagram in Fig. 3 for N f = 4 and using the Litim regulator 121,122 . As frequently noted in the literature, the pure Gross-Neveu model forms an invariant subspace in this limit: Setting g Th = 0, we can solve Eq. (18a) independently of and consistently with Eq. (18b). This fact is used in many studies of the Gross-Neveu model to restrict the theory space to just the Gross-Neveu vertex, even though, in principle, also the Thirring vertex has the same symmetries as the irreducible Gross-Neveu model; in the case of the reducible Gross-Neveu model, there are even three additional vertices compatible with the reduced symmetries of the model 42 . With regard to the flow equations (13) of the momentum-dependent model, we already observe that this property no longer holds true in the general momentum-dependent case, because there are two diagrams of order g 2 GN contributing to the flow of g Th . We will come back to this point in Sect. VII.
The phase diagram with its fixed points and flows towards the long-range physics of the pointlike model has already been discussed extensively in Ref. 82. For our purposes, we merely recall that the phase diagram exhibits four fixed points as shown in Fig. 3. In addition to the trivial Gaussian fixed point, there are two critical fixed points with one relevant eigenperturbation, which we will refer to as the Gross axis nonetheless. In addition, there is one more fixed point having two relevant perturbations and named B in the following.
B. Large-N f limit As a second instructive limit, we study the fixed-point structure of the momentum-dependent RG flow of the vertices in the limit N f → ∞, employing the s-channel approximation. We observe that the diagram functionals (7) and (8) scale with 1/N f and hence can be neglected in this limit. In fact, only the s2 diagrams in the flow equations for g GN and g Th remain finite. Furthermore, the anomalous dimension (14) vanishes exactly because Σ 2 [g; η] ≡ 0 in the s-channel approximation. As a result, the flow equations (13a) and (13b) decouple, and both take the form with a 1 = 1, a 2 = 0 for g GN and a 1 = −1, a 2 = 1 for g Th .
Here, we have defined being the integrated threshold kernel (5) for s-channel processes. Note that this quantity depends only on the Mandelstam variable s = (p 3 + p 4 ) 2 . The substantial technical simplification of the large-N f limit arises from the fact that the integral no longer involves the coupling function g. Therefore, the associated fixed-point equation (∂ t g * = 0) is an ordinary differential equation. As its solution is a full function rather than a "fixed point", we refer to such solutions as fixed functions in the following. In the present case, the general solution is , g ∞ = const. (21) Here g ∞ represents an a priori arbitrary constant, implying that we obtain a one-parameter family of potential fixed point solutions. The condition that g * (s) should not acquire singularities for s > 0 is fulfilled for values g ∞ > 0 (g ∞ < 0) for g GN * (g Th * ), because 1 0 K ( −1 1 K) is strictly positive (negative). In fact, even the value g ∞ = 0 leads to a well-defined coupling function for all s ≥ 0, but its asymptotic large-s behavior changes qualitatively: we find g * (s) ∼ s if g ∞ = 0, whereas g * (s) ∼ 1/ √ s if g ∞ > 0 (g ∞ < 0). For the Litim regulator 121,122 the remaining integral in Eq. (21) can even be evaluated analytically, yielding a solution in closed form. Since the resulting expression is rather extensive, we simply plot the solutions in Fig. 4. The arrangement of axes is similar to Fig. 3, with the projected third axis exhibiting the s dependence. The B fixed function is obtained by combining the nontrivial solutions for g GN * and g Th * .
The integration constant g ∞ can be fixed by including an additional piece of information from finite-N f results: for this, we first note that the solution function (21) is not analytic: in a power series expansion around s = 0, we also encounter half-integer powers of s. This peculiarity is caused by the kernel a1 a2 K(s) regardless of the choice of regulator as is discussed in Appendix A. It turns out that the lowest-order nonanalyticity ∼ √ s is absent for finite-N f solutions. This motivates to impose a "maximum regularity condition" also on the large-N f solution: we require the constant g ∞ to be fixed such that as many derivatives of g * in s = 0 exist as possible. This regularity condition implies g ∞ = 1 16 for g GN * and g ∞ = − 1 32 for g Th * , which corresponds also to the values used for the plots in Fig. 4.
The resulting eigenperturbations are for the Gaussian and nontrivial fixed functions, respectively, with C being a normalization constant in both cases. In order for the perturbations to be well-defined in s = 0, the spectrum must be bounded from above by Re θ ≤ −1 for the Gaussian and Re θ ≤ 1 for the nontrivial fixed points. If we additionally demand for maximum regularity, the critical exponents become quantized with spacing 2, i.e. θ = −1, −3, −5, . . . for perturbations around Gaussian and θ = 1, −1, −3, . . . for perturbations around nontrivial solutions. In fact, the importance of a proper choice of the Hilbert space for spanning the fixed functions is also well studied for the analysis of fixed potentials in field space for Wilson-Fisher type fixed points 100,123,124 and beyond [125][126][127][128][129] . We observe a similar feature here for fixed functions in momentum space in the form of the maximum regularity criterion.
Since we can perturb in the g GN and g Th directions of theory space independently, the resulting spectrum of any of the four fixed functions is obtained by merging one spectrum for either direction, leaving the Gaussian fixed point with zero relevant exponents, the GN and Th fixed points with one relevant exponent θ 1 = 1 each, and the B fixed point with two relevant exponents θ 1 = θ 2 = 1. Incidentally, this is also the spectrum found in the partially bosonized language 50,74 as well as a generalized fermionic truncation scheme 32 .

VII. FIXED FUNCTIONS AND SPECTRA FOR FINITE FLAVOR NUMBERS
With the information from the bosonized version and two limiting cases of the fermionic formulation, we are now equipped to study the momentum-dependent system for finite N f . The situation is considerably more complex because all four types of diagram functionals contribute to the flow equations (13), and all involve integrals of the unknown functions g GN , g Th , and η. In view of the complications posed by these nonlinear integro-differential equations, we continue to rely on the s-channel approximation, and resort to numerical methods to compute fixed functions in theory space and their associated spectra.

A. Solution strategy
We first consider the fixed point equations, i.e. ∂ t g GN = ∂ t g Th = 0 in Eq. (13).
The integrodifferential structure of the equation and the nonanalyticity of the threshold kernel constitute rather challenging premises for any numerical computation. Notwithstanding, the pseudo-spectral method using Chebyshev rational functions 130-132 has proved to be a versatile and robust tool to obtain a solution of arbitrary accuracy, at least in principle. As an example from the FRG context, it has been used successfully to obtain high-accuracy solutions of effective potentials for various model systems before 68,90,91,95,[133][134][135][136] . We will therefore expand the couplings and the anomalous dimension function as a series of Chebyshev rational functions R n (n = 0, 1, . . .). These are obtained from the Chebyshev polynomials of the first kind T n by a compactification of the semi-infinite interval [0, ∞), The parameter L is arbitrary, but fixed. It should typically reflect the order of the length scales of the function to be modeled 131 ; for our purposes, L = 5 turned out to yield satisfying convergence. By definition, the Chebyshev rational functions inherit useful properties such as analyticity, boundedness, and orthogonality from their polynomial parents. Most importantly, they form a basis set of L 2 ([0, ∞)) with measure dx L/x/(x+L), and the Chebyshev series is guaranteed to converge on the entire domain if the expanded function is square-integrable 137 . This should be contrasted, for instance, with Taylor series, which converge only within a disc up to the closest singularity in the complex plane. Following this approach, we write the coupling functions g GN , g Th and the anomalous dimension function η as truncated series of Chebyshev rational functions, (25) Note that we parametrize the argument for the coupling functions as √ s rather than s. This is motivated by the asymptotic behavior as s → ∞: If the fixed point coupling functions are bounded, the quantum terms in Eq. (13) are suppressed by at least 1 √ s compared to the scaling terms due to the asymptotics of the threshold kernel (5). Similarly, the anomalous dimension function (14) vanishes in the limit p 2 → ∞. Consequently, we find that g * (s) ∼ 1 √ s as s → ∞ for both g GN and g Th . The ansatz (25) is therefore better adapted to the asymptotic behavior than a plain s parametrization. This turns out to be crucial for achieving convergence of the numerical solutions. Moreover, we already observed in the previous section that the threshold kernel generates half-integer powers of s, so that an expansion of g * of the above form appears natural. By contrast, the asymptotic behavior of the anomalous dimension η * , following from Eq. (14) together with the asymptotics of the couplings, is given by η * (p 2 ) ∼ 1 p 2 for p 2 → ∞. This justifies the parametrization of Eq. (25).
The remaining task is to determine the expansion coefficients g GN,n , g Th,n , and η n using the pseudo-spectral method 131 . Plugging the series expansions (25) into the relations (13) and (14) and demanding the equations to hold on a discrete set of collocations points {s m } n g max m=0 and {p 2 m } n η max m=0 leaves us with a nonlinear system of algebraic equations for precisely these coefficients. To solve this system, we use a Newton-Raphson method 138,139 to iteratively reduce the residuals starting from some initial guess for the coefficients. For the collocation points, we choose the roots of the next-order expansion function R nmax+1 . This ensures the pseudo-spectral method to coincide with the spectral method if the inner products in the latter are evaluated numerically using an optimal Gaussian quadrature rule 131 .
For the concrete computation, there is one subtlety which remains to be addressed: Even after a restriction to a pure s-channel dependence, the flow equations (13) for the coupling functions do not uniquely determine the remaining s dependence of the couplings. For a given value of s, the directions of the external momenta in three dimensional spacetime are not uniquely fixed. In order to compute the corresponding integrals, we have to make a choice for these directions, and the result for the s-dependent coupling does generically depend on this choice. Our strategy now is to perform the computa-tions for a set of distinct geometric configurations of momenta, monitoring the variation of the resulting solution functions upon changing the configuration. In fact, this serves as a self-consistency check for the s-channel approximation: if the different configurations show large variations, the contributions from the t and u channels must be sizable.
The set of momentum geometries consists of the four different momentum configurations shown in Fig. 5, which are expected to display possible variations rather prominently. In each of them, the magnitude of all four momenta is the same, p 2 i = p 2 = s 2+2 cos (p3,p4) . In the parallel configuration in panel (a), we have t = u = 0, leading to the "purest" s-dependence. The orthogonal configurations in panels (b) and (c) accentuate one of the t-or u-channels against the other: In the t-favored configuration (b), we have t = s and u = 0, whereas in the u-favored one (c), t = 0 and u = s. In the symmetric configuration in panel (d), the angle between any two momentum vectors is the same, i.e. cos (p i , p j ) = − 1 In order to compute the spectra of fixed point solutions, we proceed similarly to the large-N f case discussed above by perturbing around a fixed function solution with g GN (s) = g GN, * (s) + e −θt ε GN (s) and g Th (s) = g Th, * (s)+e −θt ε Th (s). Contrary to the large-N f case, however, we cannot consider perturbations in the g GN and g Th directions independently because the equations do not decouple. Expanding to first order in ε, we arrive at the linearized flow equations In order to solve these equations and extract the eigenvalues θ, we again use a pseudo-spectral approach and expand the perturbations in terms of Chebyshev rational functions. It is instructive to first examine the asymptotic behavior of the eigenperturbations for s → ∞. The diagram functionals are again all suppressed by at least the decay of the threshold kernel in this limit, so that only the scaling terms remain. Consequently, we find ε(s) ∼ s a , where the asymptotic power a and the critical exponent θ are related by the asymptotic scaling relation This implies that a and θ balance each other: The more irrelevant the perturbation ε (smaller θ), the faster it grows as s → ∞ (larger a). Hence, in order to probe the irrelevant part of the spectra beyond Re θ = −1, we have to allow for asymptotically growing eigenperturbations.
Since the Chebyshev rational functions are bounded, we multiply by a polynomial in s, motivating the ansatz ε(s) = (s + L) amax n ε max n=0 ε n R n (s) .
For the choice of the value of a max , we followed two different approaches. For the first one, we set a max to a constant numerical value, thus keeping the maximum growth of the perturbations fixed. The advantage of this method is that it reduces Eqs. (26) to an algebraic eigenvalue problem once the equations are evaluated on a collocation grid. As a disadvantage, the required expansion order n ε max grows as we want to reach deeper into the irrelevant part of the spectrum. Since the possible large-s scaling ε(s) ∼ s a of the ansatz (28) ranges from a = a max − n ε max to a = a max , it can only capture perturbations with exponents θ such that Re θ ∈ [−2a max − 1, 2n ε max − 2a max − 1] according to Eq. (27). In order to retain the relevant part of the spectrum, we therefore have to increase n ε max when increasing a max . As the strongly irrelevant exponents are not of particular interest anyway, this disadvantage is typically not an issue.
Our second approach uses a scheme with a running asymptotics, where a max = −1−θ 2 is adapted to the eigenvalue θ according to the scaling relation (27). This method is computationally much more demanding because the eigenvalue equations (26) become nonlinear. In practice, it turns out that both methods yield identical exponents within the validity limits of our approximations obtained from varying the momentum configuration and/or the regularization scheme.
For the actual computations of fixed functions, eigenperturbations and critical exponents, we have implemented the procedures described above in a C++ program. The multidimensional numerical integrations were carried out using the Cuhre algorithm of the CUBA library [140][141][142][143] , which is a deterministic, high-precision integration scheme featuring globally adaptive subdivisions. As for the regularization, we mostly used an exponential regulator with shape function Algebraic eigenvalue equations were solved using the Eigen3 library 144 .

B. Fixed point structure
The overall fixed point structure resembles the pointlike and large-N f limits, but there are a few modifications that deserve attention. We show a selection of fixed functions for N f = 12 (upper row) and N f = 2 (lower row) in Fig. 6. The actual solutions for the four momentum configurations of Fig. 5 are all plotted as solid lines. To visualize the deviation arising from the variation of momentum configurations, we additionally display elliptic tubes around the mean of these configurations with principle axes given by the standard deviations. As in the limiting cases discussed above, we can distinguish four different solutions, the endpoint of which for s → 0 approximately coincides with the fixed-point positions in the pointlike limit. Correspondingly, we associate the fixed functions with the Gaussian, the (irreducible) Gross-Neveu, the Thirring, and the "B" universality classes. The coupling and anomalous dimension functions decay for s → ∞ or p 2 → ∞, respectively. For relatively large flavor numbers, e.g. N f = 12 in Fig. 6a-b, the variation with momentum configuration is very small for all fixed points and the solutions are close to the infinite-N f solution (cf. Fig. 4). The s → 0 limit of the fixed function is also quantitatively close to the fixed-point value in the pointlike limit for the Gross-Neveu universality class, whereas the Thirring and the "B" fixed points show some larger deviations. The agreement with the large-N f limit is a manifestation of the fact that the s-channel approximation is closed in this limit. For smaller values of N f , the GN fixed function remains remarkably stable when varying the momentum configuration, whereas the deviations increase significantly for the Th and B fixed points. In fact, we cannot even make a firm statement about the existence of the B fixed point in the current approximation, because convergence of the Newton-Raphson iteration slowed down significantly here when N f < 2. In the symmetric configuration, we have even been unable to find a solution.
A plot of the leading eigenvalues of the perturbations around the three nontrivial fixed functions is presented in Fig. 7. For each eigenvalue, we show its minimum and maximum value for different momentum configurations. These values essentially coincide for the GN fixed point, and they vary only moderately and for small N f at the Th fixed point. Both of them have one relevant exponent larger than zero as in the pointlike and large-N f limits. The B fixed point shows stronger variations and apparently also a qualitative change of behavior: a third exponent appears to become relevant in the u-favored orthogonal configuration for N f = 2, 3, casting doubt on the reliability of our results in this parameter region. In any case, the large-N f limit is rediscovered consistently, as the eigenvalues approach the values found in Section VI B.
We will now analyze all three nontrivial fixed points in more detail.

C. Gross-Neveu fixed point
In the pointlike approximation, we found that the Gross-Neveu vertex forms an invariant subspace of theory space: the coupling g GN does not induce the flow of any of the other symmetry-compatible vertices 42 . In particular, no flow of g Th in Eq. (18b) is induced on the Gross-Neveu axis. The flow equations of the fully momentumdependent model (13) tell us that this does no longer hold in general. There are t-and u-channel processes of order g 2 GN entering the flow of g Th , which simply happen to cancel in the pointlike limit. Consequently, the Gross-Neveu coupling does not parametrize an invariant subspace any longer, because we cannot find a consistent general fixed function g GN, * after setting g Th ≡ 0 in Eqs. (18). This should serve as a word of caution when restricting the theory space to invariant subspaces in the pointlike limit. Nevertheless, the violation of this invariance is reassuringly mild in the Gross-Neveu universality class, at least in the s-channel approximation. For instance in Fig. 6, there is no visible difference. To illustrate how far the fixed functions actually leave the Gross-Neveu subspace, we compute the maximum value of g Th * and relate it to g GN * at the same s value. This ratio, g Th * (smax) gGN * (smax) with s max such that |g Th * (s)| ≤ |g Th * (s max )| for all s, is plotted in Fig. 8. It reaches up to 4% for N f = 1 and decays like 1/N f with the flavor number. Hence the restriction to g GN alone remains a satisfactory approximation for all integer N f .
In fact, our approach, which concentrates on resolving the momentum structure of fermionic four-point vertices, should not be expected to be adequate to produce highprecision estimates for critical exponents. The latter certainly require an accurate treatment of order parameter fluctuations. From large-N f computations, it is clear that local order parameter self-interactions contribute to the dynamics near criticality. In the fermionic description, these would correspond to 8-fermion operators which are not included in the present study. For instance, large-N f expansions predict that the relevant eigenvalue θ 1 approaches the large-N f limit from below, θ 1 1 − 0.270 N f , for N f → ∞, whereas our results exhibit an approach to θ 1 | N f →∞ = 1 from above. Incidentally, the large-N f result θ 1 = 1 can be proved to hold for any critical fixed point in the pointlike limit 42,147 . We observe that this correct limit is also obtained with momentum dependent vertices.
In order to assess the relevance of the present momentum-resolving approximation, we take a closer look at the critical exponents in the regime of small flavor number N f , as displayed in Fig. 7a. For instance in the case N f = 1, we find θ 1 1.23 (3), η ψ 0.028 (7), where the uncertainties reflect the variation of the result upon choosing different regulators. This can be compared to most recent four-loop 4 − expansion (incl. Borel resummation) 71 , yielding θ 1 1.114 (33), η ψ 0.102 (12), or similarly to FRG results including partial bosonization: θ 1 1.075(4), η ψ 0.0645. We observe an agreement on the 10. . . 15% level for θ 1 , but a substantially larger deviation for η ψ .
The picture is similar for the graphene-like case N f = 2, where we obtain θ 1 1.09(1), η ψ 0.012 (2), to be compared to the four-loop 4− expansion (incl. Borel resummation), θ 1 0.993 (27), η ψ 0.043 (12), and the FRG with partial bosonization, θ 1 0.994(2), η ψ 0.0276. These results suggest that the class of operators considered in this work do have a satisfactory overlap with those that dominate the leading critical exponent, whereas the anomalous dimension is not well captured. Incidentally, various results from Monte-Carlo simulations and conformal bootstrap show larger deviations (already among each other, as well as with respect to the quoted expansion and FRG findings).

D. Thirring fixed point
As for the Th fixed point solution, the variation with momentum configuration is much larger and becomes significant for small flavor numbers. Fig. 6 demonstrates that the s-channel approximation becomes rather crude towards smaller flavor numbers and cannot capture all essential dependencies. As expected, the outermost Th solutions in Fig. 6 correspond to the two orthogonal configurations where the influence of either of the remaining two Mandelstam channels is maximally exposed against the other; the parallel and symmetric configurations lie in between. The same holds true for the anomalous dimensions, where we observe considerable variation between the configurations for small momenta. From top to bottom, the curves in both plots belong to the t-favored orthogonal, symmetric, parallel, and u-favored orthogonal configurations, respectively. We conclude that the momentum dependence becomes increasingly important as N f becomes smaller, indicating that neither the pointlike limit nor the s-channel approximation can be expected to resolve all features. Nevertheless, the qualitative picture is consistent, and the additional momentum dependence does not modify the topology of the phase diagram. In particular, our results confirm the existence of the Thirring fixed point and thus of a nonperturbatively renormalizable continuum quantum field theory for a wide range of N f values. This viewpoint is corroborated by our results for the critical exponents displayed in Fig. 7b. Deviations between the momentum configurations become visible only for small N f and are stronger in the sub-leading exponents, whose accuracy is naturally lower. The leading exponent, however, varies only little down to N f = 2, where the differences affect the second decimal place and are presumably still smaller than those from regulator variations for fixed momentum configuration. Exponents for various flavor numbers are listed in Tab. I with error ranges corresponding to the configurational variation bands in Fig. 6. For the anomalous dimension η ψ = η * (0), the differences between the momentum configurations are larger. The general trend is again a decreasing value with increasing N f , eventually converging to η ψ = 0 as N f → ∞. The anomalous dimension for N f = 1, however, cannot reliably be extracted since η ψ flips sign for some momentum configurations in our approximation.
Critical exponents of the Thirring model have also been computed using the FRG with dynamical bosonization for scalar as well as vector boson channels 74 . The values found this way for the relevant exponent θ 1 are generally smaller than our results, e.g., θ 1 = 0.4 for N f = 2 or θ 1 = 0.92 for N f = 4, but the differences reduce as N f is increased. Also both FRG results agree on the finding that the large-N f limit θ 1 = 1 is approached from below. In view of our analysis in Section V, a better resolution of the vector channel suggests to include the kinetic term squared, (χ / p s χ) 2 , in future studies. This vertex is naturally entailed in the dynamically bosonized description 74 at the same level as the Thirring interaction.
While this analysis of the fixed-point properties provides information about the high-energy behavior of the model and its possible UV completion as an asymptotically safe quantum field theory, many studies in the literature focus on the long-range IR behavior and possible low-energy phases. Here, a variety of earlier studies suggested the existence of a critical flavor number N f,crit below which spontaneous breaking of the chiral/flavor symmetry can occur at sufficiently strong coupling [72][73][74][75][76][77][80][81][82][83] . As discussed in Refs. 42 and 74, the concrete realization of chiral symmetry in these calculations may have strongly influenced the critical quantities, potentially leading to a different universality class. For instance, the phase transitions observed with lattice realizations using staggered fermions [78][79][80]83,148,149 , exhibited chiral phase transitions for smaller flavor numbers, e.g., for N f = 2 reporting a critical exponent θ 1 in the range 1.2 . . . 1.4. Only recently, lattice simulations have been performed that realize the chiral U(2N f ) symmetry using domain-wall 69,84 or SLAC 85,86 fermions. The results using domain-wall fermions indicate that 0 < N f,crit < 2 with some preference for N f,crit > 1 by the data 84 . Latest results of Ref. 86 using SLAC fermions with an analytic continuation of the parity-even theory to arbitrary real N f provide clear evidence that the critical flavor number is smaller than unity, N f,crit < 1. The latter implies that no phase transition exists for integer N f .
The present exploratory study should not be expected to give accurate estimates for the critical flavor number. The observation of an increasing relevance of momentum dependence at small N f indicates that a larger operator basis is needed to extract the critical quantities reliably. Nevertheless, we can use the range of uncertainties of the present calculation to estimate a range of possible values for N f,crit . For this, we use the observation made in Refs. 74 and 82 that the Th fixed point moves from the Thirring subspace towards an NJL-type subspace as N f decreases. This suggests to interpret the presence or absence of chiral symmetry breaking as a competition between the NJL and Thirring interactions, where the former supports but the latter suppresses the formation of a chiral condensate (as suggested by the large-N f limit 72,73 ). In the irreducible representation, the NJL TABLE I. Critical exponents with error bounds given by the variation with the external momentum configuration. This provides a lower bound of the error only and does not include deviations between different regularization schemes. In addition, the differences between the distances to the NJL and Thirring subspaces at s = 0 and s = ∞ are given.  (1) 0.59(1) 0.27 (25) vertex is given by the scalar flavor nonsinglet 74 By a Fierz transformation, we can relate the coupling function g NJL (p i ) to the Gross-Neveu and Thirring couplings and obtain with the above specified convention for the external momenta. This implies that the NJL subspace is characterized by g GN = g Th . In the pointlike limit, the fixed-point coordinates of the Thirring fixed point fulfill this condition for N f = 1.75 74 . As a measure for the distance of a point (g GN , g Th ) in theory space to the NJL subspace, we use which is a function of the external momenta or, in our case, the Mandelstam variable s. In a similar fashion, we can define the distance from the Thirring subspace as d Th := |g GN |. We assess the competition between the Thirring and NJL interactions by comparing these distances in Fig. 9. More precisely, we plot the difference d NJL − d Th normalized by the location of the fixed point as a function of s: As before, we indicate the variation between the different momentum configurations by a shaded area. Negative values of ∆ NJL−Th mean that the solution is closer to the NJL subspace, whereas for positive values it is closer to the Thirring subspace. The dashed lines in the figure show the values of this difference as found in the pointlike limit. In this case, the value N f = 5.01 where d NJL = d Th The s-channel momentum dependence reveals that the fixed point solutions tend to be closer to the Thirring subspace for small values of s. This suggests N f,crit to be smaller than the naive pointlike estimate and thus also of the more elaborate result based on the derivative expansion 74 . Extracting an improved estimate is inhibited by the variations of ∆ NJL−Th with the momentum configuration which increase substantially as s → ∞. As ∆ NJL−Th is a relative measure, this increase can be attributed in part to the fact that g GN , g Th → 0 as s → ∞. Explicit estimates and variations for ∆ NJL−Th at s = 0 and s = ∞ are also listed in Tab. I. Still, we can use d NJL d Th , or ∆ NJL−Th 0, as an indicator for chiral symmetry breaking to occur. Requiring this condition for some s and all momentum configurations would lead to 3 < N f,crit < 4. If we demanded that it holds at s = 0 for some momentum configuration, we would find 2 < N f,crit < 3. Requiring it for all s and all momentum configurations gives N f,crit < 2. The result of this last and most strict criterion is compatible with those of recent lattice simulations 84,86 .
In the literature, it has been conjectured that the chiral-symmetry breaking phase transition of the 2+1 dimensional Thirring model and that of QED 3 could be related 76,148 or even in the same universality class 150 . If so, our result based on the most strict criterion is also compatible with recent lattice simulations 151,152 of QED 3 .

E. B fixed point
Contrary to the GN and Th fixed points, convergence of the fixed point couplings worsens significantly for small flavor numbers at the B fixed point. This manifests itself in two ways: On the one hand, we needed to fine-tune the initial guesses provided to the fixed point solver much more thoroughly to achieve convergence of the Newton-Raphson iteration. On the other hand, the residuals of the obtained solution functions between the collocation points increase strongly at fixed order of the Chebyshev expansions when N f is lowered. These trends are absent at the GN and Th fixed points.
At the same time, we notice increasing deviations between the different momentum configurations, similar to the Th fixed point. Furthermore, the B fixed point solutions appear to merge into the Th solutions for small flavor numbers as s → ∞ (see, e.g., Fig. 6c). In some configurations, the Gross-Neveu coupling g GN crosses over from g GN > 0 to g GN < 0 as s → ∞, experiencing a root at some value of s. However, there is no clear pattern regarding the positions of these roots. From asymptotic considerations (cf. Section VII A), we know that all four fixed points eventually converge to g GN = g Th = 0 as s → ∞. The observations here naively suggest that there may be a fusion of the Th and B fixed points for large flavor number and momentum transfer. With the increase of residuals for small N f in mind, however, this merger is likely to be an artifact indicating that the solver gets stuck at two different fixed points for small and large values of s.
Hence, while the existence of the B fixed point is wellestablished for large flavor numbers, it may either vanish and wander off to infinity, or diversify, potentially developing a continuum of solutions for small N f . Clarifying the situation will again require a finer resolution of the coupling functions' momentum dependence.
Likewise, the spectrum at the B fixed point (Fig. 7c) shows greater uncertainties compared to the Th (let alone the GN) fixed point. The red shading for N f = 1...2 is to indicate that the large residuals and the lack of convergence in the symmetric configuration cast doubts on whether the obtained functions are actual solutions. In any case, the fixed point has two relevant directions for large N f as in the pointlike limit and the bosonized formulation 74 . At N f = 2, 3, an additional third relevant exponent occurs in the u-favored orthogonal configuration. This could hint at further fixed points because the corresponding eigenflows must eventually approach an IR limit. Again, though, the extra fixed point might also lie at infinity. Nevertheless, since for these small flavor numbers the residuals of the corresponding fixed point solutions begin to increase, it remains unclear whether this extra relevant direction is physical.

VIII. CONCLUSIONS
We have studied a class of effective theories of Dirac materials in 2+1 dimensional spacetime, featuring N f gapless (reducible) Dirac fermion flavors on the classical level with a U(2N f ) flavor ("chiral") symmetry group. This class contains the (irreducible) Gross-Neveu model as well as the Thirring model classically defined in terms of local four-fermion interactions. These are widely studied similar models with partly rather different long-range properties.
We have focused on the fixed-point structure of these theories, exploring the momentum dependence of the full propagator and the interaction vertices for the first time. For this, we have used the FRG-framework and spanned the momentum dependence in terms of Mandelstam variables. In many regimes, a reduction to the s channel turns out to be self-consistent for quantitatively addressing the fixed-point properties. Using pseudo-spectral methods, we have constructed the momentum-dependent fixed-point solutions for the couplings at the different fixed-points of the model class. These fixed points reflect the various universality classes of this set of models which serve to define UV-complete interacting quantum field theories (asymptotically safe theories). More specifically, our results confirm the existence of the Gross-Neveu universality class for all flavors as well as the Thirring universality class for a wide range of flavors. In the large-N f limit, we are able to construct the momentum-dependent fixed-point solutions analytically. Despite the nontrivial momentum dependence which indicates the presence of an infinite set of operators at the fixed point, both models feature only one RG-relevant direction, making the fixed points a candidate for being associated to a 2ndorder quantum phase transition potentially visible in the long-range properties.
In fact, for the Gross-Neveu model, the fixed point can directly be linked to a chiral phase transition to a gapped phase with an Ising-type order parameter reflecting the Z 2 -symmetry breaking by a fermion condensation for all values of N f . The model is, however, not in the Ising universality class, as gapless fermionic modes contribute near the phase transition. Our momentum-dependent results exhibit a remarkable degree of stability of the Gross-Neveu-type interactions, lending strong support to previous quantitative FRG analyses based on derivative expansions and focusing on the Gross-Neveu channel.
The Thirring model has been suggested to be connected to a chiral phase transition of NJL-type with symmetry breaking pattern U(2N f ) → U(N f ) ⊗ U (N f ) which, however, does not occur at large N f . Whereas our momentum-dependent results also show a large degree of stability for larger flavor number, instabilities do appear towards smaller N f . This indicates that the small N f region requires a more careful resolution of the momentum structure of the interactions than previously anticipated. In particular, a determination of the critical flavor number N f,crit below which a strong-coupling quantum phase transition can occur is likely to be affected by this more complex structure of the Thirring model. Our explorative studies of the momentum structure in this regime indicate that N f,crit is smaller than previously anticipated, potentially so small that chiral symmetry does not break spontaneously in this model for any integer value of N f . While a precise estimate remains difficult, our results are compatible with recent lattice simulations that account for exact chiral symmetry providing evidence for N f,crit being near 84 or, in fact, below 86 N f = 1.
In addition to the free theory, the model class also contains another interacting fixed point (fixed point B) having two relevant directions for large N f . This fixed point shows the largest degree of instability with respect to a variation of momentum configurations. Hence, its properties and even its existence become questionable towards smaller N f , leaving us with the Gross-Neveu and Thirring universality classes as the relevant cases in this model set of maximal flavor symmetry.