Non-Gaussian fixed points in fermionic field theories without auxiliary Bose-fields

The functional equation governing the renormalization flow of fermionic field theories is investigated in $d$ dimensions without introducing auxiliary Bose-fields on the example of the Gross-Neveu and the Nambu--Jona-Lasinio model. The UV safe fixed points and the eigenvectors of the renormalization group equations linearized around them are found in the local potential approximation. The results are compared carefully with those obtained with partial bosonisation. The results do not receive any correction in the next-to-leading order approximation of the gradient expansion of the effective action.


Introduction
The Functional Renormalization Group (FRG) has developed into an important investigation tool of the large distance behavior of strongly interacting quantum field theories [1][2][3][4][5][6]. In particular, the emergence of bound states/condensates can be studied very efficiently by introducing the corresponding composite fields with appropriate quantum numbers into the set of operators from which the low-energy effective action builds up [7,8]. Typically, these fields are introduced at short distances without the proper kinetic term through some appropriate "Dirac-δ" functionals equating composites formed from the "microscopic" fields with the "observable" fields. Formally, this construction associates with the new fields vanishing wavefunction renormalization constants Z composite (k = ) = 0 at the defining ultraviolet scale .
If there are dynamical objects in the corresponding channel, their wavefunction renormalization should grow away from zero when one reaches the compositeness scale. Below this momentum scale a local effective field, not revealing any a e-mail: jakovac@caesar.elte.hu b e-mail: patkos@galaxy.elte.hu c e-mail: posfay.peter@gmail.com internal structure should represent them, possessing its own kinetic term. This expectation was checked in the auxiliary field formulation of the O(N )-model [8], where the composite field introduced at the "microscopic" scale via Hubbard-Stratonovich transformation, that is, with no kinetic term, became at low scale a propagating dynamical degree of freedom on its own. Similarly, large anomalous scaling corrections were shown in the Yukawa coupling of boson-fermion models when searching for new non-Gaussian (interacting) fixed points [9].
Ultraviolet stable non-Gaussian fixed points provide promising alternative for consistent UV-completion of perturbatively non-renormalizable theories, like the quantized Einstein-Hilbert gravitational theory [10]. An asymptotic safety scenario could circumvent the triviality problem of the Higgs sector of the Standard Model [11,12]. Such ideas were put forward for the first time for the UV-completion of quantum electrodynamics [13][14][15]. Another actual issue of interest is to restrict the mass spectra of excitations in effective models of particle physics, like the Higgs sector of the Standard Model [16,17]. The study of analogous questions in simpler theories helps to develop appropriate methods of investigations.
A compelling example for the above scenario is represented by theories with four-fermion coupling like the Gross-Neveu [18] or the Nambu-Jona-Lasinio [19,20] model. One might introduce into these theories bosonic fields corresponding to certain fermionic bilinears through the delta functions δ( a... μ... (x) − (ψ(x) a... μ... ψ(x))), where a... μ... is a conveniently chosen matrix with a set of internal (a) and Lorentz (μ) indices. At the "microscopic" scale these fields do not have any dynamics. In successful searches for non-Gaussian fixed points substantial running of the wavefunction normalization of the composites has been observed and exploited [8,9].
In a recent paper we have proposed a scheme where one can explore the fixed point structure of fermionic models in the framework of the Wetterich equation without introducing auxiliary variables [21]. Such approach has been introduced earlier in connection with the dynamical breakdown of chiral symmetry in gauged Nambu-Jona-Lasinio models [22][23][24][25]. Fermionic evolution equations were developed for QCD by Meggiolaro and Wetterich [26] truncated at the four-Fermi level. Various three-dimensional theories have been investigated recently including into the effective action the full set of four-Fermi operators taking into account Fierz-relations among them [27,28]. There is continuous interest in fermionic theories also in condensed matter physics, where the expansion of the effective action in powers of the Fermi fields is usually also truncated at the level of four-Fermi interactions, sometimes including also momentum dependent (non-local) six-Fermi vertices [29]. Our paper constructed a rather general framework for the application of the Renormalization Group method to purely fermionic relativistic field theories without this limitation. For the fermionic "potential" of the Gross-Neveu and the chiral Nambu-Jona-Lasinio model fully explicit evolution equations were constructed. We focus in the present paper on mapping the fixed point pattern of these theories and compare our results with those obtained with approaches employing partial bosonization.
The running of the couplings starts slightly below the compositeness scale, where one can treat the composite objects discussed above as elementary (point-like) and introduce an effective potential depending on arbitrary powers of the invariants directly formed from the fermion background. At high enough momentum, above the compositeness scale higher powers of the invariants should get smeared into nonlocal combinations of the Fermi fields, reconciling in this way the Grassmannian nature of these variables with the arbitrary powers apparently present in the coarse grained potential.
Although there exist a considerable number of works in the literature dealing with the local fermionic potential approximation [21][22][23][24][25], and in our earlier publication [21] we have already worked out a rather general framework for the treatment of this formalism, it is worth to describe in a mathematically more accurate and less intuitive way how the Local Potential Approximation (LPA) is introduced in the fermionic case and what kind of approximations lie in the background.
First of all we have to pin down that the fermionic effective action contains an arbitrary power of the fermionic variables. Indeed, k [¯ , ] is the generator of the proper multifermion vertices at scale k: These proper vertices are not zero for any n. An example of the one-loop contribution in a theory with four-fermion vertices (like the Gross-Neveu or the Nambu-Jona-Lasinio model) can be seen in Fig. 1.
The assumption behind the LPA, both for bosons and fermions, is that, for the running of the effective action, the most important contributions come from that kinematical regime, where the vertex varies much faster in spacetime than the propagators. This means that any diagram contributing to the running which contains a non-local proper vertex (n) k can be approximated by the contribution where the vertex is concentrated to a single point, cf. Fig. 2. The point-like vertex assumption is good if the value of diagrams a and b on Fig. 2 are numerically close. Then we can replace the generator of the first diagram, which is the complete effective action, by the generator of the vertices of the second diagram, which is n>1 U n (ψ(x)ψ(x)) n . (Actually, also the dependence on other fermion bilinears, compatible with the Lorentz and the internal invariance of the theory, is allowed. The restriction to the scalar combination does not restrict the generality of our arguments.) When we seek a formal derivation, we make the approximation that the value of the proper vertex where L is the aforementioned compositeness, or localization length scale. This means that all the fields are localized effectively within a small neighborhood of x 1 , we call it V x 1 , its volume we denote by V . According to the LPA we assume that the fields are slowly varying on the scale L, then the proper vertex can be substituted by an average value¯ The factorization of the average value out of the non-local vertex must be done carefully, to avoid the appearance of a formal zero due to the fermionic nature of the variables. The key observation is that for any power n. Thus the approximate formula for the effective action which corresponds to the diagram on the right of Fig. 2 is obtained by first putting the neighboringψ and ψ fields to the same point, and then factoring out the average value of the proper vertex. So we can write As an effective way of writing we use the notation in the limit V → 0: This defines a (quasi) local potential for the fermionic fields. The physical conclusion is that instead of point-like fermion fields we should work with fermion bilinears averaged on patches, if we wish to study the point-like limit of higher order 2n-fermion couplings. The gradient expansion of the effective action is a series expansion in increasingly non-local terms. It represents a unique hierarchy only if the field content is not enlarged by introducing propagating composite fields. Our aim is to study the gradient expansion in terms of the original Fermi fields and compare with the results of a different truncation of the derivative expansion arising from the introduction of propagating auxiliary fields. Since the two cases have different kinetic parts in the Lagrangian, one might expect different convergence rates in the search for interacting UV-stable fixed point theories.
In the present paper we shall demonstrate that some relevant results demonstrated earlier in the auxiliary Bose field formulation in the next-to-leading order (NLO) of the gradient expansion of the Wetterich equation can be obtained also in pure fermionic LPA, provided we keep all powers in the fermionic potential. In particular, we study the effect of the anomalous dimension of the wavefunction normalization parameter of the fermions and find that it vanishes in the ground state, demonstrating this way the stability of our result at NLO of the gradient expansion.
The paper is organized as follows. In Sect. 2 we reformulate the version of the Wetterich equation derived in [21] for the Gross-Neveu model in a space-dependent fermionic background. It is projected on a constant background in Sect. 3, where the dependence of the effective potential on a non-zero scalar composite condensate is calculated. Here we present also the results of a fully analogous study of the Nambu-Jona-Lasinio model. In Sect. 4 it is shown that there is no anomalous scaling for fermionic fields, giving more robustness to the results obtained in LPA. In the Conclusions we compare our results with previous investigations. In particular, we show that the N f = ∞ results do reproduce all features of the d = 3 non-Gaussian asymptotically safe fixed point.

Wetterich equation for the N f -flavor Gross-Neveu model in inhomogeneous background
The Ansatz which corresponds to the next-to-leading order (NLO) of the gradient expansion of the Euclidean effective action , taking into account the scale dependence of the wavefunction renormalization of the defining Fermi fields is the following: where we have written out explicitly the bispinor index α and the flavor index l; the quantity I (and with it the quantum action k ) is invariant under the global discrete chiral symmetry transformation The operator content of the potential part of (5) is not complete. In principle all Lorentz and chiral invariant quartic combinations should have been included. The number of independent variables is then reduced by the Fierz-relations [27,28,30]. Such an Ansatz can be called Fierz-complete. It was established in [21] that in the complete expression of the right hand side of the Wetterich equation only the invariant I (x) appears, although in separate piecewise contributions also other invariants (namely, (ψ(x)γ m ψ(x)) 2 ) are generated. This observation is also true for the one-flavor Nambu-Jona-Lasinio model to be discussed below. In view of these findings our study, though not being Fierz-complete, is selfconsistent.
It is worth to discuss in some detail the physical range of variation of the invariant variable I and the characteristics of the potential U . This can be done explicitly for N f = ∞ with the help of the auxiliary field formulation. First, we analyze the case U = g 2 /(2N f )I . Its action is rewritten with the auxiliary field σ (x) as The model at N f = ∞ is solved by finding the saddle point of the effective action arising after integrating over the fermions [31]. A phase transition occurs into the broken symmetry phase at some g 2 cr < 0. For g 2 < g 2 cr the auxiliary field σ has a non-zero, real expectation value M, which determines also the size of the fermionic condensate: Since g 2 is negative, the auxiliary potential −N f ρ/g 2 > 0 can be interpreted as a physically stable potential of the σfield. This solution is matched with the mean-field potential of the original model by requiring Substituting the saddle point value of (ψψ) , one recognizes that I varies along the positive axis, and its potential energy is bounded from above in the broken symmetry phase. This argument is made more general by replacing in the auxiliary Then the saddle point equation is (ψψ) = σ U aux (ρ), and the matching of the N f = ∞ potentials leads to the meanfield relation This relation implies that a stable power-like asymptotic behavior −U aux ∼ a n ρ n , a n > 0 corresponds to an asymptotic behavior −(2n − 1)a n I n for U GN . It is natural to expect that the established characteristics of U (I ) is carried over also to the finite N f case. The derivation of the Wetterich equation for the effective action with x-dependent Fermi fields (ψ(x), ψ(x)) and their transposed doublers (ψ T (x),ψ T (x)) closely follows the steps presented in our previous publication [21]. We start with a form where the traces with respect to the bispinor and flavor indices have already been done: Here (ψ G kŨ ψ) stands forψ(x) α j G αβ, jl k (x, y)Ũ (y)ψ β l (y) and summation is understood over all discrete indices. A similar detailed expression corresponds to (ψ T G k , the flavor-diagonal, infrared regularized propagators, are given as and Here F k (γ m ∂ m ) is a non-local functional built with γ m ∂ m which freezes efficiently out the propagation modes with wave numbers below the actual normalization scale k. For its Fourier transform there are several propositions which will appear explicitly below. In these expressions one also introduces Below when discussing the scale dependence of Z k , we shall use the short hand notation

The local potential approximation and its fixed points
One projects the Wetterich equation on the local potential by substituting into its right hand side constant Fermi fields (ψ 0 , ψ 0 ). After performing the operations indicated on the right hand side of Eq. (10) one finds an expression which depends only on the scalar invariant (ψ 0 ψ 0 ) and was given explicitly with optimized infrared regulator [32] in Eq. (31) of [21]. This equation is transformed in two steps into a form convenient for finding the fixed points and the scaling exponents characterizing the behavior of different operators in its neighborhood. First, one introduces the following dimensionless rescaled variables taking into account the anomalous scaling of the wavefunction renormalization parameter ln Z k ∼ −η ln k: We search for a non-Gaussian fixed point solution of this equation in the LPA, where one sets η = 0. In order to make our treatment easier to follow, we consider a second rescaling related to the large-N f scaling of the different quantities: where Q d = S d /(d(2π) d ) contains the surface S d of the d-dimensional unit sphere. These two steps lead to the following evolution equation for y k (x): where the prime denotes the derivative with respect to x and ∂ t = ∂/∂(ln k).
Since these are the coefficients of the Taylor expansion of the fermionic potential which have physical significance, providing the point-like limit of the 2n-fermion vertices, it is adequate to search for the fixed point of this RG-equation in form of a power series: One finds the following equation for l n * : where the function F is a nonlinear expression of the coefficients with indices lower than n. It is easy to solve it after one finds the non-zero solution of the equation for l 1 * : This equation was already given in [21] for the N f = ∞ case and was shown to coincide with the result of [9] obtained with an Ansatz truncated at n = 1. It is worthwhile to emphasize that the non-Gaussian fixed point exists in the physical range of l 1 * only for N f > 1/2. The apparent singularities at N f = 1/2 inherited from the denominator of l 1 * in later formulas do not have any physical meaning. Using the value of l 1 * we can determine higher order coefficients, too. l 1 * = 0 yields l i * = 0 (the Gaussian fixed point), while in the non-Gaussian fixed point the l i * coefficients have non-trivial value. The explicit expression for l 2 * , for instance, reads .
The evolution of the non-Gaussian fixed point with changing dimensionality was analytically determined at N f = ∞ in [9]. It was found that the fixed point coordinates of the couplings λ * 2n were all proportional to d − 4, therefore the non-Gaussian fixed point merges with the Gaussian one when d → 4. When taking Eq. (19) at N f = ∞ one finds the following equation for l ∞ 2 * (l ∞ 1 * = −(d − 2)/4): This would require in d = 4 F(l ∞ 1 * ) = 0, which is not fulfilled. Therefore there is no non-Gaussian solution for d = 4 at N f = ∞, nor in our treatment.
Keeping, however, the ∼ 1/N f term in (19) for finite N f a non-Gaussian fixed point with finite "coordinates" persists up to d ≤ (8N f − 12)/(2N f − 5) as one can see from the denominator of (21). For finite N f the existence of an upper critical dimension was not discussed in [9]. Our fermionic LPA, analytic for arbitrary values of N f is not sufficient to settle this question. The momentum dependence of the four-Fermi coupling probably strongly influences the conclusion.
Going to the lower critical dimension one might remark that for d → 2 the non-Gaussian fixed point merges with the Gaussian one in agreement with [9].
One can reconstruct from the power series also the potential using the following procedure. The first two terms of the right hand side of Eq. (17) determine the large-x asymptotics of y * : y * as ∼ x d/ (2(d−1)) . In d = 3 the power equals 3/4 and one quickly can check the last two terms on the right hand side of (17) asymptotically vanish consistently. The Taylor series should sum asymptotically into this power law.
Being a polynomial, however, they cannot converge uniformly to x d/(2(d−1)) . This results in a wild oscillation of the different terms observed also in Ref. [9]. This behavior is illustrated in Fig. 3, where the power series of the potential are presented as obtained with n max = 5, 10, 15, 20 for N f = 2, d = 3, respectively.
One notes the uniform behavior of the series in a finite and symmetric neighborhood of the origin. This problem is cured by factoring out (1 + x 2 ) d/(4(d−1)) from the power series of the potential, which is insensitive to the sign of x. This factor has the correct asymptotics, while behaves polynomially for small x values. We know that the ratio of the original power series and of the asymptotic factor necessarily approaches a constant for large x. This can be achieved using symmetric Padé approximants to this ratio. Finally we obtain the fixed point potential with the following expression: Here the function Padé N N refers to the (N , N ) Padé approximant generated from the 2N th order Taylor series of the expression in the squared bracket.
Using the Padé approximation may be dangerous since the polynomials in its numerator and denominator may produce artificial zeroes and poles, respectively. However, after separating the correct asymptotics as described above, we gained the experience that most choices for N > 12 led to a smooth and uniformly converging sequence of potentials. The variation of the fixed point potential with N f is illustrated in Fig. 4, where the Padé approximants are displayed for various N f values together with the exact numerical solution of the N f = ∞ case, again for d = 3. The flow equation emerging from right hand side of (17) for N f = ∞ can be solved with the method of characteristics [33,34]. Its simple numerical implementation for the onevariable problem of the fixed point potential provides a test for the reliability of the method sketched above. The relative difference between the exact fixed point potential and the Padé approximants starts as negligible near the origin and saturates for x > 10 below 4 %. In Fig. 3 we also display the resulting fixed point potential together with the power series approximants. One recognizes the asymptotic series nature of the Taylor series which coincides on shorter and shorter interval with the exact solution.
Still, one can make use of this expansion in solving the linearized eigenvalue problem which by the structure of (19) is of a lower triangle matrix form. Each power represents there-fore an eigendirection with the following scaling exponents for the corresponding couplings: In d = 3 there is a single relevant direction (n = 1, GN 1 = 1), irrespective of the value of N f . This result is compatible with the N f = 8 Monte-Carlo simulation of [35], but it deviates from the behavior of the short series of 1/N fexpansion [35,36], which display increasing N f -dependence below N f ∼ 20. The exponent calculated with numerical solutions of the Wetterich equation as applied to the auxiliary field formulation stays rather close to the unit value for N f = 3, 4 [9,37]. The relevant eigenvalue of the N f = 1 case is significantly different [38].
The present spectra of exponents in the N f → ∞ limit exactly reproduces the result obtained in the partially bosonized representation [9]. The extra Yukawa coupling h k , is not running there (∂ t h k = 0) due to the non-trivial scaling exponent of the auxiliary field. For the remaining (irrelevant) operator set one can establish a clear correspondence by the scaling exponents between the 2n-fermion couplings l n of our treatment and the n-boson vertices λ 2n of the auxiliary field formulation. The values found from (24) approach the limiting N f = ∞ values more steeply than the corresponding exponents determined in the auxiliary field formulation [9]. By the correspondence the renormalization flow pattern in the coupling space l n is easily mapped onto the flow in the λ 2n -space of Ref. [9] around the non-Gaussian fixed point. This correspondence gives some support to treatψ 0 ψ 0 = ρ as a true bosonic field [25], but it is cleaner to think in terms of the "average" correspondence where V is the volume defined by the compositeness scale (cf. the discussion of the fermionic effective potential in the introduction). Before proceeding to the investigation of the effects which the wavefunction renormalization might exert on the above results we shortly summarize the results of a rather analogous LPA analysis performed in the N f = 1 Nambu-Jona-Lasinio model: One has to go through the same steps as for the GN-model, starting from Eq. (51) of Ref. [21]. The quantities I NJL and U NJL have the same canonical dimensions like the corresponding quantities of GN-model, furthermore one can scale out also the phase space factor Q d rather similarly, one can introduce the scaled variables: and arrive at the scaled RGE: After introducing the dimensionless variables, one finds with a Taylor-series search a non-Gaussian fixed point for the 2nfermion couplings. All higher function can be given in terms of the fixed point value of the n = 1 coupling The fixed point equation for l n * has the same structure as for the GN-model: The corresponding system of linearized flow equations for the deviation of the point-like limit of the 2n-fermion couplings from their fixed point values, denoted as δl n = l n −l n * has again triangular form One finds for the scaling exponents when using the value of l 1 * : The single relevant exponent 1 = d − 2 is the same as for the Gross-Neveu model, the remaining irrelevant part of the spectra is different. Here again one can compare our results with earlier investigations of the pure fermion representation, truncated at n = 1. The relevant scaling exponent agrees with the result found for d = 4 in [28] and also the fixed point value of l 1 * is the same if one takes into account our slightly different conventions in defining the spinor variables. As we mentioned before our investigation is not Fierzcomplete, only Fierz-consistent. The study of Ref. [28] was extended to include also vectorial (Thirring-type) interaction which led also to another non-trivial fixed point, though the number of relevant operators and the corresponding scaling exponent coincides with our finding when specified to d = 4. It has to be noted that for the non-Gaussian fixed point of the NJL-model we did not find any signal for an upper critical dimension.
Although we work with the pure fermionic theory, it is also possible to estimate the scaling of an indirectly defined effective Yukawa interaction. The Yukawa interaction "pulls apart" the point-like four-fermion interaction inserting a scalar propagator between pairs of the external fermion legs, therefore at zero momentum a relation exists between the coupling constants. Using the bosonic propagator at zero momentum, it reads h 2 /m 2 σ = 1 , where h is the Yukawa coupling, m σ is the effective boson mass. For the scaling of the Yukawa coupling we should take into account that where η σ is the anomalous dimension of the σ field). Then we find h 2 ∼ k − h 2 with h 2 = 1 + η σ − 2. This agrees with Eq. (40) of Ref. [9]. To access the bosonic anomalous dimension we need the bosonic wavefunction renormalization and use the scaling Z σ ∼ k −η σ . With the tentative assignment σ ∼ψψ, the bosonic dynamics should come from the insertion of the operator Z σ [∂ m (ψ(x)ψ(x))] 2 . This is a possible way of extending the fermionic treatment, but in the present formulation we do not have this operator; thus we have to set it zero for the bosonic anomalous dimension. Therefore now the conjectured scaling exponent for the Yukawa term is h 2 = 1 −2.
Since in the present model η ψ = 0 (cf. next section), we find

Wavefunction renormalization in the Gross-Neveu model
The projection of the Renormalization Group equation on the wave function renormalization constant is given by the equation The task is to substitute into the right hand side of this equation the three terms on the right hand side of the RGE (10). Diagrammatically, one has contributions from the set of Feynman-diagrams illustrated on Fig. 1, just two of the legs are not static. It is clear that there is a non-trivial dependence on the external momentum since these diagrams are overwhelmingly not tadpole-type. One observes that the result of the operations prescribed in (33) still depends on the background field. Generally one chooses its homogeneous value characterizing the ground state, that is, the minimum of the effective fixed point potential. This principle dictates us in the present case by the global features of U (I ) established in Sect. 2 to choose I 0 = (ψ 0 ψ 0 ) 2 = 0. The experience with various model investigations shows that the anomalous dimension η k = − ln Z k is proportional to the invariant of the theory, and therefore in the symmetric phase η = 0 [39]. Still one has to put I 0 = 0 only after carefully checking that it does not lead in the relevant integrals to infrared divergences, since the mass term in the propagators is proportional to this quantity.
In studies of pure fermionic formulation truncated at low powers of the invariants, the anomalous fermionic wavefunction exponent η ψ was found to vanish both in the GN-and the NJL-models [27,28]. With the non-truncated Ansatz the computation on the right hand side of (33) becomes quite tedious. Some of its details are worth to be presented, which follows below for the GN-model.
We start with TrlogG −1 k , and we promptly use the fact that G −1 k is diagonal in flavor. Its contribution can be expressed aŝ (we use an abbreviated notation for the functional derivative).
Since the regularized kinetic parts of g −1 and g (T )−1 do not depend on the fermion fields, a straightforward calculation gives for the different terms in the above integrands eventually evaluated on a constant ψ 0 background the following expressions: It is obvious that when substituting these expressions into the appropriate parts of (34) one encounters either γ α 1 α 2 m δ α 2 α 1 = 0 orψ α 1 0 γ α 1 α 2 m ψ α 2 0 . This latter is not included into the Ansatz (5), therefore we drop it also on the right hand side of the Wetterich equation. The same analysis goes through for Trlog G (T )−1 k , therefore even before setting I 0 to zero one recognizes that the first two terms of (10) do not contribute to the running of Z k .
For the evaluation of the contribution from the last term of (10) one can write an expression structurally identical to (34): After the tedious but straightforward computation of the derivatives one substitutes the constant spinorial background and exploits the fact that on such background the propagators g and g T are translationally invariant and also The Fourier transforms of the infrared regularized propagators g and g (T ) on a constant ψ 0 background read where r k F ( p) is the regularizing modification of the kinetic term and P F ( p 2 ) = p 2 (1 +r k F ( p)) 2 . Throughout the calculation we use the linear regulator [32]: The propagators determine also the Fourier transform of Q(x − y): Also here we omit all terms which would be proportional to the vectorial condensateψγ m ψ or give zero after the multiplication by γ α 1 α 2 (see above). When expanding the occurring integrals to linear order in the external momentum q one encounters expressions proportional to ∂ P F ( The coefficients of all integrals are proportional to some power of I 0 . The presence of ∂ P F ( p 2 )/∂ p 2 excludes the infrared region from the integration. In the integrals where the integrand is proportional to ∂r F ( p 2 )/∂ p 2 the p 2 -dependence of the rest of the integrands comes from the infrared regularized propagators. These terms, however, in the infrared region are frozen to constants, therefore the infrared contribution to the integral is usually of the form dpp d−1 ( pq) 2 / p 3 . It is regular for d > 2. Therefore in both types of integrals one can safely send the coefficients to zero.
Finally, there are also contributing integrals which are of the general form Since the term 1 + r F ( p) ∼ p −1 in the infrared region where the propagators are p-independent, these integrals are infrared regular for d + 1 > l, which is true for all occurring cases. Again, one is allowed to set in the coefficients of these integrals I 0 = 0. The whole rather tiresome discussion of some 15 integrals leads to the short conclusion that in the present formulation of the fermionic FRG:

Discussion
The existence of a non-Gaussian fixed point in the Gross-Neveu model in d = 3 with the single relevant (infrared repulsive) operator I = (ψψ) 2 around which theories with an infinite number of fermion flavors (N f = ∞) and with four-fermion coupling can be consistently renormalized has been established quite some time ago by investigating the ultraviolet behavior of the four-point function of the theory both with [31] and without [40,41] introducing an auxiliary field σ (x) ∼ψ(x)ψ(x). The present investigation confirmed the UV-safe behavior of the three-dimensional (d = 3) four-fermion models relying on the analysis of the functional renormalization group equations derived without introducing any auxiliary Bose fields. The spectra of scaling exponents appearing for the Gross-Neveu model in (24) for N f = ∞ fully reproduces the exponents found in [9] using the auxiliary formulation of the model. In the case of the latter approach the running of the auxiliary field renormalization (η σ = 4 − d) is essential. Our finding of zero anomalous dimension of the fermions even for finite N f is in qualitative agreement with the rather small η ψ found in the auxiliary field reformulation [9].
The existence of a non-Gaussian fixed point with a single infrared unstable direction appears for d > 2. Another important issue is to see if there is an upper dimension d max for the existence of the non-Gaussian fixed point. According to the analytic solution of the auxiliary field formulation of the ERG equations for N f = ∞ [9] the Gaussian and the non-Gaussian fixed points merge in d = 4. In our approach at N f = ∞ the fixed point value of the four-fermion coupling (l 2 * ) is pushed to infinity, which confirms that there is no consistent non-Gaussian fixed point in d = 4. For finite N f we find an N f -dependent value for d max > 4, clearly indicating the necessity to include momentum dependence into the effective action already on the four-fermion level.
In d = 3 we confirm the existence of a non-Gaussian fixed point with just one relevant eigendirection in the coupling space for all values of N f = 1, 2, .... The spectrum we find slightly deviates from those established numerically in [9]. The non-trivial influence from the anomalous dimensions of the Fermi field and, in particular, of the auxiliary field could be the reason for this difference. In particular this is so by introducing a kinetic term, which intuitively corresponds to the kinetic term of the auxiliary scalar field Z σ [∂ m (ψ(x)ψ(x))] 2 ; also in the present formulation effects related to the anomalous dimension of the composites might show up. These effects eventually should result in an improvement of the spectra of scaling exponents and the estimate for d max . By an intuitive correspondence it also provides an estimate for the scale dependence of the squared Yukawa coupling of the auxiliary field formulation.
The analogous results obtained for the Nambu-Jona-Lasinio model suggests that the existence of non-Gaussian fixed point(s) in d = 3 could be a generic feature of models with four-fermion invariants.
In conclusion, we find rather encouraging the level of agreement we found in analyzing the fixed point structure of the two model systems with and without bosonic composite fields. This fact hints at prospective efficient usage of the technique developed in [21] without the introduction of any scalar auxiliary field also in more complicated models.